加測陀螺邊的井下導線控制網抗差總體最小二乘算法

司亞君，何燕蘭

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

1 引 言

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

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

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

2 平差模型

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

(1)

(2)

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

(3)

(4)

(5)

(6)

(7)

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

圖1井下導線控制網

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

(8)

應用最小二乘準則：

vTPv=min

(9)

(10)

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

(11)

(12)

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

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

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

為對觀測向量與系數矩陣中含有的隨機性誤差同時實現最小化約束，克服系數矩陣與觀測向量中含有粗差對參數求解的影響，基于迭代選權的思想，應用總體最小二乘準則求解變量中含有誤差的井下導線網模型待估參數：

(13)

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

(14)

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

=min

(15)

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

(16)

(17)

(18)

(19)

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

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

(20)

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

(21)

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

(22)

(23)

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

(22)

(23)

(24)

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

4 算例與分析

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

圖2 某礦井導線控制網

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

井下導線控制網觀測數據 表1

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

不同準則求解的控制點坐標及其精度 表2

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

混入粗差的觀測值 表3

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

不同算法求解的精度 表4

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

5 結 論

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

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

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

[3] 亢瑞紅，郭廣禮，胡洪. Helmert定權在井下平面控制導線網平差中的應用[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] 張國良，朱家鈺，顧和和. 礦山測量學[M]. 徐州：中國礦業(yè)大學出版社，2001.

[7] 崔希璋，於宗儔，陶本藻等. 廣義測量平差[M]. 武漢：武漢大學出版社，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] 陶葉青，高井祥，姚一飛. 基于中位數法的抗差總體最小二乘估計[J]. 測繪學報，2016，45(3)：297～301.

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