孫文舟,殷曉冬,暴景陽,曾安敏

1. 海軍大連艦艇學(xué)院軍事海洋與測(cè)繪系，遼寧 大連 116018； 2. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079； 3. 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054

高精度確定海底控制點(diǎn)的三維坐標(biāo)是建立海底大地控制網(wǎng)的關(guān)鍵環(huán)節(jié)[1]。船載GNSS結(jié)合水下聲學(xué)測(cè)距是實(shí)施水下聲標(biāo)(控制點(diǎn))交會(huì)定位的有效手段[2-10]，影響這種方法定位精度最主要的因素是聲速在時(shí)間和空間上復(fù)雜變化所引起的測(cè)距誤差。采用圓走航的方式可以保證水面觀測(cè)點(diǎn)相對(duì)于海底控制點(diǎn)對(duì)稱分布[11]，對(duì)提高控制點(diǎn)水平坐標(biāo)的精確度有重要幫助，然而垂直方向的解依然較差且不穩(wěn)定。

為了提高海底控制點(diǎn)垂直解的精確度，國內(nèi)外學(xué)者進(jìn)行了不同的嘗試。文獻(xiàn)[12]提出了利用觀測(cè)歷元之間差分的方法，這種方法可以有效地消除聲速誤差長(zhǎng)周期項(xiàng)的影響，但是短周期項(xiàng)依然存在。文獻(xiàn)[13]利用傳播時(shí)間建立了聲速誤差模型，提出了在500 m水深情況下，不依賴聲速剖面的水下靜態(tài)目標(biāo)三維坐標(biāo)快速求解方法。文獻(xiàn)[14]認(rèn)為聲速誤差項(xiàng)為常數(shù)，通過南海的試驗(yàn)驗(yàn)證了該方法相比于傳統(tǒng)方法的有效性。近年來，壓力傳感器的測(cè)深精度不斷提高，又有學(xué)者提出了深度約束下的海底控制點(diǎn)坐標(biāo)確定方法[15-17]，進(jìn)一步提高了三維坐標(biāo)解的穩(wěn)健性和精度。

上述方法都是對(duì)聲速變化引起的測(cè)距誤差中可以參數(shù)化的部分進(jìn)行建模改進(jìn)。根據(jù)文獻(xiàn)[2,18]對(duì)水下聲學(xué)測(cè)距誤差的研究，由內(nèi)波引起的測(cè)距誤差形態(tài)復(fù)雜，難以用少量的參數(shù)表達(dá)，而引入?yún)?shù)過多又往往會(huì)導(dǎo)致法方程的病態(tài)甚至秩虧。因此，本文嘗試?yán)冒雲(yún)?shù)平差模型進(jìn)行控制點(diǎn)坐標(biāo)的解算，以期提高垂直解的精確度。

1 海底控制點(diǎn)坐標(biāo)半?yún)?shù)平差模型解算方法

1.1 聲速變化引起的測(cè)距誤差分析

影響聲速變化的主要因素是太陽輻射和海水運(yùn)動(dòng)[19]。聲速剖面的變化包括日變化、季節(jié)變化以及附加之上的隨機(jī)擾動(dòng)，聲速剖面的季節(jié)變化主要與太陽直射點(diǎn)在南北回歸線之間周期性的運(yùn)動(dòng)有關(guān)，而聲速剖面日變化主要與太陽高度角的日變化有關(guān)，由外力擾動(dòng)引起的海洋內(nèi)波造成水體微團(tuán)垂向運(yùn)動(dòng)，使得上下層海水混合引起溫度的變化從而導(dǎo)致海水聲速變化[20]。文獻(xiàn)[2]對(duì)聲學(xué)測(cè)距誤差的研究表明，由聲速變化引起的測(cè)距系統(tǒng)誤差項(xiàng)包含兩部分，長(zhǎng)周期項(xiàng)與短周期項(xiàng)，其中長(zhǎng)周期誤差項(xiàng)與潮汐的影響有關(guān)，而短周期誤差項(xiàng)則主要與海洋內(nèi)波現(xiàn)象有關(guān)。之后文獻(xiàn)[18]在夏威夷島附近進(jìn)行了相關(guān)試驗(yàn)，再一次驗(yàn)證了這個(gè)結(jié)論，并且認(rèn)為聲速主要受溫度變化影響，鹽度變化對(duì)其影響微乎其微[18]。

為進(jìn)一步分析文獻(xiàn)[2,18]的試驗(yàn)結(jié)論，研究聲速對(duì)測(cè)距誤差的影響規(guī)律，本文進(jìn)行了如下試驗(yàn)，首先利用文獻(xiàn)[21]的標(biāo)準(zhǔn)方程構(gòu)建了背景聲速剖面[21-22]

(1)

式中，C1是聲道軸處的聲速；z1是聲道軸的深度；ε=7.4×10-3是擾動(dòng)系數(shù)；B=1300為聲道厚度尺度；z是水層深度。

利用表層溫度的變化規(guī)律以及熱傳導(dǎo)方程構(gòu)建聲速時(shí)空變化的部分[23]

Δc=Δc0e-βzsinωt-βz-φ

(2)

式(2)可表示聲速剖面在不同時(shí)間不同深度的聲速日變化值，通過式(1)和式(2)可構(gòu)造一個(gè)穩(wěn)定變化的聲速剖面。

在海底控制點(diǎn)坐標(biāo)解算的過程中往往采用某一固定的聲速剖面，用上述方法構(gòu)建的聲速剖面簇的平均值作為固定聲速剖面，計(jì)算在水深為3000 m，聲速采樣間隔為5 m的條件下，聲速變化所引起的測(cè)距誤差。為了不引入入射角偏差而導(dǎo)致的誤差，仿真試驗(yàn)采用垂直波束，設(shè)表層溫度日變化的振幅為0.2℃，每小時(shí)采樣一次，共采樣24次，表層海水聲速變化與測(cè)距誤差的計(jì)算結(jié)果如圖1所示。

圖1 表面聲速變化與測(cè)距誤差曲線Fig.1 Surface sound velocity variation and ranging error curve

溫度的日變化只影響表層海水聲速，圖1(a)是對(duì)表層海水等間隔取樣得到的海水聲速平均值在時(shí)間上的變化，圖1(b)是各采樣時(shí)間所對(duì)應(yīng)的測(cè)距誤差的變化。從圖中可以看出，測(cè)距誤差的變化規(guī)律與聲速的變化規(guī)律基本相似，周期相同但相位存在差異。

文獻(xiàn)[2,18]的試驗(yàn)結(jié)果顯示，測(cè)距誤差長(zhǎng)周期項(xiàng)的變化與潮汐的變化規(guī)律相近。這種現(xiàn)象是由于潮波的振動(dòng)產(chǎn)生了日潮內(nèi)波或半日潮內(nèi)波，在潮內(nèi)波的影響下海水產(chǎn)生垂向運(yùn)動(dòng)，從而導(dǎo)致了上層溫度較高的海水與下層溫度較低的海水進(jìn)行熱量交換，引起了海水溫度的變化，溫度變化是引起海水聲速變化最主要的因素，其周期與潮內(nèi)波的周期相同。測(cè)距誤差短周期項(xiàng)的變化與海洋內(nèi)波有關(guān)，與潮內(nèi)波的作用機(jī)理相同，海洋內(nèi)波造成了海水垂向運(yùn)動(dòng)，引起了海水溫度的變化進(jìn)而導(dǎo)致聲速的變化。相比于聲速長(zhǎng)周期項(xiàng)變化，短周期項(xiàng)變化較為復(fù)雜，因此與其對(duì)應(yīng)的測(cè)距誤差短周期項(xiàng)隨時(shí)間變化同樣具有復(fù)雜的形式，難以進(jìn)行參數(shù)化建模。

1.2 圓走航距離交會(huì)定位模型

水下控制點(diǎn)定位的觀測(cè)方程通常表示為

ρi=fXi,Xo+δρdi+δρvi+εi

(3)

式中，ρi是第i時(shí)刻換能器到海底應(yīng)答器的距離觀測(cè)值；f(Xi,Xo)是兩者之間的直線距離，Xi是第i時(shí)刻換能器的位置，Xo是海底應(yīng)答器的位置；δρdi是由應(yīng)答器電路延遲引起的系統(tǒng)性誤差，是可計(jì)算的改正量；δρvi是聲速變化及聲線彎曲引起的系統(tǒng)性誤差，為主要誤差；εi是第i時(shí)刻測(cè)量的隨機(jī)噪聲。

若δρdi已通過外部設(shè)備進(jìn)行改正，則式(3)線性化展開得[14]

(4)

(5)

若進(jìn)行了n次觀測(cè)，可將式(5)寫成矩陣形式

V=BX+δρv-L

(6)

由上文分析可知，由聲速變化引起的測(cè)距誤差變化復(fù)雜，尤其是測(cè)距誤差的短周期項(xiàng)，難以對(duì)其進(jìn)行參數(shù)化的建模，若將每次觀測(cè)的系統(tǒng)誤差δρv都作為未知參數(shù)，則此時(shí)未知參數(shù)的個(gè)數(shù)為n+3，而觀測(cè)方程的個(gè)數(shù)為n，因此式(6)無法得到未知參數(shù)的唯一解。本文采用半?yún)?shù)模型進(jìn)行控制點(diǎn)三維坐標(biāo)的解算。

1.3 控制點(diǎn)坐標(biāo)解算的半?yún)?shù)模型

根據(jù)式(4)，可寫出如下矩陣形式的觀測(cè)方程

L=Bx+S+Δ

(7)

式中，S=[δpv1δpv2…δpvn]T，是由聲速變化引起的測(cè)距系統(tǒng)誤差的n維未知向量。式(7)即為海底控制點(diǎn)坐標(biāo)解算的半?yún)?shù)模型[24-25]。

式(7)的誤差方程是

(8)

由于未知參數(shù)多于方程的個(gè)數(shù)，因此無法直接利用最小二乘原理進(jìn)行求解，需要選擇新的平差準(zhǔn)則[22]

VTPV+αSTRS=min

(9)

式中，α是光滑因子，對(duì)兩參數(shù)V和S起平滑作用；R是一個(gè)給定的正定矩陣，本文中選定R=GTG，其中

同時(shí)，若令所有歷元聲速變化引起的測(cè)距誤差和為0，則又可設(shè)置一個(gè)限制條件

eTS=0

(10)

利用式(8)—式(10)可構(gòu)造拉格朗日函數(shù)

2γeTS

(11)

根據(jù)所構(gòu)造的拉格朗日函數(shù)即可推導(dǎo)出式(8)—(9)的解，這里直接給出推導(dǎo)結(jié)果

(12)

將計(jì)算得到的坐標(biāo)改正數(shù)加上所設(shè)定坐標(biāo)初始值即可獲得海底控制點(diǎn)的三維絕對(duì)坐標(biāo)。

(13)

式中，δpvl是測(cè)距誤差的長(zhǎng)周期項(xiàng)；δpv0是長(zhǎng)周期項(xiàng)測(cè)距誤差的振幅；T是測(cè)距誤差的長(zhǎng)周期項(xiàng)的周期；φ是初始相位。對(duì)長(zhǎng)周期項(xiàng)參數(shù)化處理后，短周期項(xiàng)仍然認(rèn)為不可參數(shù)化，同樣根據(jù)式(12)計(jì)算坐標(biāo)改正數(shù)，加上坐標(biāo)初始值獲得控制點(diǎn)三維絕對(duì)坐標(biāo)。

2 仿真試驗(yàn)分析

為了驗(yàn)證半?yún)?shù)平差模型解算海底控制點(diǎn)坐標(biāo)的有效性，本文進(jìn)行了如下仿真試驗(yàn)：水深設(shè)置為3000 m，應(yīng)答器布放于海底，測(cè)量船以應(yīng)答器為中心，半徑為3000 m的圓航跡航行，速度為6節(jié)(約3 m/s)，每10 s實(shí)施一次水聲測(cè)距。解算的坐標(biāo)系采用笛卡兒坐標(biāo)系，應(yīng)答器為原點(diǎn)，x軸指向東，y軸指向北，z軸指向天頂。測(cè)距誤差的仿真采用與文獻(xiàn)[12]相同的方法，表示為

(14)

式中，誤差共包含4項(xiàng)，第1項(xiàng)和第2項(xiàng)分別是聲速變化引起的短周期項(xiàng)和長(zhǎng)周期項(xiàng)，第3項(xiàng)是測(cè)區(qū)相關(guān)性誤差，第4項(xiàng)為隨機(jī)性誤差。參照文獻(xiàn)[2]在北太平洋的試驗(yàn)和文獻(xiàn)[18]在夏威夷島附近的試驗(yàn)，測(cè)距誤差中的長(zhǎng)周期項(xiàng)與潮汐變化規(guī)律近似，周期為12 h，振幅約20 cm，由內(nèi)波引起的測(cè)距誤差短周期項(xiàng)的周期從幾十分鐘到幾個(gè)小時(shí)不等，振幅約為12 cm，因此令c1=12 cm、c2=20 cm、c3=2 cm，Tw=20 min，隨機(jī)誤差滿足方差為5 cm的高斯分布，t0是初始時(shí)刻，t是任意時(shí)刻，x是t時(shí)刻測(cè)量船的三維坐標(biāo)，x′是海底應(yīng)答器的三維坐標(biāo)。除此之外，設(shè)置潮汐的周期為12 h，振幅為5 m，海面波浪的周期為20 s，振幅為2 m，測(cè)量船水平方向的定位精度為5 cm，垂直方向定位精度為10 cm。

試驗(yàn)中將本文的方法與其他兩種方法進(jìn)行對(duì)比，方法1是直接對(duì)觀測(cè)值利用最小二乘進(jìn)行解算，方法2是采用文獻(xiàn)[12]所提出的差分算法進(jìn)行解算，本文方法為方法3。

第1次仿真試驗(yàn)的采樣數(shù)為4320，采樣時(shí)長(zhǎng)剛好為測(cè)距誤差長(zhǎng)周期項(xiàng)的一個(gè)周期，試驗(yàn)中分別考慮存在內(nèi)波和不存在內(nèi)波兩種情況，表1給出了不同方法解算所得控制點(diǎn)坐標(biāo)偏差的統(tǒng)計(jì)參數(shù)。

表1 不同方法解算所得控制點(diǎn)坐標(biāo)偏差的統(tǒng)計(jì)參數(shù)

對(duì)比水平方向和垂直方向坐標(biāo)偏差的精確度(MSE)可以發(fā)現(xiàn)，聲速誤差對(duì)控制點(diǎn)垂直解的精確度影響較大，而對(duì)水平解的精確度影響小，這與圓走航測(cè)距定位消除水平方向誤差的設(shè)計(jì)相吻合，也被國內(nèi)外研究所證實(shí)。對(duì)比兩個(gè)試驗(yàn)采用差分算法計(jì)算的垂直坐標(biāo)偏差的精確度，可驗(yàn)證文獻(xiàn)[12]的結(jié)論，即差分可以消除測(cè)距誤差中長(zhǎng)周期項(xiàng)的影響，但是短周期項(xiàng)的影響依舊存在。當(dāng)存在內(nèi)波時(shí)，水平方向依舊可以達(dá)到較高精確度，但是垂直方向解的精確度會(huì)明顯變差。將3種方法進(jìn)行橫向比較，當(dāng)不存在內(nèi)波時(shí)，差分算法是3種方法中垂直解精確度最高的，高于半?yún)?shù)模型和最小二乘法，而在實(shí)際的海洋環(huán)境中，更普遍的情況是存在海洋內(nèi)波，此時(shí)，差分算法垂直解的精確度降低，而最小二乘法與半?yún)?shù)模型解的精確度較高。

第2次試驗(yàn)的采樣數(shù)為720，采樣時(shí)長(zhǎng)為測(cè)距誤差長(zhǎng)周期項(xiàng)周期的1/6，試驗(yàn)中同樣分別考慮存在內(nèi)波與不存在內(nèi)波兩種情況，不同方法解算控制點(diǎn)坐標(biāo)偏差的統(tǒng)計(jì)參數(shù)見表2。

表2 不同方法解算所得控制點(diǎn)坐標(biāo)偏差的統(tǒng)計(jì)參數(shù)

將試驗(yàn)2與試驗(yàn)1垂直坐標(biāo)偏差的精確度進(jìn)行對(duì)比可得，無論是否存在內(nèi)波，當(dāng)觀測(cè)時(shí)長(zhǎng)不是測(cè)距誤差長(zhǎng)周期項(xiàng)的整數(shù)倍時(shí)，最小二乘法垂直解的精確度都會(huì)降低。差分方法的結(jié)論與前一試驗(yàn)一致，整體精確度降低，這主要與觀測(cè)數(shù)據(jù)的數(shù)量減少有關(guān)。半?yún)?shù)模型的垂直解則依舊可以得到較高的精確度。圖2—圖5中測(cè)距誤差是指采用文獻(xiàn)[12]的方法仿真得到的真實(shí)測(cè)距系統(tǒng)誤差，而誤差估值是指半?yún)?shù)模型在不同試驗(yàn)條件下對(duì)每次測(cè)距誤差的估值，可以看出半?yún)?shù)模型可以較準(zhǔn)確地描述測(cè)距誤差的變化規(guī)律，對(duì)比圖2和圖4、圖3和圖5，可以發(fā)現(xiàn)在內(nèi)波存在試驗(yàn)中，半?yún)?shù)模型對(duì)測(cè)距誤差估值的結(jié)果明顯發(fā)散，精度明顯下降，表現(xiàn)為圖3中測(cè)距誤差與誤差估值的中誤差分別為0.17和0.22 cm；圖5中測(cè)距誤差與誤差估值的中誤差分別為0.10和0.18 cm。對(duì)比表2中半?yún)?shù)模型在存在和不存在內(nèi)波兩種條件下垂直坐標(biāo)偏差的精確度，可以得出內(nèi)波存在時(shí)垂直解誤差的精確度由6.43增加到25.53，證明內(nèi)波會(huì)影響半?yún)?shù)模型解算的精確度，但依舊小于最小二乘法的253.37和差分算法的389.20。從表2中還可以看出，兩種條件下，雖然最小二乘法垂直解的精度高于半?yún)?shù)模型，但由于其準(zhǔn)確度較低，所以其精確度遠(yuǎn)小于半?yún)?shù)模型。第2次試驗(yàn)對(duì)測(cè)距誤差長(zhǎng)周期項(xiàng)進(jìn)行了參數(shù)化處理，結(jié)果顯示，無內(nèi)波時(shí)對(duì)振幅的估計(jì)為20.02 cm，而當(dāng)內(nèi)波存在時(shí)，對(duì)振幅的估計(jì)為15.92 cm，相對(duì)于所設(shè)定的值20 cm，估值的準(zhǔn)確性明顯降低，對(duì)應(yīng)垂直方向坐標(biāo)解的誤差也增大，這個(gè)結(jié)果表明，當(dāng)測(cè)距誤差同時(shí)包含參數(shù)項(xiàng)和非參數(shù)項(xiàng)時(shí)，非參數(shù)項(xiàng)誤差會(huì)影響對(duì)參數(shù)項(xiàng)估值的準(zhǔn)確性。

圖2 半?yún)?shù)模型對(duì)測(cè)距誤差的估值(無內(nèi)波)Fig.2 Ranging error estimated by semi-parametric model (no internal waves)

圖3 半?yún)?shù)模型對(duì)測(cè)距誤差的估值(有內(nèi)波) Fig.3 Ranging error estimated by semi-parametric model (with internal waves)

圖4 半?yún)?shù)模型對(duì)測(cè)距誤差的估值(無內(nèi)波)Fig.4 Ranging error estimated by semi-parametric model (no internal waves)

圖5 半?yún)?shù)模型對(duì)測(cè)距誤差的估值(有內(nèi)波) Fig.5 Ranging error estimated by semi-parametric model (with internal waves)

在控制點(diǎn)三維坐標(biāo)解算的過程中，最小二乘法和差分算法要想達(dá)到最優(yōu)解，都需要較為苛刻的條件，最小二乘法要求觀測(cè)的時(shí)長(zhǎng)為測(cè)距誤差長(zhǎng)周期項(xiàng)的整數(shù)倍，而差分算法要求不存在短周期項(xiàng)的測(cè)距誤差，前者顯然實(shí)際工程應(yīng)用價(jià)值較低，因?yàn)闇y(cè)量船連續(xù)圓走航12 h甚至24 h難度大，而后者所需要不存在內(nèi)波現(xiàn)象的海洋環(huán)境條件更加難以滿足。因此半?yún)?shù)模型的優(yōu)勢(shì)則明顯突出出來，它不但可以取得較為精確的坐標(biāo)解算值，而且對(duì)客觀環(huán)境的依賴度相對(duì)較低。

3 結(jié) 論

針對(duì)傳統(tǒng)方法確定海底控制點(diǎn)三維坐標(biāo)垂直解精確度差的問題，本文提出了利用半?yún)?shù)平差模型處理測(cè)距系統(tǒng)誤差的方法，通過分析和試驗(yàn)得出了以下結(jié)論：

(1) 測(cè)距誤差的周期性變化與聲速周期性變化有關(guān)。海水通過渦動(dòng)熱傳導(dǎo)和內(nèi)波所引起對(duì)流熱傳導(dǎo)，將上層海水熱量向下層海水傳遞，引起了溫躍層聲速的變化，從而產(chǎn)生了測(cè)距誤差周期性的變化規(guī)律。

(2) 由海洋內(nèi)波等隨機(jī)海洋現(xiàn)象引起的測(cè)距誤差短周期項(xiàng)(非參數(shù)項(xiàng))不會(huì)影響控制點(diǎn)坐標(biāo)水平解的精確度，但會(huì)影響控制點(diǎn)坐標(biāo)垂直解的精確度。

(3) 利用半?yún)?shù)平差模型可以有效處理由聲速變化引起的測(cè)距系統(tǒng)誤差，從而得出優(yōu)于傳統(tǒng)方法的海底控制點(diǎn)三維坐標(biāo)解算結(jié)果。