楊 磊 金諸斌 劉祖斌 盧奐采

(浙江工業(yè)大學特種裝備制造與先進加工技術教育部重點實驗室 杭州310014)

0 引言

隨著人類對海洋探索的日益重視,水下聲學定位系統(tǒng)應用日益廣泛,目前已研制成功的有長基線定位系統(tǒng)、短基線定位系統(tǒng)、超短基線定位系統(tǒng),以及近年來提出的變基線定位系統(tǒng)[1]等。但是,聲速受到溫度、鹽度、深度變化的影響導致它是非均勻分布的,因此聲線在海水中的傳播路徑通常也是彎曲的并且由于海面和海底的反射以及折射使得水下聲線傳播產生多徑效應[2]。如果在水聲定位系統(tǒng)中使用平均聲速,則會導致系統(tǒng)定位精度較低,特別是當首達聲線為非直達聲情形時,會嚴重降低系統(tǒng)的定位精度。針對上述問題,近年來科研人員提出了眾多的聲速修正方法。

基于統(tǒng)計學的聲速修正方法得益于人工智能算法的進步,近年得到了較大的發(fā)展。文獻[3]針對長基線定位系統(tǒng),利用粒子群算法求解最優(yōu)參數計算有效聲速。針對短基線定位系統(tǒng),文獻[4]將區(qū)域劃分自適應粒子群優(yōu)化算法用于尋找有效聲速,結果表明該算法可有效提升系統(tǒng)性能。文獻[5]提出利用期望最大化方法來估計水下兩點間的有效聲速以提高單信標系統(tǒng)定位精度。文獻[6]首先將聲速剖面進行分段線性處理,并默認聲線以直線形式傳播,隨后通過數學解析式來描述聲線的傳播過程,進而實現聲速修正。文獻[7]在文獻[6]研究的基礎上認為可以結合聲速反演的方法來進一步提高聲線修正的精度?；谏渚€模型的聲線修正方法由于更加符合聲線在水下傳播的實際過程,應用更加廣泛?？紤]到水中存在的噪聲影響以及時變的水聲環(huán)境,文獻[8]對基于射線模型下直達聲區(qū)域內水平距離和時延估算的克拉美羅界展開了研究,這對于基于射線模型的聲速修正方法精度對比有重要意義。文獻[9]將水下聲線的傳播形式進行分類后,提出通過射線追蹤的方式,建立一個包含不同聲傳播路徑的時延和水平距離信息表格供水下定位使用。Bellhop 是水聲中重要計算工具,文獻[10]直接通過使用該工具計算首達聲線后進行聲速修正。文獻[11]利用二分迭代法計算直達聲區(qū)域的首達聲線后,再得到對應的有效聲速,并通過遺傳算法對區(qū)域有效聲速表進行稀疏后提高實際使用效率。文獻[12]在非直達聲區(qū)域,利用二分迭代法計算兩點間所有本征聲線后,以最短時延得到首達聲線后再計算出對應的有效聲速。利用迭代法解算聲線初始掠射角是一種有效的方法,但是由于聲線在水下傳播路徑的復雜性,常常導致迭代法失效[13],這也是本文討論的重點。

受聲速剖面結構和聲源與接收點間距離影響,兩點間首達聲線會存在直達、反射和折射3 種形式,而不再拘泥于直達聲,此時需要通過計算兩點間所有本征聲線后以最短時延為標準確定首達聲線。但是,在面對聲線傳播經歷多跨度和包含反轉點時,利用迭代法解算聲線初始掠射角失效的現象在文獻中卻鮮有涉及。為了更加精準解算本征聲線初始掠射角,本文提出依據反轉點和多跨度對初始掠射角區(qū)間進行分類,并在包含反轉點的單個跨度角度區(qū)間內采取先預處理后迭代的方式。最后,利用該算法針對南海某地負梯度聲速剖面展開仿真研究,計算結果表明該算法在非直達聲下依然可以精準計算任意一點首達聲線和有效聲速值。

1 聲線修正方法

1.1 有效聲速法

有效聲速法[14]是基于水下聲學射線模型提出的一類聲速修正方法,它利用等效直線代替聲線傳播的曲線軌跡,不僅降低了聲速修正的計算復雜度且有著較高的精度,因此被廣泛應用于水下定位系統(tǒng)中。

在水下聲學射線模型中,聲線的傳播路徑如圖1所示,其過程滿足斯涅耳定律:

圖1 聲線傳播路徑示意圖

其中,c(z) 是深度為z處的聲速值,θ(z) 是深度為z處的掠射角,zS和zP分別是聲源和接收點處深度,k0是滿足斯涅耳定律的聲線常數。聲線在水下傳播的水平距離h和時延t的一般表達式分別為

在不考慮遠距離條件下,水下聲速變化可以采取豎直分層介質模型,而忽略水平變化的影響[15]。為了精確逼近聲線傳播的軌跡,將聲速剖面沿深度方向進行分層,層內默認為恒定聲速梯度,則每層的聲線軌跡為一段圓弧。對聲源和接收點間聲速剖面分為N層,層內聲速梯度為

其聲線傳播的水平距離和時延公式可寫為

有效聲速可由首達聲線傳播的斜距與時延的比值計算得到:

其中,R表示兩點間的真實斜距,t表示利用射線模型計算的首達聲線傳播時延。為了精準計算首達聲線的傳播時延,必須要通過搜索兩點間所有本征聲線后,以最短時延為標準獲取首達聲線的時延信息[16]。

1.2 多跨度

隨著水平距離的不斷增大,聲線會經歷多次反射或反轉形成多跨度后到達接收點。圖2 是南海某區(qū)域深海聲道,首達聲線經過海底反射后到達接收點,此時跨度數是0。圖3 中聲速剖面是文獻[17]提及的一類雙線性聲速剖面,此時首達聲線到達形式是經過2 次海面和海底反射后到達接收點,經歷的跨度數為2。

圖2 首達聲線是經過海底反射的形式(跨度為0)

圖3 首達聲線是經過海面和海底反射的形式(跨度為2)

圖4 是聲源深度在zs=950 m 和接收點深度在zd=380 m 時,聲線經歷不同跨度時傳播的水平距離和初始掠射角的變化曲線。圖中跨度為0、1、2 的角度區(qū)間分別對應[14.31 °,16.73 °]、[16.74 °,28.56°]、[28.57°,38.71°],且每個區(qū)間內的變化曲線是單調的,但是相鄰2 個變化曲線是間斷的,因此采用二分迭代法計算聲線初始掠射角時,若沒有依據不同的跨度對應的角度區(qū)間選取迭代初值,可能會導致迭代區(qū)間受限而無法得到正確的結果。在圖4 中,當聲源和接收點間的真實水平距離為13 300 m時,跨度為0 對應的區(qū)間內并不存在本征聲線;水平距離為10 000 m 時,3 個跨度都存在對應的本征聲線,但需要分類計算。

圖4 跨度為0、1、2 時聲線初始掠射角和水平距離關系

1.3 反轉點

水下聲線傳播形式中的折射會導致聲線反轉并產生反轉點[18],它的特點是聲線在此處的掠射角為0 °,且反轉點深度隨初始掠射角的改變而變化。由斯涅耳定律可得:

式中z0、θ0為聲源處的深度和掠射角。如圖5 是采用圖3 中的聲速剖面,計算得到的首達聲線是經歷2 次反轉和海面反射后到達接收點;圖6 是采用圖4中的聲速剖面,計算得到的首達聲線是經歷反轉和海底反射后到達接受點。

圖5 首達聲線是經過兩次反轉和海面反射的形式

圖6 首達聲線是經過反轉和海底反射的形式

由斯涅耳定律可知,聲線從聲源處出發(fā)經過不同路徑傳播到達接收點,若存在c(z)＞c(z0),且初始掠射角在合適范圍內就會出現反轉點,此處把包含反轉點的反轉區(qū)間臨界角稱為臨界反轉初始掠射角。圖7 和8 是以南海和北極某地的聲速剖面為例,通過改變聲源深度計算得到的臨界反轉初始掠射角。由圖可知,當聲源位于聲速剖面最大點對應的深度處(即c(z0)= max[c(z)])時,臨界反轉初始掠射角為0°,這就表明此時并不存在反轉區(qū)間,無論初始掠射角多大,都不會出現聲線反轉的情形。因此,臨界反轉初始掠射角與聲源位置處聲速呈負相關變化關系。

圖7 南海某地聲速剖面下臨界反轉初始掠射角隨聲源深度的變化圖

圖8 北極某地聲速剖面下臨界反轉初始掠射角隨聲源深度的變化圖

反轉點對計算本征聲線的影響主要有如下2點。

(1)當聲線中存在上(下)反轉點時,聲線在深度方向所能達到的最高(低)點就是反轉點。當初始掠射角在較小區(qū)間時聲線始終無法到達接收點層,但是隨著角度的增大則存在一個最小的能到達接收點層的初始掠射角,如圖9 中的θL,將其定義為最小到達初始掠射角。因此實際計算本征聲線時,可以忽略小于最小到達初始掠射角的角度區(qū)間。

圖9 最小到達初始掠射角示意圖

(2)對于包含反轉點的聲線傳播形式而言,結合式(5)和圖9 可知,由于聲線的反轉點深度隨初始掠射角的變化而變化,使得聲線傳播水平距離不僅是初始掠射角θ0的函數,還受聲速梯度的影響。圖10 是采用圖4 中的聲速剖面,當聲源zs=900 m、接收點zd= 380 m、跨度為0 時,包含反轉點(角度區(qū)間[6.71 °,14.30 °])和不包含反轉點(角度區(qū)間[14.31 °,16.73 °])的聲線傳播水平距離隨初始掠射角變化曲線。由圖可知,包含反轉點的變化曲線不僅呈現出非單調變化趨勢,同時曲線也是非光滑的(這是由于對聲速剖面分層后默認層內聲速梯度是線性所導致的,文獻[19]也對分層效應帶來的誤差進行了分析)。當兩點間水平距離如圖中所示為10 000 m 時,若直接使用二分迭代法解算本征聲線初始掠射角,會出現算法不收斂的情形。由于在水下聲線傳播過程中反轉點只出現在部分初始掠射角范圍內,為了快速、精準地計算本征聲線,只對包含反轉點的聲線計算過程進行預處理即可。因此,有必要將聲線掠射角范圍按照聲線傳播形式中是否包含反轉點進行分類。

圖10 反轉點對聲線初始掠射角與水平傳播距離變化關系的影響

2 改進的非直達聲下聲速修正方法

2.1 本征聲線時延和水平距離的計算

如圖11 所示,聲線傳播形式可由子跨度進行描述[20],其中S1表示從聲源出發(fā)經過海面反射或者上反轉點折射后,回到聲源層的子跨度;S2表示從聲源層到達接收點層的子跨度;S3表示從接收點層出射后經過海底反射或下反轉點折射后,回到接收點的子跨度。根據具體到達形式,本征聲線可分為如下4 種到達形式:第1 類是直達聲;第2 類是經過海面反射或者反轉點折射后到達接受點;第3 類是經過海底反射或者反轉點折射后到達接收點;第4類是經過海面反射或者反轉點折射后,再經過海底反射或者反轉點折射后到達接收點。

圖11 聲線傳播的子跨度

結合式(5)可計算4 種到達形式的聲線傳播水平距離:

一個完整的跨度聲線所傳播的水平距離為

對于存在反轉點的聲線到達形式,應該首先確定反轉點所在的深度位置。當聲源和接收點位置確定后,利用式(9)可以確定不同初始掠射角的聲線對應的反轉點處的聲速c(z反轉)。然后在聲速剖面中尋找c(z反轉) 所處的層間位置,結合式(9)可得反轉點深度:

綜上,任意類型的本征聲線的水平傳播距離可以表示為

傳播時延的計算可以類比聲線傳播水平距離的計算方法,并結合式(6)得到式(14)。

其中,k表示聲線所經歷的完整跨度數,T表示完整跨度所傳播的時延,ti表示其中一種到達形式下去除完整跨度后部分路徑所傳播的時延。

2.2 非直達聲下有效聲速計算

非直達聲下有效聲速的計算前提是獲取首達聲線,因此整個算法從搜索本征聲線開始。聲線初始掠射角向上和向下的角度區(qū)間分別是[-90 °,0]和[0,90 °],斯涅耳定律中使用的是余弦函數,因此在算法中統(tǒng)一使用[0,90 °]計算即可。在給定聲速剖面和聲源位置后,應首先判斷聲線傳播路徑是否存在反轉點。依據1.3 節(jié)可知,當聲源位于聲速剖面最大值對應深度時(c(z0)= max[c(z)]),臨界反轉初始掠射角為0 °,那么聲線不存在反轉點;當聲源不在聲速剖面最大值對應的深度處時,則聲線可能包含反轉點。此時,應假設以接收點層為反轉點深度,計算得到最小到達初始掠射角θ1:

那么在[0,θ1] 初始掠射角區(qū)間內發(fā)出的聲線都無法到達接收點,計算本征聲線時可以忽略。再以聲速剖面最大值max[c(z)] 層為反轉點所在層,得到臨界反轉初始掠射角θ2:

那么在[θ2,90°] 初始掠射角區(qū)間內必然不會再次出現反轉點。綜上所述,可以將聲線初始掠射角范圍初步分為3 個區(qū)間:[0,θ1]、[θ1,θ2]、[θ2,90°]。[0,θ1] 是不存在本征聲線的初始掠射角區(qū)間,[θ1,θ2] 是聲線存在反轉點的區(qū)間,[θ2,90°] 區(qū)間內的聲線不會存在反轉點。

圖12 總結了該有效聲速算法的邏輯流程。其中有反轉點算法(算法1)和無反轉點算法分別是針對聲線包含反轉點和不包含反轉點情形下提出的解算本征聲線初始掠射角方法。有反轉點算法對初始掠射角完成兩次分類后,還需要按照步進角Δθ遍歷計算單個跨度內每個初始掠射角的聲線水平傳播距離(步驟1～6),再將兩點間真實水平距離減去該結果,尋找到差值的零點位置后,將該零點的角度區(qū)間作為迭代角度區(qū)間解算本征聲線初始掠射角(步驟7～22)。此外,二分迭代法中并沒有確定迭代的方向(算法中默認為遞增),因此當差值實際變化為遞減時,需要調換迭代法中的上下限(步驟21)。無反轉點算法可以直接省去有反轉點算法中的步驟5～22,在步驟4 確定聲線初始掠射角范圍后,直接利用二分迭代法(算法2)計算不同情況即可。算法中步進角Δθ的選擇可以依據聲速剖面分層的情況來確定,通常當聲速剖面的相對梯度較大時,建議選擇較小的步長。

圖12 有效聲速計算流程圖

3 算例驗證

本文算例采用了Argo 系統(tǒng)中南海某處1 月平均聲速剖面,將聲源放置在不同位置來計算區(qū)域有效聲速,并對有效聲速空間分布展開分析,驗證本算法的準確性。算例的水平距離范圍是[0,2000] m、深度距離范圍為[0,1000] m,計算步距為Δr=Δz=10 m。

圖13 至圖16 展示了聲源在zs=1025 m 下的區(qū)域有效聲速分布和區(qū)域某點首達聲線到達形式。結合圖13 和圖14 可知,有效聲速的空間分布隨著水平距離增大可以分為3 個區(qū)域,首先經歷平穩(wěn)區(qū)間,然后經歷增大區(qū)間,最后經歷躍變區(qū)間。需要特別說明的是,隨著水平距離的增大,上述的變化歷程會不斷重復。在平穩(wěn)區(qū)間內,首達聲線均為直達形式,此時采用傳統(tǒng)的迭代法即可計算出首達聲線;在增大區(qū)間內,首達聲線在大聲速區(qū)域內經歷了反轉,此時若沒有考慮反轉點對迭代法解算聲線掠射角的影響,則可能導致無法得到正確結果;在躍變區(qū)間,首達聲線為經過反轉和海底反射后到達,且經歷的跨度數為1,這表明在跨度為0 的角度區(qū)間并不存在首達聲線,此時不僅需要考慮反轉點的影響,還需要將初始掠射角依據跨度進行分類,分別計算后得到首達聲線。

圖13 有效聲速空間分布圖

圖14 不同深度下有效聲速隨水平距離變化示意圖

圖15 接收點在zr=200 m、R=800 m 時,兩點間的本征聲線圖(跨度為0)

圖16 接收點在zr=200 m、R=18 000 m 時,兩點間的本征聲線圖(跨度為1)

圖17 和18 是聲源位于水面(zs=0 m)時計算的區(qū)域有效聲速分布圖。此時,由于聲源位于聲速剖面最大值處,那么根據第2 節(jié)分析可知聲線必然不會出現反轉現象。但是,隨著水平距離的增大,首達聲線從直達聲變化為經海底反射后到達(圖19),再變化為經海底反射和海面反射后到達的形式(圖20),聲線經歷的跨度數從0 增加到1,需要根據跨度數對初始掠射角進行分類后再解算本征聲線的初始掠射角。

圖17 有效聲速空間分布圖

圖18 不同深度下有效聲速隨水平距離變化示意圖

圖19 接收點在zr=900 m、R=10 000 m 時兩點間本征聲線圖(跨度為0)

圖20 接收點在zr=900 m、R=12 000 m 時兩點間本征聲線圖(跨度為1)

上述2 個案例計算中,首達聲線都是從直達聲開始,隨著水平距離的增大,如果滿足存在反轉點條件,則會經過反轉后到達;如果不滿足,則會經歷反射后到達;最后都會經歷多次反射或者反轉后到達的情況。計算結果表明,該算法可以準確地解算聲線存在反轉點和多跨度時的本征聲線初始掠射角,并且在非直達聲情形下得到任意一點的有效聲速值。

4 結論

基于射線模型的有效聲速法是通過迭代法解算兩點間所有本征聲線后,再利用首達聲線信息計算出有效聲速值的方式實現聲速修正的。但是,當聲線存在多跨度和反轉點時,由于聲線水平傳播距離隨初始掠射角變化曲線存在分段函數和非單調的現象,使得直接使用迭代法無法得到需要的結果。為了解決上述問題,本文提出了一種非直達聲下的聲速修正方法。

首先,文中指出了反轉點對解算本征聲線的影響,提出了最小到達掠射角和臨界反轉初始掠射角2 個概念。前者指的是可以到達接收點層的最小初始掠射角,小于該角度的區(qū)間在實際計算時可以直接忽略;后者指的是在小于該角度的區(qū)間內,聲線均包含反轉點。并且依據這2 個角度對初始掠射角區(qū)間([0,90 °])進行了初步的分類。與此同時,考慮到聲線包含反轉點時即使在單個跨度對應的角度區(qū)間內,聲線傳播水平距離隨初始掠射角變化曲線為非單調的,算法首先將初始掠射角區(qū)間按照Δθ步進角計算所有角度對應的水平傳播距離,再與兩點間真實水平距離做差后尋找曲線的零點,并且根據實際做差的符號變化來調整迭代法中的上下限。通過這種先預處理再迭代的方式,有效保證了迭代法在聲線存在反轉點時依然可以準確地計算出首達聲線。此外,依據在不同跨度下聲線傳播水平距離隨初始掠射角變化曲線的特點,為了避免迭代區(qū)間受到限制從而無法獲取準確的結果,有必要將初始掠射角區(qū)間再依據不同的跨度數進行二次分類,隨后在不同跨度對應的角度區(qū)間內分別使用迭代法計算。

最后2 個案例的計算結果表明,本文提出的算法可以在聲線存在反轉點和多跨度條件下準確解算出首達聲線,并得到對應的有效聲速值。該方法具有解算速度快且精度較高的優(yōu)點,對于提高水聲定位系統(tǒng)的技術指標有一定的實際價值,因此后續(xù)會進行真實的水下定位實驗,并對不同海洋水下環(huán)境下的應用進行對比研究。