方一涵 胡文凱 于素君

(①中國(guó)地質(zhì)大學(xué)(武漢)自動(dòng)化學(xué)院，湖北武漢430074；②東方地球物理公司裝備服務(wù)處裝備研究中心，河北涿州 072751)

0 引言

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

近年來，學(xué)者對(duì)海洋長(zhǎng)基線聲學(xué)定位算法進(jìn)行了相關(guān)研究，并取得了一定成效?；诰嚯x交會(huì)定位基本方程，方守川[5]、趙建虎等[6]、劉慧敏等[7]提出了深度約束的交會(huì)定位方法，提高了在野外聲學(xué)觀測(cè)值組成的幾何圖形條件較差的情況下的定位精度，解決了控制點(diǎn)坐標(biāo)存在垂直解不穩(wěn)定問題。但其主要適用于海洋工程領(lǐng)域?qū)５卓刂泣c(diǎn)的定位，并不適用于海底勘探實(shí)時(shí)聲學(xué)定位作業(yè)。方守川[5]和趙爽等[8]所述的差分定位方法，解決了作業(yè)船航跡對(duì)稱及聲場(chǎng)入射角不大的情況下的無(wú)需聲速剖面的定位問題，提高了水下聲學(xué)定位的精度，但它們均是針對(duì)野外數(shù)據(jù)觀測(cè)值具有雙邊對(duì)稱幾何圖形的情況。劉慧敏等[9]提出了顧及聲線彎曲的淺海多目標(biāo)水聲定位算法，在聲速測(cè)量不準(zhǔn)且大入射角觀測(cè)數(shù)據(jù)占比較大的情況下，通過改善觀測(cè)方程數(shù)學(xué)模型顯著提高定位精度。但其假設(shè)條件是相同的入射角具有近似的聲線彎曲誤差，這就需定位船圍繞海底應(yīng)答器全方位地采集聲學(xué)信號(hào)或雙邊觀測(cè)。

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

結(jié)合實(shí)際生產(chǎn)作業(yè)過程中聲學(xué)數(shù)據(jù)通過聲學(xué)定位系統(tǒng)周期性的采集，其形成的觀測(cè)值是定位船單一航行方向采集，并可能存在野值和粗差等不合格的觀測(cè)值，從而影響實(shí)時(shí)定位精度的情況?；谝陨蠁栴}，本文提出了一種基于抗差估計(jì)的海底勘探實(shí)時(shí)聲學(xué)定位方法，首先給定觀測(cè)值的可能范圍并對(duì)其進(jìn)行野值剔除處理，然后基于概率統(tǒng)計(jì)假設(shè)進(jìn)行粗差探測(cè)，最后構(gòu)建極值函數(shù)確定觀測(cè)值權(quán)陣，通過迭代計(jì)算實(shí)現(xiàn)實(shí)時(shí)聲學(xué)定位，從而滿足放樣導(dǎo)航船修正航線的需要和保證檢波器放樣精度要求。

1 傳統(tǒng)實(shí)時(shí)定位原理

海底勘探定位船一般配備全球?qū)Ш蕉ㄎ幌到y(tǒng)(Global Navigation Satellite System，GNSS)、測(cè)深儀、聲學(xué)定位系統(tǒng)及實(shí)時(shí)定位軟件。其中，GNSS為船載聲學(xué)探頭提供絕對(duì)坐標(biāo)，船載聲學(xué)探頭走航過程中與海底的應(yīng)答器進(jìn)行聲波實(shí)時(shí)測(cè)距，通過不同時(shí)刻、不同船位距離交會(huì)方式即可確定海底應(yīng)答器的位置。在海底地震勘探作業(yè)中，考慮到潮流影響，為提高放纜船或節(jié)點(diǎn)船放樣導(dǎo)航的精度和作業(yè)效率，一般都會(huì)采用聲學(xué)定位船跟隨其后進(jìn)行實(shí)時(shí)聲學(xué)定位的作業(yè)方法(圖1)。

圖1 實(shí)時(shí)聲學(xué)定位示意圖

在這種工作模式下，定位船的實(shí)時(shí)定位采集軟件根據(jù)設(shè)定好的定位采集間隔，對(duì)在其有效工作范圍內(nèi)的海底應(yīng)答器進(jìn)行呼叫，應(yīng)答器返回應(yīng)答信號(hào)。

定位船進(jìn)行定位作業(yè)時(shí)，在不同時(shí)刻和位置會(huì)發(fā)射聲波，其可獲得接收聲波的海底應(yīng)答器與船載聲學(xué)探頭之間的距離觀測(cè)值。類似于測(cè)邊網(wǎng)距離交會(huì)方式[10]，可構(gòu)建如下觀測(cè)方程，采用最小二乘法估計(jì)實(shí)時(shí)解算海底應(yīng)答器的空間位置[5,10]

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

(1)

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

(2)

式中：ρvz,i是聲學(xué)觀測(cè)值，為在i時(shí)刻定位船v聲學(xué)探頭到水下應(yīng)答器z之間的觀測(cè)距離；f(pz,i,pt,i)為i時(shí)刻聲學(xué)探頭t初始位置pt,i到水下應(yīng)答器位置pz,i的幾何距離函數(shù)；δρi為測(cè)距誤差；δρt,i為時(shí)間延遲產(chǎn)生的系統(tǒng)誤差；δρv,i為不同位置聲速結(jié)構(gòu)變化引起的系統(tǒng)誤差；ερ,i為隨機(jī)誤差。

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

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

(3)

式中：avz,i為i時(shí)刻f(pz,i,pt,i)對(duì)pz,i的一次偏導(dǎo)系數(shù)；bvz,i為i時(shí)刻f(pz,i,pt,i)對(duì)pt,i的一次偏導(dǎo)系數(shù)；dpt,i為定位船聲學(xué)探頭在i時(shí)刻的空間位置誤差。若定位船定位網(wǎng)絡(luò)利用船體姿態(tài)修正后，對(duì)最后應(yīng)答器定位的影響可忽略，即bvz,idpt,i≈0。式(3)中的各項(xiàng)具體形式分別為

(4)

(5)

f(pz,0,pt,i)=

(6)

(7)

其中pz，0=(xz,0,yz,0,zz,0)為應(yīng)答器的初始位置，一般可由測(cè)線理論點(diǎn)坐標(biāo)(或放樣點(diǎn)的實(shí)際坐標(biāo))給定，fi為幾何距離函數(shù)。

在最小二乘法估計(jì)的計(jì)算模型中，若定位船在n個(gè)時(shí)刻采樣了n個(gè)航跡點(diǎn)，到每一個(gè)水下應(yīng)答器就有n個(gè)觀測(cè)距離，誤差方程為

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

(8)

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

(9)

(10)

(11)

(12)

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

(13)

根據(jù)最小二乘法估計(jì)，可得

(14)

式中：P為對(duì)角權(quán)陣，在實(shí)際中可取為相應(yīng)觀測(cè)距離的倒數(shù)；li為i時(shí)刻觀測(cè)值ρ與幾何距離函數(shù)f的差值；n為總的觀察值。

2 基于抗差估計(jì)的計(jì)算方法

考慮到野外數(shù)據(jù)中可能存在的野值或粗差，為了得到最優(yōu)且可靠的海底應(yīng)答器實(shí)時(shí)聲學(xué)定位的估計(jì)解，分三個(gè)步驟實(shí)現(xiàn)定位船實(shí)時(shí)聲學(xué)定位的抗差估計(jì)，即通過野值剔除處理、粗差探測(cè)和構(gòu)建極值函數(shù)確定觀測(cè)值權(quán)陣，通過迭代計(jì)算求得最優(yōu)解的抗差估計(jì)。

2.1 野值剔除處理

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

2.2 粗差探測(cè)

經(jīng)過野值剔除處理后，所得到的觀測(cè)數(shù)據(jù)序列中可能還存在粗差。如果不處理這些可能存在的粗差，將導(dǎo)致聲學(xué)定位計(jì)算結(jié)果不可靠。根據(jù)式(9)整理可得如下方程

(15)

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

(16)

式中R=I-AN-1ATP

式(16)兩邊取數(shù)學(xué)期望可得

E(V)=RE(Δ)

(17)

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

(18)

(19)

通過以上粗差探測(cè)法遍歷觀測(cè)值數(shù)組中的元素，循環(huán)一次可發(fā)現(xiàn)可能存在的粗差，剔除之，重新進(jìn)行平差計(jì)算，再次構(gòu)建統(tǒng)計(jì)量進(jìn)行探測(cè)直至不再發(fā)現(xiàn)為止。

2.3 抗差估計(jì)

定位船上的聲學(xué)觀測(cè)值是依據(jù)聲學(xué)應(yīng)答的時(shí)間間隔依次進(jìn)入系統(tǒng)的。經(jīng)過野值處理和可能存在的粗差探測(cè)處理，獲得了一組比較干凈的觀測(cè)數(shù)據(jù)。為了確保最終的計(jì)算結(jié)果是最優(yōu)且可靠的，還需設(shè)計(jì)一種能夠抗差的最優(yōu)估計(jì)方法。

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

(20)

(21)

(22)

(23)

式中c為不等于零的正常數(shù)。再重新組成誤差方程組，兩種權(quán)函數(shù)分別進(jìn)行迭代計(jì)算

(24)

(25)

式中k為迭代次數(shù)，w取w1或w2。

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

圖3 權(quán)函數(shù)圖像

由于權(quán)因子是|vi|的非增函數(shù)，隨|vi|的增大而減小。進(jìn)一步說明殘差的絕對(duì)值越大時(shí)，權(quán)因子越小，對(duì)估計(jì)值的貢獻(xiàn)度越小。這兩種權(quán)函數(shù)迭代計(jì)算方法的最優(yōu)估值計(jì)算過程方便、簡(jiǎn)單，構(gòu)建權(quán)因子過程參與的變量較少，適于實(shí)時(shí)的聲學(xué)定位計(jì)算。

上述對(duì)權(quán)函數(shù)的迭代計(jì)算中，若|w(k+1)-w(k)|

3 方法驗(yàn)證

為進(jìn)一步驗(yàn)證上述方法的正確性，分別從模擬仿真和野外實(shí)際應(yīng)用情況測(cè)試驗(yàn)證算法。

3.1 模擬仿真實(shí)例

根據(jù)實(shí)時(shí)計(jì)算的特點(diǎn)利用聲學(xué)定位系統(tǒng)廠商所提供的模擬仿真軟件，模擬了聲學(xué)定位船的GPS系統(tǒng)、測(cè)深儀及聲學(xué)定位設(shè)備在野外實(shí)時(shí)聲學(xué)定位的數(shù)據(jù)采集過程。模擬仿真軟件模擬運(yùn)行時(shí)，定位船由南向北航行(圖4)。在此過程中啟動(dòng)定位采集軟件發(fā)送目標(biāo)地址組號(hào)的應(yīng)答器的呼叫命令并采集編號(hào)為16-2應(yīng)答器的聲學(xué)應(yīng)答數(shù)據(jù)，并保存計(jì)算機(jī)內(nèi)存中。

圖4 模擬聲學(xué)定位觀測(cè)值分布

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

通過本文算法的粗差探測(cè)過程，發(fā)現(xiàn)第4序號(hào)數(shù)據(jù)存在粗差，在進(jìn)行最終的估計(jì)之前，予以剔除。對(duì)剩下的觀測(cè)值重新組成誤差方程組進(jìn)行最優(yōu)估計(jì)計(jì)算，分析比較了傳統(tǒng)的最小二乘法和本文抗差估計(jì)方法的定位解(表1)。

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

表1 不同估計(jì)定位解

3.2 野外實(shí)際計(jì)算應(yīng)用

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

經(jīng)本算法野值剔除處理驗(yàn)證發(fā)現(xiàn)：地址號(hào)為1、3、4、6、8應(yīng)答器的觀測(cè)值不存在野值；2號(hào)應(yīng)答器存在3個(gè)野值；5號(hào)應(yīng)答器存在11個(gè)野值；7號(hào)應(yīng)答器存在36個(gè)野值；9號(hào)應(yīng)答器存在58個(gè)野值。由圖5中紅色線段可知，在野外作業(yè)環(huán)境中大入射角的觀測(cè)值是野值的可能性更大。

再進(jìn)行粗差探測(cè)驗(yàn)證得出，所有的應(yīng)答器觀測(cè)值在野值剔除后仍存在粗差(圖5)。地址號(hào)為1的應(yīng)答器觀測(cè)值存在3個(gè)粗差；2號(hào)應(yīng)答器存在7個(gè)粗差；3號(hào)應(yīng)答器存在3個(gè)粗差；4號(hào)應(yīng)答器存在2個(gè)粗差；5號(hào)應(yīng)答器存在7個(gè)粗差；6號(hào)應(yīng)答器存在6個(gè)粗差；7號(hào)應(yīng)答器存在8個(gè)粗差；8號(hào)應(yīng)答器存在6個(gè)粗差；9號(hào)應(yīng)答器存在3個(gè)粗差。

圖5 不同地址號(hào)應(yīng)答器的定位觀測(cè)值分布、野值剔除及粗差探測(cè)結(jié)果

然后用三種算法計(jì)算得到最優(yōu)估計(jì)結(jié)果(圖6)。由圖6a中可見，三種算法所得海底應(yīng)答器的位置趨勢(shì)基本相同，但兩種權(quán)函數(shù)計(jì)算的位置比最小二乘法更相互趨近；由圖6b可見，用兩種權(quán)函數(shù)計(jì)算的點(diǎn)位中誤差結(jié)果優(yōu)于傳統(tǒng)最小二乘法估計(jì)解。

圖6 三種算法應(yīng)答器的計(jì)算結(jié)果比較

4 結(jié)論

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

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

(2)基于統(tǒng)計(jì)假設(shè)的粗差檢驗(yàn)方法，有效地發(fā)現(xiàn)了可能存在的粗差并剔除，能夠保留一組不含粗差的較精確的觀測(cè)值；

(3)采用抗差自適應(yīng)權(quán)陣構(gòu)建的估計(jì)方法，選取兩種有界函數(shù)作為極值函數(shù)，得到的抗差估計(jì)解比傳統(tǒng)最小二乘法提高了內(nèi)符合精度，達(dá)到了滿意效果。