變異函數(shù)模型參數(shù)的加權(quán)總體最小二乘回歸法

趙英文 王樂洋，2，3

1 東華理工大學(xué)測繪工程學(xué)院，南昌市廣蘭大道418號，330013

2 流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點實驗室，南昌市廣蘭大道418號，330013

3 江西省數(shù)字國土重點實驗室，南昌市廣蘭大道418號，330013

變異函數(shù)能同時描述區(qū)域變量的隨機性和結(jié)構(gòu)性［1］，在研究數(shù)據(jù)的分布特性和數(shù)據(jù)插值等方面發(fā)揮著重要的作用。變異函數(shù)模型參數(shù)估計方法中，人工擬合法采用肉眼來確定變異函數(shù)模型的參數(shù)，效率低且可靠性不高。文獻［2］提出用加權(quán)回歸多項式法來擬合變異函數(shù)模型的參數(shù)，但沒有解決參數(shù)正負(fù)號問題。線性規(guī)劃法［3－4］、最小二乘法［5－6］和加權(quán)最小二乘法［7－9］雖然提高了參數(shù)的計算效率，但只認(rèn)為變異函數(shù)值含有誤差，沒有考慮到變異函數(shù)模型中距離值的隨機誤差。文獻［10］首次提出總體最小二乘概念，從數(shù)值分析的角度解決系數(shù)矩陣也含有誤差的平差問題。當(dāng)系數(shù)矩陣含有誤差時，最小二乘解是有偏的，而總體最小二乘解是無偏的［11－12］。考慮到距離值的隨機誤差，文獻［13］提出用總體最小二乘法求解變異函數(shù)模型參數(shù)，認(rèn)為變異函數(shù)值和距離值是等精度的，把變異函數(shù)值的權(quán)陣作為行尺度矩陣左乘觀測向量和系數(shù)矩陣，然后用SVD 分解法進行解算。

本文將加權(quán)總體最小二乘回歸法引入到變異函數(shù)模型參數(shù)估計中。以冪函數(shù)模型為例，通過協(xié)方差傳播律［14］發(fā)現(xiàn)分組后變異函數(shù)值和距離值是不等精度的，并給出距離值的定權(quán)方法。再結(jié)合熵權(quán)法和點對數(shù)法進行參數(shù)的迭代求解，最后通過模擬數(shù)據(jù)和實測數(shù)據(jù)驗證加權(quán)總體最小二乘回歸法的合理性和優(yōu)越性。

1 變異函數(shù)模型及參數(shù)估計

1.1 變異函數(shù)模型

離散型變異函數(shù)γ（h）可由式（1）求解：

式中，h為兩點間的分隔距離，N（h）為距離為h的點對數(shù)目，Z（xi，yi）和Z（（xi，yi）＋h）分別為位置（xi，yi）和（xi，yi）＋h處的區(qū)域變化量。

根據(jù)不同的形狀和結(jié)構(gòu)，變異函數(shù)模型可分為球狀模型、高斯模型、指數(shù)模型、冪函數(shù)模型和孔穴效應(yīng)模型等［1］。冪函數(shù)模型為：

式中，M為常系數(shù)，α為冪指數(shù)。由文獻［1］知，α必須小于2；若α≥2，則冪函數(shù)不再是變異函數(shù)。

1.2 變異函數(shù)模型參數(shù)估計

變異函數(shù)模型事先是未知的，一般的做法是利用已知樣本點計算出的距離值h和變異函數(shù)值γ（h）擬合出相應(yīng)的變異函數(shù)模型。實際中點位分布散亂，需在計算出任意點對間的距離后對距離進行分組。給定一個容許誤差T，落在（h－T，h＋T）內(nèi)的距離值個數(shù)為N（h），則該組距離值為所有距離值的均值，即。

以冪函數(shù)模型為例，對式（2）兩邊取對數(shù)：

令

式中，y∈Rm×1，A∈Rm×n，x∈Rn×1，m為分組后距離的個數(shù)，n為未知參數(shù)個數(shù)，冪函數(shù)模型n取2。

考慮到觀測值的誤差，冪函數(shù)模型線性化后的回歸模型簡記為：

式中，y表示變異函數(shù)值，A為包含距離的系數(shù)陣，x為冪函數(shù)模型參數(shù)。

式中，e為向量y的誤差向量，EA為系數(shù)矩陣A的誤差矩陣，為參數(shù)最佳估值，vec（·）為矩陣按列拉直運算，為單位權(quán)方差，QE為A的協(xié)因數(shù)陣，Qe為y的協(xié)因數(shù)陣。

如果向量y和系數(shù)矩陣A的元素等精度，根據(jù)總體最小二乘平差準(zhǔn)則求得參數(shù)的總體最小二乘解：

式中，e為向量y的誤差向量，EA為系數(shù)矩陣A的誤差矩陣，為參數(shù)最佳估值，vec（·）為矩陣按列拉直運算，為單位權(quán)方差，QE為A的協(xié)因數(shù)陣，Qe為y的協(xié)因數(shù)陣。

如果向量y和系數(shù)矩陣A的元素等精度，根據(jù)總體最小二乘平差準(zhǔn)則求得參數(shù)的總體最小二乘解：

總體最小二乘解法主要有SVD 分解法和迭代解法［12］。顧及到系數(shù)矩陣A中含有常數(shù)列，文獻［15］提出一種顧及系數(shù)矩陣部分含有誤差的加權(quán)總體最小二乘迭代解法，把算法中的協(xié)因數(shù)陣變?yōu)閱挝魂嚕玫皆撈讲顔栴}的總體最小二乘解：

如果顧及y和A中元素不同的精度和貢獻程度，變異函數(shù)模型參數(shù)估計問題成為加權(quán)總體最小二乘平差問題。文獻［16］提出一種加權(quán)總體最小二乘算法，根據(jù)加權(quán)總體最小二乘平差準(zhǔn)則：

該算法的迭代過程為［16］：

1）計算加權(quán)最小二乘解作為迭代初始值：

3）迭代終止，得到解xi。

1.3 權(quán)的確定

1.3.1 向量y對應(yīng)權(quán)陣Py的確定

Py的確定方法主要有熵權(quán)法［7］和點對數(shù)法［9］。依據(jù)文獻［7］，把分組后的距離h和其對應(yīng)的個數(shù)N（h）作為影響因子構(gòu)造評價矩陣R，?。剑踙N（h）］，得到變異函數(shù)值的權(quán)Wi。文獻［9］提出的點對數(shù)法定權(quán)公式為：

1.3.2 系數(shù)矩陣A對應(yīng)權(quán)陣PA的確定

假設(shè)觀測值z的隨機誤差，令zi和zj分別等于式（1）中的Z（xi，yi）和Z（（xi，yi）＋h），由一個點對求得變異函數(shù)值γ（h）的方差D1：

顧及到式（1），由N（h）個點對求得變異函數(shù)值γ（h）的方差D（γ（h））：

分組后距離h由N（h）個符合分組條件的hij先求和再除以N（h）得到，則距離h的方差為：

由式（15）和（17），且一般情況下m個距離h對應(yīng)的N（h）不完全相同，所以分組后變異函數(shù)值和距離值是不等精度的。從而，距離的權(quán)可定義為方差的倒數(shù)，即，W（h）為m維列向量。

在計 算 中，取Py＝diag（Wi），PA＝diag

上述提到的權(quán)適用于冪函數(shù)模型，而對于其它的變異函數(shù)模型，可在模型線性化后，按照文獻［7］和［9］的方法確定權(quán)陣Py，依據(jù)本文的推導(dǎo)思路確定權(quán)陣PA。因為不同的函數(shù)模型對應(yīng)的系數(shù)矩陣會含有不同的距離表現(xiàn)形式，所以權(quán)陣PA需依據(jù)具體的函數(shù)模型確定。

2 算例分析

2.1 模擬算例

模擬一組規(guī)則分布的坐標(biāo)（x，y），并計算任意兩點間的距離。假設(shè)變異函數(shù)模型為冪函數(shù)模型，通過給定的冪函數(shù)γ（h）＝0.2h1.5代入h計算得到γ（h）。由MATLAB 產(chǎn)生一組均值為0、標(biāo)準(zhǔn)差初始值為0.009、步長為0.005、終止值為0.049的正態(tài)分布隨機誤差序列，根據(jù)式（15）和（17）分別在h和γ（h）加上相應(yīng)的誤差，共產(chǎn)生9組不同標(biāo)準(zhǔn)差下的數(shù)據(jù)，每組模擬200次。采用4種方法進行參數(shù)估計，其中WLS1 和WTLS1的定權(quán)方法是熵權(quán)法，WLS2 和WTLS2 的定權(quán)方法是點對數(shù)法，TLS采用文獻［15］方法，WTLS采用文獻［16］方法，‖ΔX‖為參數(shù)M和α的估值與真值之差的2范數(shù)。估計結(jié)果見表1。

表1 不同方法的估計結(jié)果Tab.1 The results estimated by different methods

表1中結(jié)果為9種不同標(biāo)準(zhǔn)差下的均值。由表1可知，加權(quán)總體最小二乘法所得的參數(shù)殘差范數(shù)均小于其他3種方法，采用點對數(shù)法定權(quán)的結(jié)果優(yōu)于熵權(quán)法定權(quán)的結(jié)果。由圖1可知，當(dāng)所加誤差較小時，幾種方法所得結(jié)果比較接近。隨著所加誤差的增大，加權(quán)總體最小二乘法所得結(jié)果更接近于參數(shù)的真值，說明加權(quán)總體最小二乘法在參數(shù)估計方面具有更高的精度和合理性。

2.2 高程異常插值算例

區(qū)域一數(shù)據(jù)來自于文獻［17］，其GPS控制網(wǎng)由17個同精度GPS水準(zhǔn)點構(gòu)成，高程異常值變化平緩。選取5個點作為已知點，剩余12個作為檢核點。區(qū)域二數(shù)據(jù)來自于文獻［18］，其GPS控制網(wǎng)由24個同精度GPS水準(zhǔn)點構(gòu)成，高程異常值變化較大。選取10 個點作為已知點，剩余14個作為檢核點。點位分布如圖2。

圖2 點位分布圖（左圖為區(qū)域一，右圖為區(qū)域二）Fig.2 Distribution of points（area 1on the left side，area 2on the right side）

對上述兩組數(shù)據(jù)分別計算距離值和變異函數(shù)值，在距離分組后選取冪函數(shù)模型作為變異函數(shù)模型，參數(shù)估計方法同上，使用RMS 和STD 作為評價指標(biāo)。

RMS為均方根預(yù)報誤差：

STD 為預(yù)報殘差標(biāo)準(zhǔn)差：

式中，Zi、Zi′分別為k個檢核點的觀測值和通過Kriging插值計算得到的檢核點估計值，ri是檢核點預(yù)報殘差，是檢核點預(yù)報殘差的均值。結(jié)果如表2、表3所示。

表2 區(qū)域一不同方法的檢核點預(yù)報結(jié)果Tab.2 The predicted results of check points calculated by different methods in area 1

表3 區(qū)域二不同方法的檢核點預(yù)報結(jié)果Tab.3 The predicted results of check points calculated by different methods in area 2

圖3 區(qū)域一檢核點預(yù)報值與觀測值差值曲線圖Fig.3 Difference of check points between the predicted values and the observed ones in area 1

圖4 區(qū)域二檢核點預(yù)報值與觀測值差值曲線圖Fig.4 Difference of check points between the predicted values and the observed ones in area 2

由表2、表3可知，無論RMS還是STD 都表明，加權(quán)總體最小二乘法的精度最高，總體最小二乘法和加權(quán)最小二乘法次之，最小二乘法效果最差。加權(quán)總體最小二乘法使兩個區(qū)域的變異函數(shù)模型參數(shù)估計精度分別提高70%和60%左右，均方根預(yù)報誤差分別減少5mm 和98mm。圖3和圖4中，WLS和WTLS對應(yīng)于表中的WLS2和WTLS2，可知加權(quán)總體最小二乘法對應(yīng)殘差分布曲線較其他兩種方法變化更平緩，更接近橫坐標(biāo)軸。兩個區(qū)域的數(shù)據(jù)均表明，使用點對數(shù)法對應(yīng)的加權(quán)最小二乘法所得到的結(jié)果優(yōu)于按熵權(quán)法對應(yīng)的結(jié)果，但是兩種定權(quán)方法對應(yīng)的加權(quán)總體最小二乘法結(jié)果卻十分接近。算例中每組數(shù)據(jù)的兩種WTLS法得到的冪指數(shù)α均大于2，考慮到§1.1提到的冪函數(shù)模型適用條件，為保持模型的變異函數(shù)特性，在進行Kriging插值過程中，重新把α取值為1.999 999，這也與文獻［7］的做法一致。兩種WTLS法得到的冪函數(shù)模型常系數(shù)不同、冪指數(shù)相同，這就使得到的RMS 和STD差別很小。數(shù)據(jù)一結(jié)果的精度維持在mm 級，數(shù)據(jù)二結(jié)果的精度由dm 級提高到了cm 級，加權(quán)總體最小二乘法在這兩種情形下都能有效地提高參數(shù)的估計精度。以上結(jié)果證明，把加權(quán)總體最小二乘法引入到變異函數(shù)領(lǐng)域進行參數(shù)估計是可行和有效的。

3 結(jié) 語

本文在考慮距離值誤差的基礎(chǔ)上，通過協(xié)方差傳播律發(fā)現(xiàn)分組后的變異函數(shù)值和距離值是不等精度的，并給出距離值的定權(quán)方法，再結(jié)合熵權(quán)法和點對數(shù)法兩種變異函數(shù)值的定權(quán)方法，把加權(quán)總體最小二乘回歸法引入到變異函數(shù)模型參數(shù)估計中。模擬數(shù)據(jù)和實測數(shù)據(jù)證明了加權(quán)總體最小二乘回歸法的可行性和有效性，相對于最小二乘法、加權(quán)最小二乘法和總體最小二乘法，加權(quán)總體最小二乘回歸法能得到更高精度的變異函數(shù)模型參數(shù)估值。本文僅對冪函數(shù)模型進行了算例討論，而對于其他變異函數(shù)模型的適用性算例驗證以及距離誤差對函數(shù)模型的影響機制、進一步提高加權(quán)總體最小二乘法的參數(shù)估計精度和解算效率，還有待于研究。

［1］徐建華.現(xiàn)代地理學(xué)中的數(shù)學(xué)方法［M］.北京：高等教育出版社，2002（Xu Jianhua.Mathmatica Methods in Contemporary Geography［M］.Beijing：Higher Education Press，2002）

［2］王仁鐸，胡光道.線性地質(zhì)統(tǒng)計學(xué)［M］.北京：地質(zhì)出版社，1988（Wang Renduo，Hu Guangdao.Linear Geostatistics［M］.Beijing：Geological Press，1988）

［3］矯希國，劉超.變差函數(shù)的參數(shù)模擬［J］.物化探測技術(shù)，1996，18（2）：157－161（Jiao Xiguo，Liu Chao.Estimation of Variation Parameter［J］.Computing Techniques for Geophysical and Geochemical Exploration，1996，18（2）：157－161）

［4］李玲，何濤，張武，等.變異函數(shù)線性化的統(tǒng)一參數(shù)估計方法研究［J］.長江大學(xué)學(xué)報：自然科學(xué)版（理工卷），2010，7（2）：127－129（Li Ling，He Tao，Zhang Wu，et al.Study on the Unity Parameter Estimation Method of Linear Variogram［J］.Journal of Yangtze University：Nat Sci Edit，2010，7（2）：127－129）

［5］李明，高星偉，文漢江，等.Kriging方法在GPS水準(zhǔn)擬合中的應(yīng)用［J］.測繪科學(xué)，2009，34（1）：106－107（Li Ming，Gao Xingwei，Wen Hanjiang，et al.The Application of Kriging Method in GPS Leveling Fitting［J］.Science of Surveying and Mapping，2009，34（1）：106－107）

［6］郭泉河，李秀海.不同變異函數(shù)的泛Kriging法的GPS高程擬合結(jié)果［J］.黑龍江工程學(xué)院學(xué)報：自然科學(xué)版，2011，25（4）：26－28（Guo Quanhe，Li Xiuhai.Application of Universal Kriging Technology with Different Semivariation Function Models to GPS Height Anomaly Fitting［J］.Journal of Heilongjiang Institute of Technology，2011，25（4）：26－28）

［7］潘家寶，戴吾蛟，章浙濤，等.變異函數(shù)模型參數(shù)估計的信息熵加權(quán)回歸法［J］.大地測量與地球動力學(xué)，2014，34（3）：125－128（Pan Jiabao，Dai Wujiao，Zhang Zhetao，et al.Parameter Estimation of Variogram Model by Using Information Entropy Weighted Regression［J］.Journal of Geodesy and Geodynamics，2014，34（3）：125－128）

［8］嚴(yán)華雯，吳健平.加權(quán)最小二乘法改進遺傳克里金插值方法研究［J］.計算機技術(shù)與發(fā)展，2012，22（3）：92－95（Yan Huawen，Wu Jianping.Reasearch on Genetic Algorithm Kriging Optimized by Weight Least Square［J］.Computer Technology and Development，2012，22（3）：92－95）

［9］曾懷恩，黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究［J］.測繪工程，2007，16（5）：5－13（Zeng Huaien，Huang Shengxiang.Research on Spatial Data Interpolation Based on Kriging Interpolation［J］.Engineering of Surveying and Mapping，2007，16（5）：5－13）

［10］Golub G H，Loan C V.An Analysis of the Total Least－Squares Problem［J］.SIAM Journal on Numerical Analysis，1980，17（6）：883－893

［11］王樂洋.總體最小二乘解性質(zhì)研究［J］.大地測量與地球動力學(xué)，2012，32（5）：48－52（Wang Leyang.Research on Properties of Total Least Squares Estimation［J］.Journal of Geodesy and Geodynamics，2012，32（5）：48－52）

［12］王樂洋，許才軍.總體最小二乘研究進展［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2013，38（7）：850－856（Wang Leyang，Xu Caijun.Progress in Total Least Squares［J］.Geomatics and Information Science of Wuhan University，2013，38（7）：850－856）

［13］Felus Y A，Schaffrin B.A Total Least－Squares Approach in Two Stages for Semivariogram Modeling of Aeromagnetic Data［C］.IAMG2005，Toronto，2005

［14］王樂洋，魯鐵定.總體最小二乘平差法的誤差傳播定律［J］.大地測量與地球動力學(xué)，2014，34（2）：55－59（Wang Leyang，Lu Tieding.Propagation Law of Errors in Total Least Squares Adjustment［J］.Journal of Geodesy and Geodynamics，2014，34（2）：55－59）

［15］Fang X.Weighted Total Least Squares Solutions for Applications in Geodesy［D］.Hanover：Leibniz University of Hanover，2011

［16］Jazaeri S，Amiri－Simkooei A R，Sharifi M A.Iterative Algorithm for Weighted Total Least Squares Adjustment［J］.Survey Review，2014，46（334）：19－27

［17］朱衛(wèi)東，李全海.基于標(biāo)準(zhǔn)化動量BP神經(jīng)網(wǎng)絡(luò)的GPS高程轉(zhuǎn)換［J］.大地測量與地球動力學(xué)，2010，30（1）：123－125（Zhu Weidong，Li Quanhai.Conversion of GPS Height Based on Standardization Momentum BP Neural Network［J］.Journal of Geodesy and Geodynamics，2010，30（1）：123－125）

［18］黎劍.區(qū)域GPS高程異常擬合及建模方法研究［D］.昆明：昆明理工大學(xué)，2013（Li Jian.Research on Regional GPS Height Anomaly Fitting and Modeling［D］.Kunming：Kunming University of Science and Technology，2013）