司亞君，何燕蘭

(江蘇省地質(zhì)測繪院，江蘇 南京 211102)

1 引 言

由于受作業(yè)環(huán)境的限制，導(dǎo)線測量仍然是井下平面控制測量的主要作業(yè)方法，其成果精度的高低是影響井上、井下測量基準(zhǔn)一致性的重要因素。突發(fā)災(zāi)害發(fā)生時，井上、下測量基準(zhǔn)的統(tǒng)一對確定事故發(fā)生地點、救生通道的順利打通有重要意義。同時，導(dǎo)線測量成果也是礦區(qū)“一張圖”建設(shè)、數(shù)字礦山建設(shè)的基礎(chǔ)數(shù)據(jù)[1，2]。為獲得較高精度的礦區(qū)測量數(shù)據(jù)，在應(yīng)用現(xiàn)代測量手段、提高測量儀器的硬件性能，以獲取更高精度的礦山測量外業(yè)數(shù)據(jù)的同時，研究人員一直在探討與改進礦山測量數(shù)據(jù)處理的數(shù)學(xué)模型與方法[3～5]。數(shù)學(xué)模型的有效性與適用性是測量數(shù)據(jù)處理精度的關(guān)鍵，傳統(tǒng)加測陀螺邊的井下導(dǎo)線控制網(wǎng)的數(shù)據(jù)處理方法是在最小二乘準(zhǔn)則(least squares，LS)的基礎(chǔ)上，建立高斯-馬爾可夫模型(Gauss-Markov model，G-M model)求解參數(shù)[6]。當(dāng)觀測向量含有的誤差滿足正態(tài)隨機分布時，最小二乘估計為最優(yōu)、無偏估計。受儀器觀測精度、礦區(qū)作業(yè)環(huán)境的限制以及數(shù)據(jù)處理模型線性化等因素的影響，使得觀測向量中含有誤差的同時，模型的系數(shù)矩陣中也含有誤差，并且可能含有粗差；此時，應(yīng)用最小二乘估計求解參數(shù)不再具有最優(yōu)、無偏的特征[7]。因此，在最小二乘準(zhǔn)則基礎(chǔ)上進行井下導(dǎo)線網(wǎng)平差計算的傳統(tǒng)方法具有一定的缺陷。

為解決導(dǎo)線網(wǎng)數(shù)據(jù)處理模型的觀測向量與系數(shù)矩陣中同時含有誤差的問題，將總體最小二乘準(zhǔn)則(total least squares，TLS)引入導(dǎo)線網(wǎng)測量數(shù)據(jù)處理?？傮w最小二乘是由Adcock為解決變量中含有誤差(error-in-variables，EIV)的問題于1877年提出[8]；1980年，Golub和Van Loan提出可行性較強的基于奇異值分解算法[9](singular value decomposition，SVD)，總體最小二乘準(zhǔn)則在圖像處理、信號處理、化學(xué)工程等專業(yè)開始有廣泛應(yīng)用[10～12]；2008年，Burkhard Schaffrin和Andreas Wieser提出適合處理加權(quán)條件下的基于拉格郎日函數(shù)的總體最小二乘算法[13]，總體最小二乘受到測繪工作者更多的關(guān)注與研究[14～16]。

在建立加測陀螺邊井下導(dǎo)線控制網(wǎng)間接平差模型的基礎(chǔ)上，為實現(xiàn)對觀測向量與系數(shù)矩陣中含有的誤差同時進行最小化的約束，基于拉格郎日函數(shù)，給出模型參數(shù)的總體最小二乘迭代算法。為克服觀測向量與系數(shù)矩陣中含有粗差對模型參數(shù)求解的影響，應(yīng)用Huber權(quán)函數(shù)[17]，對觀測向量與系數(shù)矩陣的權(quán)矩陣進行迭代定權(quán)，降低含有粗差的觀測值對參數(shù)求解的影響。應(yīng)用某礦區(qū)兩井定向的導(dǎo)線控制網(wǎng)數(shù)據(jù)對提出的算法進行驗證，并與傳統(tǒng)算法在精度上進行比較。

2 平差模型

設(shè)兩相鄰控制點的水平距離與方位角分別為Si、ai，兩相鄰控制點的平面坐標(biāo)為(xiyi)、(xi+1yi+1)，則：

(1)

(2)

為應(yīng)用線性估計模型求解參數(shù)，須將式(1)，式(2)分別線性化：

(3)

(4)

(5)

(6)

(7)

上諸式中，lsi、lβi、lai分別為邊長觀測值、角度觀測值、陀螺邊的方位角初值減去其觀測值。

圖1井下導(dǎo)線控制網(wǎng)

現(xiàn)有一井下導(dǎo)線控制網(wǎng)(如圖1所示)，有n個待求控制點，共觀測(n+1)條邊的距離、(n+2)個水平角、加測t個陀螺邊。以待求控制點平面坐標(biāo)改正數(shù)為變量，建立間接平差模型：

(8)

應(yīng)用最小二乘準(zhǔn)則：

vTPv=min

(9)

(10)

式中P為觀測向量的權(quán)矩陣，單位權(quán)方差與參數(shù)的協(xié)方差矩陣：

(11)

(12)

式(11)中，r為多余觀測量。

應(yīng)用最小二乘準(zhǔn)則求解參數(shù)的前提是偶然誤差僅存在于觀測向量中。根據(jù)加測陀螺邊的導(dǎo)線網(wǎng)平差模型(式(8))，模型的系數(shù)矩陣由觀測元素的函數(shù)組成，同時，在對式(1)、式(2)進行線性化的過程中，也會引入舍入誤差。因此，導(dǎo)線網(wǎng)模型的系數(shù)矩陣中也含有誤差；此外，受井下觀測環(huán)境與硬件性能等因素的限制，觀測量中會混入粗差。此時，應(yīng)用最小二乘準(zhǔn)則進行井下導(dǎo)線網(wǎng)平差不再具有無偏、最優(yōu)的性質(zhì)。

3 模型參數(shù)的抗差總體最小二乘算法

為對觀測向量與系數(shù)矩陣中含有的隨機性誤差同時實現(xiàn)最小化約束，克服系數(shù)矩陣與觀測向量中含有粗差對參數(shù)求解的影響，基于迭代選權(quán)的思想，應(yīng)用總體最小二乘準(zhǔn)則求解變量中含有誤差的井下導(dǎo)線網(wǎng)模型待估參數(shù)：

(13)

式中，P1為觀測向量的權(quán)矩陣；b為系數(shù)矩陣B中含有的誤差矩陣，P2為其權(quán)矩陣。并且：

(14)

式中，vec(b)表示對矩陣b進行列向量化；Q1、Q2分別為觀測向量與系數(shù)矩陣的協(xié)因數(shù)陣。根據(jù)拉格郎日函數(shù)，建立如下模型：

=min

(15)

式中，λ為拉格郎日乘數(shù)。要使式(15)最小，則函數(shù)對各變量的偏導(dǎo)數(shù)等于零，即：

(16)

(17)

(18)

(19)

根據(jù)式(16)～式(19)，建立模型參數(shù)的抗差迭代算法[16]。

(1)計算模型參數(shù)的初值

(20)

(2)根據(jù)參數(shù)的初值(式(20))，計算參數(shù)的總體最小二乘解

(21)

式中，Q0為m×m階矩陣，Qb為n×n階矩陣，同時Q0?Qb=Q2(?為矩陣的克羅內(nèi)克積)。協(xié)因數(shù)陣Q1、Q2根據(jù)其權(quán)矩陣計算。應(yīng)用式(16)、式(17)、式(22)計算v、b的初值，根據(jù)給定的限值判斷觀測向量與系數(shù)矩陣中是否含有粗差，并根據(jù)Huber權(quán)函數(shù)迭代定權(quán)[17]

(22)

(23)

令Pi=kPi-1，則觀測向量與系數(shù)矩陣的協(xié)因數(shù)陣可以重新確定，并參與以下的迭代計算。需要指出的是，自適應(yīng)調(diào)節(jié)因子矩陣k(式(22))也須參與以下的迭代計算。

(22)

(23)

(24)

算法的基本步驟為：①在最小二乘準(zhǔn)則下求解參數(shù)的最小二乘解，并據(jù)此計算參數(shù)的總體最小二乘解的初值；②根據(jù)參數(shù)總體最小二乘解的初值以及由拉格郎日函數(shù)導(dǎo)出的相應(yīng)條件，在Huber權(quán)函數(shù)的基礎(chǔ)上，判斷觀測值中是否含有粗差，并重新計算模型的權(quán)陣及其協(xié)因數(shù)陣；③根據(jù)重新確定的協(xié)因數(shù)陣，迭代求解參數(shù)的估計直到收斂于給定的域值；需要注意的是，在迭代過程中，協(xié)因數(shù)陣應(yīng)根據(jù)Huber權(quán)函數(shù)參與迭代計算，以不斷降低含有粗差的觀測值對參數(shù)求解的影響。

4 算例與分析

應(yīng)用某礦的井下導(dǎo)線控制測量數(shù)據(jù)，對論文提出的算法進行驗證。儀器的水平角觀測精度為4″，陀螺經(jīng)緯儀的定向精度為10″，激光測距離儀的測距精度為 2 mm+2 ppm×D。導(dǎo)線網(wǎng)如圖2所示。

圖2 某礦井導(dǎo)線控制網(wǎng)

已知控制點A的平面坐標(biāo)為(xA=50.789 m，yA=8 057.566 m)，控制點B的平面坐標(biāo)為(xB=7.81 m，yB=8 108.035 m)，角度觀測值、已歸算為平面距離的距離觀測值、加測的陀螺邊方位角觀測值如表1所示。

井下導(dǎo)線控制網(wǎng)觀測數(shù)據(jù) 表1

分別應(yīng)用最小二乘準(zhǔn)則與總體最小二乘準(zhǔn)則求解待估控制點的坐標(biāo)，控制點的坐標(biāo)與求解精度(點位中誤差σD)如表2所示。

不同準(zhǔn)則求解的控制點坐標(biāo)及其精度 表2

為應(yīng)用抗差總體最小二乘算法對加測陀螺邊的井下導(dǎo)線網(wǎng)進行平差，將部分觀測數(shù)據(jù)中加入粗差，分析含有粗差的觀測值對不同算法的擾動性影響?；烊氪植畹挠^測值如表3所示。

混入粗差的觀測值 表3

分別應(yīng)用最小二乘估計(LS)、穩(wěn)健估計[16](rubust estimation，RLS)、總體最小二乘估計[13](TLS)、以及本文的抗差總體最小二乘估計(RTLS)對混入粗差的觀測值與其他觀測值組成的導(dǎo)線網(wǎng)進行平差，點位中誤差與單位權(quán)中誤差如表4所示。

不同算法求解的精度 表4

表4的計算結(jié)果表明，含有粗差的觀測值對導(dǎo)線網(wǎng)平差有顯著影響。由于最小二乘算法與總體最小二乘算法不具有抗差性，精度呈現(xiàn)一定的衰減。相比較而言，基于穩(wěn)健估計的最小二乘算法與總體最小二乘算法統(tǒng)計精度較高；由于能夠同時對觀測向量與系數(shù)矩陣中含有的誤差實現(xiàn)最小化的約束，本文的抗差總體最小二乘算法求解精度高于穩(wěn)健最小二乘算法，適用于井下加測陀螺邊的導(dǎo)線網(wǎng)平差。

5 結(jié) 論

在對加測陀螺邊的井下導(dǎo)線控制網(wǎng)建立間接平差模型的基礎(chǔ)上，對由于受觀測環(huán)境與儀器觀測精度的限制，導(dǎo)致數(shù)據(jù)處理模型的觀測向量與系數(shù)矩陣中都會含有誤差，指出應(yīng)用最小二乘準(zhǔn)則處理井下導(dǎo)線控制網(wǎng)數(shù)據(jù)只是對觀測向量中存在的誤差進行最小化約束。為克服傳統(tǒng)數(shù)據(jù)處理方法存在的不足，實現(xiàn)對模型的觀測向量與系數(shù)矩陣中含有的誤差進行最小化的約束，同時顧及粗差對模型參數(shù)的影響，在應(yīng)用拉格郎日函數(shù)的基礎(chǔ)上，建立了基于Huber權(quán)函數(shù)的總體最小二乘抗差迭代算法。結(jié)合某礦的兩井定向數(shù)據(jù)，驗證算法在進行加測陀螺邊的井下導(dǎo)線網(wǎng)平差的有效性，并從點位中誤差與單位權(quán)中誤差兩個角度與其他算法進行了比較。結(jié)果表明，本文算法優(yōu)于其他算法，更適用于井下導(dǎo)線控制網(wǎng)的測量數(shù)據(jù)處理。

[1] 吳立新，殷作如，鄧智毅等. 論21世紀(jì)的礦山——數(shù)字礦山[J]. 煤炭學(xué)報，2000，25(4)：337～342.

[2] 吳立新，殷作如，鐘亞平. 再論數(shù)字礦山：特征、框架與關(guān)鍵技術(shù)[J]. 煤炭學(xué)報，2003，28(1)：1～6.

[3] 亢瑞紅，郭廣禮，胡洪. Helmert定權(quán)在井下平面控制導(dǎo)線網(wǎng)平差中的應(yīng)用[J]. 金屬礦山，2010(11)：104～107.

[4] Guo Jinyun，Guo Shuyang，Liu Xin. On Zero Position Correction of Hanging Tape of Gyro-theodelite[J]. Survey Review，2008，40：129～134.

[5] Wetherelt A，Hunt P. Azimuth Determinations using an Adapted wild GAKI(an Alternative Approach)[J]. Survey Review，2004，37：592～603.

[6] 張國良，朱家鈺，顧和和. 礦山測量學(xué)[M]. 徐州：中國礦業(yè)大學(xué)出版社，2001.

[7] 崔希璋，於宗儔，陶本藻等. 廣義測量平差[M]. 武漢：武漢大學(xué)出版社，2001.

[8] Adcock R J. Note on the method of least squares [J]. Analyst，1877，4：183～184.

[9] Golub G H，Van Loan C F. An Analysis of the Total Least Squares Problem[J]. SIAM J Number，Anal，1980，17：883～893.

[10] 楊鴻森．基于總體最小二乘的紅外圖像去噪[J]. 激光與紅外，2008，38(9)：961～964.

[11] Van Huffel S，Vandewalle J. Analysis and Solution of the Nongeneric Total Least Squares Problem[J]. Siam Jourmal on Matrix Analysis & Applictions，1988，9(3)：360～372.

[12] Lemmerling P，De Moor B，van Huffel S. On the Equivalence of Constrained Total Least Squares and Structured Total Least Squares[J]. IEEE Transactions on Signal Processing，1996，44(11)：2908～2911.

[13] Burkhard Schaffrin，Andreas Wieser. On Weighted total Least-squares Adjustment for Linear Regression[J]. Journal of Geodesy，2008，82：415～421.

[14] Zhang S L，Tong X H，Zhang K L，et al. A Solution to EIV Model with Inequality Constraints and its Geodetic Applications[J]. Journal of Geodesy，2013，87：89～99.

[15] Mahboub V. On Weighted Total Least-squares for Geodetic Transformations[J]. Journal of Geodesy，2012，86：359～367.

[16] 陶葉青，高井祥，姚一飛. 基于中位數(shù)法的抗差總體最小二乘估計[J]. 測繪學(xué)報，2016，45(3)：297～301.

[17] 周江文，黃幼才，楊元喜等. 抗差最小二乘法[M]. 武漢：華中理工大學(xué)出版社，1997.