應用抗差估計進行海底勘探實時聲學定位的方法

方一涵 胡文凱 于素君

(①中國地質大學(武漢)自動化學院，湖北武漢430074；②東方地球物理公司裝備服務處裝備研究中心，河北涿州 072751)

0 引言

目前，海底勘探作業是發掘海洋油氣資源的重要技術手段[1-3]。海底勘探施工過程(涉及海底電纜和海洋節點等)中一般采用長基線聲學定位或初至波定位系統輔助確定鋪放于海底的地震檢波器的空間位置[4-7]。聲學定位系統安裝于專門的定位船上，可對鋪于海底的聲學應答器實時定位、指導釋放海底電纜或海洋節點的放樣導航船修正航線，從而保證鋪放的電纜或節點滿足放樣精度要求而廣泛應用。

近年來，學者對海洋長基線聲學定位算法進行了相關研究，并取得了一定成效。基于距離交會定位基本方程，方守川[5]、趙建虎等[6]、劉慧敏等[7]提出了深度約束的交會定位方法，提高了在野外聲學觀測值組成的幾何圖形條件較差的情況下的定位精度，解決了控制點坐標存在垂直解不穩定問題。但其主要適用于海洋工程領域對海底控制點的定位，并不適用于海底勘探實時聲學定位作業。方守川[5]和趙爽等[8]所述的差分定位方法，解決了作業船航跡對稱及聲場入射角不大的情況下的無需聲速剖面的定位問題，提高了水下聲學定位的精度，但它們均是針對野外數據觀測值具有雙邊對稱幾何圖形的情況。劉慧敏等[9]提出了顧及聲線彎曲的淺海多目標水聲定位算法，在聲速測量不準且大入射角觀測數據占比較大的情況下，通過改善觀測方程數學模型顯著提高定位精度。但其假設條件是相同的入射角具有近似的聲線彎曲誤差，這就需定位船圍繞海底應答器全方位地采集聲學信號或雙邊觀測。

目前，實際生產作業中的實時定位軟件普遍采用傳統的最小二乘法估計的定位方法。而上述這些定位方法所處理的野外聲學觀測數據，均是在基于已完成數據采集的前提下進行的，一般適用于海洋工程領域中水底控制網的定位或數據后處理和質量控制工作，這與海底勘探實時聲學定位的單一方向航行觀測的方式不同。

結合實際生產作業過程中聲學數據通過聲學定位系統周期性的采集，其形成的觀測值是定位船單一航行方向采集，并可能存在野值和粗差等不合格的觀測值，從而影響實時定位精度的情況。基于以上問題，本文提出了一種基于抗差估計的海底勘探實時聲學定位方法，首先給定觀測值的可能范圍并對其進行野值剔除處理，然后基于概率統計假設進行粗差探測，最后構建極值函數確定觀測值權陣，通過迭代計算實現實時聲學定位，從而滿足放樣導航船修正航線的需要和保證檢波器放樣精度要求。

1 傳統實時定位原理

海底勘探定位船一般配備全球導航定位系統(Global Navigation Satellite System，GNSS)、測深儀、聲學定位系統及實時定位軟件。其中，GNSS為船載聲學探頭提供絕對坐標，船載聲學探頭走航過程中與海底的應答器進行聲波實時測距，通過不同時刻、不同船位距離交會方式即可確定海底應答器的位置。在海底地震勘探作業中，考慮到潮流影響，為提高放纜船或節點船放樣導航的精度和作業效率，一般都會采用聲學定位船跟隨其后進行實時聲學定位的作業方法(圖1)。

圖1 實時聲學定位示意圖

在這種工作模式下，定位船的實時定位采集軟件根據設定好的定位采集間隔，對在其有效工作范圍內的海底應答器進行呼叫，應答器返回應答信號。

定位船進行定位作業時，在不同時刻和位置會發射聲波，其可獲得接收聲波的海底應答器與船載聲學探頭之間的距離觀測值。類似于測邊網距離交會方式[10]，可構建如下觀測方程，采用最小二乘法估計實時解算海底應答器的空間位置[5,10]

ρvz,i=f(pz,i，pt,i)+δρi

(1)

δρi=δρt,i+δρv,i+ερ，i

(2)

式中：ρvz,i是聲學觀測值，為在i時刻定位船v聲學探頭到水下應答器z之間的觀測距離；f(pz,i,pt,i)為i時刻聲學探頭t初始位置pt,i到水下應答器位置pz,i的幾何距離函數；δρi為測距誤差；δρt,i為時間延遲產生的系統誤差；δρv,i為不同位置聲速結構變化引起的系統誤差；ερ,i為隨機誤差。

在給定初始值pz,0情況下對式(1)線性化可得

ρvz,i-f(pz,0,pt,i)=avz,idpz,i+δρi+bvz,idpt,i

(3)

式中：avz,i為i時刻f(pz,i,pt,i)對pz,i的一次偏導系數；bvz,i為i時刻f(pz,i,pt,i)對pt,i的一次偏導系數；dpt,i為定位船聲學探頭在i時刻的空間位置誤差。若定位船定位網絡利用船體姿態修正后，對最后應答器定位的影響可忽略，即bvz,idpt,i≈0。式(3)中的各項具體形式分別為

(4)

(5)

f(pz,0,pt,i)=

(6)

(7)

其中pz，0=(xz,0,yz,0,zz,0)為應答器的初始位置，一般可由測線理論點坐標(或放樣點的實際坐標)給定，fi為幾何距離函數。

在最小二乘法估計的計算模型中，若定位船在n個時刻采樣了n個航跡點，到每一個水下應答器就有n個觀測距離，誤差方程為

[ρvz,i-f(pz,0,pt,i)]

(8)

根據式(8)進行迭代計算，直到最終待估計的應答器的位置改正量小于所給定的誤差即可。將式(8)寫成矩陣形式，即

(9)

(10)

(11)

(12)

li=ρvz,i-f(pz,0,pt,i)i=1,2,…,n

(13)

根據最小二乘法估計，可得

(14)

式中：P為對角權陣，在實際中可取為相應觀測距離的倒數；li為i時刻觀測值ρ與幾何距離函數f的差值；n為總的觀察值。

2 基于抗差估計的計算方法

考慮到野外數據中可能存在的野值或粗差，為了得到最優且可靠的海底應答器實時聲學定位的估計解，分三個步驟實現定位船實時聲學定位的抗差估計，即通過野值剔除處理、粗差探測和構建極值函數確定觀測值權陣，通過迭代計算求得最優解的抗差估計。

2.1 野值剔除處理

野外作業中，實際不同類型的海底底質對聲波傳播會產生不可定量描述的影響，使聲學定位過程中產生數據野值。為了獲得可靠的高精度聲學定位結果，在進行位置估計前，有必要根據如圖2所示步驟對可能存在的野值進行剔除處理。

2.2 粗差探測

經過野值剔除處理后，所得到的觀測數據序列中可能還存在粗差。如果不處理這些可能存在的粗差，將導致聲學定位計算結果不可靠。根據式(9)整理可得如下方程

(15)

=(I-AN-1ATP)Δ=RΔ

(16)

式中R=I-AN-1ATP

式(16)兩邊取數學期望可得

E(V)=RE(Δ)

(17)

由式(17)可知，若Δ只包含偶然誤差時，不存在粗差，那么E(Δ)=0，E(V)=0也成立。V是Δ的線性函數，故V和Δ的概率統計都服從正態分布規則。

(18)

(19)

通過以上粗差探測法遍歷觀測值數組中的元素，循環一次可發現可能存在的粗差，剔除之，重新進行平差計算，再次構建統計量進行探測直至不再發現為止。

2.3 抗差估計

定位船上的聲學觀測值是依據聲學應答的時間間隔依次進入系統的。經過野值處理和可能存在的粗差探測處理，獲得了一組比較干凈的觀測數據。為了確保最終的計算結果是最優且可靠的，還需設計一種能夠抗差的最優估計方法。

在野外定位系統獲得足夠多的觀測值下，經過野值剔除和粗差探測的處理后，可以獲得其前一次最小二乘法估計解

(20)

(21)

(22)

(23)

式中c為不等于零的正常數。再重新組成誤差方程組，兩種權函數分別進行迭代計算

(24)

(25)

式中k為迭代次數，w取w1或w2。

由圖3可見，權函數1和權函數2都具有以下特性：①均為|vi|的減函數；②都是有界函數；③均為連續函數(在自變量取值范圍內)；④權的值域為(0,1]。

圖3 權函數圖像

由于權因子是|vi|的非增函數，隨|vi|的增大而減小。進一步說明殘差的絕對值越大時，權因子越小，對估計值的貢獻度越小。這兩種權函數迭代計算方法的最優估值計算過程方便、簡單，構建權因子過程參與的變量較少，適于實時的聲學定位計算。

上述對權函數的迭代計算中，若|w(k+1)-w(k)|

3 方法驗證

為進一步驗證上述方法的正確性，分別從模擬仿真和野外實際應用情況測試驗證算法。

3.1 模擬仿真實例

根據實時計算的特點利用聲學定位系統廠商所提供的模擬仿真軟件，模擬了聲學定位船的GPS系統、測深儀及聲學定位設備在野外實時聲學定位的數據采集過程。模擬仿真軟件模擬運行時，定位船由南向北航行(圖4)。在此過程中啟動定位采集軟件發送目標地址組號的應答器的呼叫命令并采集編號為16-2應答器的聲學應答數據，并保存計算機內存中。

圖4 模擬聲學定位觀測值分布

由此過程采集到定位船船載探頭的空間位置、測深儀的測深數據及聲學定位系統輸出的聲學觀測值(圖4)。共得到一組(共15個)聲學定位數據，聲速為1500m/s。為驗證本文算法的有效性，模擬野外作業環境中可能存在的粗差現象，在第4個觀測值上人為賦予了粗差(圖4)。

通過本文算法的粗差探測過程，發現第4序號數據存在粗差，在進行最終的估計之前，予以剔除。對剩下的觀測值重新組成誤差方程組進行最優估計計算，分析比較了傳統的最小二乘法和本文抗差估計方法的定位解(表1)。

由表1可見，權函數1和權函數2的抗差估計比傳統的最小二乘法估計的中誤差分別提高了0.667m和0.786m，兩種權函數抗差估計的內符合精度(測量值之間的偏差)均較好。

表1 不同估計定位解

3.2 野外實際計算應用

采用野外實際項目的數據，進一步驗證算法的有效性。項目施工海域水深為15～20m，聲速為1562m/s。定位船由西向東共采集了組號為306、地址號為1～9應答器的實時聲學定位數據，其定位觀測值分布、野值剔除及粗差探測結果如圖5所示。

經本算法野值剔除處理驗證發現：地址號為1、3、4、6、8應答器的觀測值不存在野值；2號應答器存在3個野值；5號應答器存在11個野值；7號應答器存在36個野值；9號應答器存在58個野值。由圖5中紅色線段可知，在野外作業環境中大入射角的觀測值是野值的可能性更大。

再進行粗差探測驗證得出，所有的應答器觀測值在野值剔除后仍存在粗差(圖5)。地址號為1的應答器觀測值存在3個粗差；2號應答器存在7個粗差；3號應答器存在3個粗差；4號應答器存在2個粗差；5號應答器存在7個粗差；6號應答器存在6個粗差；7號應答器存在8個粗差；8號應答器存在6個粗差；9號應答器存在3個粗差。

圖5 不同地址號應答器的定位觀測值分布、野值剔除及粗差探測結果

然后用三種算法計算得到最優估計結果(圖6)。由圖6a中可見，三種算法所得海底應答器的位置趨勢基本相同，但兩種權函數計算的位置比最小二乘法更相互趨近；由圖6b可見，用兩種權函數計算的點位中誤差結果優于傳統最小二乘法估計解。

圖6 三種算法應答器的計算結果比較

4 結論

本文提出的基于抗差估計的實時聲學定位方法，經模擬仿真實驗和野外實際數據驗證表明，此方法能夠滿足目前海底勘探過程中實時聲學定位的需求，并可以得出如下結論：

(1)對聲學觀測值中可能存在的野值進行了預處理，主要剔除了那些大入射角的聲學觀測值，從而阻斷了其對最終計算結果的干擾；

(2)基于統計假設的粗差檢驗方法，有效地發現了可能存在的粗差并剔除，能夠保留一組不含粗差的較精確的觀測值；

(3)采用抗差自適應權陣構建的估計方法，選取兩種有界函數作為極值函數，得到的抗差估計解比傳統最小二乘法提高了內符合精度，達到了滿意效果。