袁興明,孫振
(1.山東工業職業學院 建筑與信息工程學院,山東 淄博 256414;2.山東理工大學 測繪工程系,山東 淄博 255049)
大地測量與導航定位中,測距定位方程為超定非線性方程,可采用非線性最小二乘理論進行求解.該方程傳統上依據泰勒級數展開取至一階項,采用線性平差估計求解,這會導致信息量的缺失和模型特征的改變[1-2].當非線性模型強度和殘差較大時,傳統的線性化平差估計并不能取得有效的結果[3-4].若此時伴隨不適定問題,會進一步降低參數估計解的精度和可靠性.固有曲率和參數效應是刻畫非線性模型線性化的數量指標,能夠評估推斷效果的優劣程度、適應條件和容許誤差,但無法評估非線性平差的收斂平穩情況[5-6].研究表明,非線性擾動主要來源于線性近似時系數矩陣的擾動、附加的截斷誤差及正交過程[7-8].
迭代估計是求解非線性無約束最優化方法常用的數值方法,即通過構造一定的迭代序列使其滿足一定的收斂條件.雖然牛頓法能夠處理非線性最小二乘問題,但由于其儲存成本等因素,更傾向于采用高斯-牛頓法[9-10].研究表明,高斯-牛頓法的適應條件為殘差小和非線性強度弱的非線性方程,可是當非線性方程含有不適定問題時,傳統的線性化平差估計和經典數值方法如高斯-牛頓法會由于觀測數據誤差擾動而產生強烈的不穩定特征,導致參數估計解不穩定,甚至求解失效[11-14].由此可知,根據實際問題的線性化程度、秩虧或者病態程度,選擇合適的非線性平差模型就顯得尤為重要.
對于測距定位問題,非線性最小二乘解是觀測向量末端以觀測權為質量質點系的重心.針對測距定位方程最小二乘解性質,來構造非線性數值算法的研究較為困難.在這方面,Xue等[15]作出了一定的突破,他基于非線性最小二乘解性質發展了簡單,且無需矩陣求逆的穩定數值方法,即重心迭代法,然而,該方法由于線性搜索因子1/tr(P)過于保守,收斂速度可能非常緩慢,具有較低的收斂效率,實用價值不高.為增加重心迭代法的實用價值,本文基于殘差最小步長準則提出了一種改進的重心迭代法,即松弛重心迭代法來提高重心迭代法的收斂效率.最后采用全球衛星導航系統(GNSS)數據和水下定位數據,驗證了該方法的主要結果.
設測距定位觀測方程的向量表達式為

式中:L為m×n的觀測向量;d(x)=∥xi?x∥2=為已知坐標至未知坐標的歐式距離;x=[x1,x2,···,xm]∈Rm為未知點構成的向量集;xi=[xi1,xi2,···,xim]∈Rm為已知坐標點構成的向量集;ε為觀測誤差.
非線性最小二乘平差實質是求泛函極小值的最優化問題,常用來處理式(1),則非線性最小二乘參數估計解為

式中,


為殘差向量.運用無約束非線性最小二乘求解式(1),非線性最小二乘目標函數即式(3)滿足如下正交條件方程:

由式(4)和式(5)聯立可得

根據式(6),文獻[15]導出了測距定位方程參數估計的重心迭代法,其迭代公式為
式(8)為重心迭代公式,式(7)實際是一種最速下降法,即選取重心迭代法不依賴于初始值的精度,具有計算簡單,無需矩陣取逆、計算海森矩陣的優點.重心迭代法實質上是一種最速下降法,收斂速度主要依賴于線性搜索因子1/tr(P)的影響.1/tr(P)是重心法求解的保守策略,易受到觀測數據個數n的影響,即隨著觀測數目增多,其線性搜索策略可能越保守,比如當權陣為單位陣,觀測數目n為20時,重心迭代法的步長為0.05,這時重心迭代法具有較低的收斂效率.為增加重心法的實用性,必須提高其收斂效率.為此,本文將對重心迭代法進行如下改進

式(9)為松弛重心迭代公式,其中,ωk為松弛參數,主要依據殘差性質來自適應調整迭代步長,提高數值解的整體效果.
文獻[16]提出了殘差最小準則確定步長的思想,驗證了方法的計算效率.為此,本文采用了殘差最小準則來確定松弛因子.殘差最小準則方法充分采用目標函數的導數信息和函數值信息,相比于傳統線搜索方法如黃金分割法和拋物線法具有一定優勢.下面給出了該方法確定松弛參數的嚴密推導公式.設第k+1次的擬合殘差為松弛因子的函數,即

式中,d(xk+1)為xk+1處的函數計算值,可將其在xk依據泰勒級數展開得

考慮到xk+1=xk+ωkd(xk),聯合式(1)和式(2),可將松弛因子函數化簡為

將式(9)代入式(12)聯合式(4)可得

將式(13)關于ω求其偏導并進行化簡可得

若要k+1次殘差最小,需滿足 ?R(ω)/?ω=0,可得松弛因子的確定公式.考慮到迭代過程,可將步長記為如下形式

為了測試本文提出的新方法在觀測病態和良態的表現效果,分別選取GNSS觀測數據和南海水下定位數據進行測試.本文設置迭代收斂條件為該條件所消耗成本遠高于重心法和松弛重心法的迭代成本,雖然本文用該條件來驗證該方法.
本文采用2017年4月28日采集的GNSS觀測數據(該數據文件包括GPS、BDS及GLONASS的三系統觀測值),采樣間隔為1 s,選取其中1個歷元北斗衛星導航系統(BDS)觀測數據進行偽距單點定位實驗.表1給出了原始偽距觀測值和誤差改正后的偽距觀測值(誤差改正包括:電離層改正、對流層改正和衛星鐘差改正等).迭代初值按照文獻[15]給出的方法進行計算.

表1 單歷元的BDS觀測數據
分別采用高斯-牛頓迭代法、重心迭代法和松弛重心迭代法來處理水下實測數據,解算結果如表2所示.由表可知,高斯-牛頓法、重心法和松弛重心法的數值收斂解相同但耗時最短.該算例表明,高斯-牛頓法相對于重心法和松弛重心法具有更好的局部收斂性質,這是因為BDS衛星距離位置待測點相差3萬多千米,非線性強度弱.松弛重心迭代法依據殘差準則自適應確定松弛因子,提高了重心迭代法的收斂效率,節省了計算成本.三種數值相比較而言,高斯-牛頓法計算成本較低,建議當距離觀測方程良好且距離觀測量非常大時,采用高斯-牛頓法.

表2 不同算法的解算結果
圖1給出了重心法和松弛重心法的點位迭代序列圖.圖中k表示迭代次數.由圖可知,松弛重心迭代法依據殘差準則自適應確定松弛因子,不會對最終數值收斂解和迭代序列穩定性產生影響,而明顯提高了重心迭代法的收斂,節省了運算成本.以上表明,松弛重心迭代法在繼承了重心迭代法迭代序列穩定,無需計算海森矩陣的優點之外,相對于重心迭代法具有更好的局部收斂性質,實用價值更高.


圖1 重心迭代算法和松弛重心迭代算法的點位迭代序列圖
本文采用水下定位實測數據進行驗證,數據主要是測量船圍繞應答器航行獲取.在1圈數據中選取15個相鄰數據進行驗證,數據坐標如圖2所示.假設觀測距離為等精度觀測,要求根據15個已知點坐標和觀測距離求待測點的坐標.
以1圈數據的單點定位結果作為真值.數據采集過程中,測量船在水平面航行,Z方向上具有較強的病態性;此外,測量船航行過程中,采集數據相隔時間較短,X方向和Y方向也會產生微弱的病態性.迭代初值按照文獻[15]給出的方法進行計算.經計算初始設計矩陣條件數為cond(N(x0))=3.97×1013.分別采用高斯-牛頓法、重心迭代法和松弛重心迭代法來處理水下實測數據,解算結果如表3所示.

圖2 水下實測數據平面坐標變化圖

表3 不同算法的解算結果
由表3可知,高斯-牛頓法易受到線性初值和方程組模型態性影響,迭代序列受到觀測數據誤差擾動而無法收斂,建議處理病態問題時一般不采用高斯-牛頓法.重心迭代法和松弛重心迭代法無需矩陣取逆,能夠避免觀測誤差擾動對取逆算子擾動的影響.松弛重心迭代法根據殘差準則來自適應確定松弛因子,主要是用來提高重心迭代法的收斂效率,并不會對數值收斂結果產生影響,其結果與重心迭代法的結果相同.由表進一步能夠看出,重心迭代法具有較低的收斂效率,其消耗時間為132.4420 s,迭代次數為28280次;松弛重心迭代法明顯提高了重心迭代法的收斂效率,時間大約可節省130 s,迭代次數大約是重心迭代法的1/8.
圖3給出了松弛重心迭代法和重心迭代法的點位迭代序列圖.由圖可知,重心迭代法和松弛重心迭代法在X方向和Y方向點位迭代序列變化較為平緩;在Z方向上,重心迭代法和松弛重心迭代法點位迭代序列整體變化較為平緩,可是分別在迭代1425次和3410次都產生了較大的擾動,這是由于短程水下測距定位方程病態性主要集中于Z方向,由于觀測數據誤差而導致迭代序列產生不穩定特征.以上結果表明,松弛重心迭代法在繼承重心迭代法迭代序列穩定、無需矩陣求逆和計算海森矩陣的優點之外,具有更好的局部收斂性質,實用價值更高.

圖3 松弛重心迭代算法和重心迭代算法的點位迭代序列圖
無論處理病態問題和良態問題,重心迭代法由于線性搜索因子 1/tr(P)過于保守,具有較低的收斂效率.本文提出的方法在繼承了重心迭代法不依賴于初始值的精度,具有計算簡單,無需矩陣求逆、計算海森矩陣的優點外,依據殘差最小準則確定松弛參數來自適應更新迭代步長,明顯提高重心迭代法的收斂速度.算例結果表明:松弛重心迭代法無論處理良態問題和不適定問題,相比于重心迭代法具有更好的局部收斂性質,應用價值更高.
高斯-牛頓法由于忽視了距離觀測方程的二階項信息,適合處理殘差小且觀測態性良態測距方程.但是當處理不適定測距方程時,高斯-牛頓法由于設計矩陣秩虧或者病態而解算失敗.本文提出的新方法在處理不適定問題時,相比于高斯-牛頓法有更好的收斂效果.為此,我們應該充分考慮松弛重心迭代法和高斯-牛頓法的優點來處理盡可能多的問題.