司亞君，何燕蘭

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

1 引 言

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

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

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

2 平差模型

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

(1)

(2)

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

(3)

(4)

(5)

(6)

(7)

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

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

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

(8)

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

vTPv=min

(9)

(10)

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

(11)

(12)

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

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

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

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

(13)

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

(14)

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

=min

(15)

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

(16)

(17)

(18)

(19)

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

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

(20)

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

(21)

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

(22)

(23)

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

(22)

(23)

(24)

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

4 算例與分析

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

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

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

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

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

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

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

混入粗差的觀測(cè)值 表3

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

不同算法求解的精度 表4

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

5 結(jié) 論

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

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

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

[3] 亢瑞紅，郭廣禮，胡洪. Helmert定權(quán)在井下平面控制導(dǎo)線(xiàn)網(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] 張國(guó)良，朱家鈺，顧和和. 礦山測(cè)量學(xué)[M]. 徐州：中國(guó)礦業(yè)大學(xué)出版社，2001.

[7] 崔希璋，於宗儔，陶本藻等. 廣義測(cè)量平差[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ì)[J]. 測(cè)繪學(xué)報(bào)，2016，45(3)：297～301.

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