趙英文 王樂(lè)洋，2，3

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

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

3 江西省數(shù)字國(guó)土重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號(hào)，330013

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

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

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

1.1 變異函數(shù)模型

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

式中，h為兩點(diǎn)間的分隔距離，N（h）為距離為h的點(diǎn)對(duì)數(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ù)。由文獻(xiàn)［1］知，α必須小于2；若α≥2，則冪函數(shù)不再是變異函數(shù)。

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

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

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

令

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

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

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

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

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

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

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

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

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

該算法的迭代過(guò)程為［16］：

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

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

1.3 權(quán)的確定

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

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

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

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

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

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

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

在計(jì) 算 中，取Py＝diag（Wi），PA＝diag

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

2 算例分析

2.1 模擬算例

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

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

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

2.2 高程異常插值算例

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

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

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

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

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

式中，Zi、Zi′分別為k個(gè)檢核點(diǎn)的觀測(cè)值和通過(guò)Kriging插值計(jì)算得到的檢核點(diǎn)估計(jì)值，ri是檢核點(diǎn)預(yù)報(bào)殘差，是檢核點(diǎn)預(yù)報(bào)殘差的均值。結(jié)果如表2、表3所示。

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

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

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

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

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

3 結(jié) 語(yǔ)

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

［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)計(jì)學(xué)［M］.北京：地質(zhì)出版社，1988（Wang Renduo，Hu Guangdao.Linear Geostatistics［M］.Beijing：Geological Press，1988）

［3］矯希國(guó)，劉超.變差函數(shù)的參數(shù)模擬［J］.物化探測(cè)技術(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ì)方法研究［J］.長(zhǎng)江大學(xué)學(xué)報(bào)：自然科學(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］.測(cè)繪科學(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é)報(bào)：自然科學(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ù)估計(jì)的信息熵加權(quán)回歸法［J］.大地測(cè)量與地球動(dòng)力學(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ìn)遺傳克里金插值方法研究［J］.計(jì)算機(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］.測(cè)繪工程，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］王樂(lè)洋.總體最小二乘解性質(zhì)研究［J］.大地測(cè)量與地球動(dòng)力學(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］王樂(lè)洋，許才軍.總體最小二乘研究進(jìn)展［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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］王樂(lè)洋，魯鐵定.總體最小二乘平差法的誤差傳播定律［J］.大地測(cè)量與地球動(dòng)力學(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)化動(dòng)量BP神經(jīng)網(wǎng)絡(luò)的GPS高程轉(zhuǎn)換［J］.大地測(cè)量與地球動(dòng)力學(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）