畢思昭 彭朝暉
1) (中國科學院聲學研究所, 聲場聲信息國家重點實驗室, 北京 100190)
2) (中國科學院大學物理學院, 北京 100049)
對于遠程聲傳播問題, 地球曲率的影響不可忽略.為分析地球曲率對遠距離聲傳播的影響, 本文提出了一種地球曲率影響下環境參數的修正方法, 該修正方法無需改動現有聲場計算模型, 具有可移植性好、計算簡便的特點.典型環境下的仿真結果表明, 由于地球曲率的影響, 在會聚區傳播中, 會聚區的位置向聲源方向偏移, 會聚區移動隨著距離變換近似為線性關系, 1000 km距離上偏移可達10 km; 而深海聲道傳播中, 到達結構在時間軸上整體前移, 并伴隨著在深度和時間軸上的擴展, 影響均隨著傳播距離的增加而逐步增大,1000 km距離上整體前移的幅度可達136 ms.
當前我國水聲學研究在深遠海逐漸開展, 認清深海大洋環境下的遠程聲傳播規律對遠程水聲通信、遠程水聲導航等具有重要的意義.
在深遠海環境下, 會聚區傳播和深海聲道傳播具有相對穩定的結構, 能夠高強度、遠距離地傳播聲信號.會聚區位置預報和深海聲道傳播到達結構預報是研究中重點關注的問題.Vadov[1]發現第一會聚區的計算距離要比實驗觀測距離通常多出1—2 km, 他將造成這種差別的原因部分歸結為傳輸條件的日變化.張鵬等[2]利用在中國南海海域收集到的一次深海聲傳播實驗數據, 研究了深海不完全聲道環境下的海底地形對海底反射會聚區位置的影響.吳麗麗[3]利用西太平洋海域長達1029 km的遠程聲傳播數據, 對深海完全聲道下的會聚區傳播損失進行了分析, 并分析了海底深層聲學結構對影區到達結構的影響.Van等[4,5]發現脈沖聲到達結構在深度上會發生擴展, 內波引起的散射很好地解釋了這一現象.張燕等[6]利用一次南海深海聲傳播實驗數據, 分析了深海聲道軸脈沖聲傳播特性,并利用簡正波理論解釋了不同距離和不同深度上脈沖到達結構的產生原因.可以看到, 地形起伏、海底地質類型、傳播路徑上的水文變化會對會聚區傳播、深海聲道傳播產生顯著影響.
當聲傳播距離達到數百甚至上千公里的時候,除了上述因素外, 地球曲率的影響同樣不可忽略.在電磁學和地震學領域, 針對地球曲率對波傳播的影響問題做了充分研究.Richter[7]將共性映射方法運用到電磁波遠距離傳播的地球曲率修正上.Pekeris[8]對電磁學中地球曲率修正方法的適用性進行了分析.Muller[9]給出了地球任意深度點源激發的體波的地球曲率近似方法, 并做出了改進[10].Biswas[11,12]分析了地球曲率對瑞利波、樂甫波的影響.在聲學領域, 徐傳秀等[13]通過坐標變化, 建立了一種拋物方程模型, 并簡要分析了地球曲率對聲傳播的影響, 結果表明地球曲率造成會聚區位置向遠距離移動.Yan[14,15]建立了考慮地球曲率影響的二維和三維的射線聲學模型, 比較了地球曲率修正前后, 本征聲線的到達時間和到達角度, 并發現,地球曲率修正后, 射線路徑向近距離移動.Etter[16]給出了地球曲率的一階修正公式, 但是并沒有給出這一公式的推導過程和基于此修正方法的地球曲率對聲傳播的影響分析.可以看到, 國內外在地球曲率對遠程聲傳播的影響認識上尚有矛盾的地方,還缺乏一種簡便易用的地球曲率修正方法, 對于地球曲率對聲傳播的傳播損失以及對到達結構的影響, 分析得也還不夠到位.
本文借鑒Richer[7]采用的共形映射方法, 推導出一種計算簡便的地球曲率影響下環境參數的修正方法, 通過仿真研究了地球曲率對會聚區傳播的傳播損失以及深海聲道傳播的到達結構的影響, 并對影響規律進行了總結.另外比較了不同地球近似模型下的修正效果.此項研究作為遠距離聲傳播的基礎理論工作, 對認清深海大洋環境下的遠程聲傳播規律, 開展各類應用工作都有重要的意義.
由于地球近似為一個球體, 對于遠程的聲傳播問題, 在距離-深度平面上的聲速等值面不是平行平面, 而是同心球面.當前的聲場計算都是在笛卡爾坐標系之下進行的, 在傳播距離較近的情況下,可以認為地球曲率的影響不存在.而當距離變遠的時候, 由于地球曲率造成的海洋波導的彎曲對聲傳播的影響將不可避免, 此時再用笛卡爾坐標系下建立的各類聲場計算模型進行聲場計算勢必會造成輸入參數與計算模型之間的失配問題.
下面對參數與模型的失配問題進行舉例說明.如圖1所示, 已知聲源的位置和聲源的深度z, 計算與聲源之間距離為r, 深度為z的接收點A處的聲壓.為了更好地說明問題, 在圖1中對海洋波導的厚度進行了放大處理.現有的方法是將參數代入到圖2所示的聲場計算模型進行計算.而實際的聲傳播的情況如圖1所示, 與聲源位于同深度的接收點A不在聲傳播的水平直達聲路徑上, 我們通過聲場計算模型計算的結果其實接近于接收點B的聲壓.

圖1 真實地球環境下的聲傳播形式Fig.1.The form of sound propagation in the real earth environment.

圖2 仿真環境下的聲傳播形式Fig.2.The form of sound propagation in the simulation environment.
針對這一失配問題, 有兩種解決思路.第一種是徐傳秀等[13]和Yan[14,15]的解決思路, 即重新構建引入地球曲率修正的二維/三維的傳播模型, 但是這種方法需要對現有的各類聲場計算模型進行重新編寫, 工程量大, 推廣起來比較困難.另一種思路是, 在不更改現有的常用聲場計算模型的情況下, 通過調整聲場計算模型的輸入參數, 使之適合實際的聲傳播情況, 從而實現地球曲率的修正.這類方法被稱為變參數媒質法[17].
如圖3所示, 假設地球為一標準球體, 地球半徑為R=6371.03km[18].A點為聲場中某個接收點, 坐標為 (X,Z) , 距離海平面的距離為z′, 聲源與接收點A之間的圓弧對應的圓心角為θ,A點距離地球球心的距離為a.與圖1類似, 圖3對海洋波導的厚度進行了放大處理.要求解的亥姆霍茲方程的形式為

其中φ為位移勢,k為介質波數,n為折射率.A點的坐標可以表示為

在圖4所示的笛卡爾坐標系中,A'點為圖3中A點映射到笛卡爾坐標系下的點, 坐標為(x,z).

圖3 極坐標系下的聲傳播形式Fig.3.The form of sound propagation in polar coordinate system.

圖4 笛卡爾坐標系下的聲傳播形式Fig.4.The form of sound propagation in Cartesian coordinate system.
設在復平面ζ和ξ上,A點、A'點的坐標可以表示為

共形映射方法可以將一些不規則或者不便用數學公式表達的區域邊界映射成規則的或已得到成熟分析的區域邊界, 是將流體力學、彈性力學、電學等學科的實際問題化繁為簡的重要方法[19].將真實的地球海洋波導中的點映射到笛卡爾坐標系所采用的映射關系為[20]

所以映射的形式為

因為

得到

因為

等價于

在共形映射中, 拉普拉斯算子的變換為[21]

在新的坐標系下, 二維的亥姆霍茲方程變為

其中

引入折射率修正參數m,

得到了映射到笛卡爾坐標系下的亥姆霍茲方程

根據折射率修正公式可以得到聲速的修正表達式為

其中c為修正后的聲速,c′為修正前的聲速.
從上面的推導可以看到, 需要調整的聲場參數有3個, 分別是從聲源到接收點的距離r, 接收點的深度z和聲速c.
對于距離r, 在真實地球環境下的r調整后, 為r=Rθ, 也就是以θ為圓心角,R為半徑的弧長.在沒有進行地球曲率修正的時候,r的計算是通過確定聲源與接收點兩點的經緯度坐標, 計算兩點間的測地線路徑的距離得到的, 修正后仍然遵循這一計算方法.
綜上, 地球曲率修正公式表示為

將深度z和聲速c修正后的表達式進行泰勒級數展開得到:

因為地球上海洋的平均海深為3734 m, 深度超過8000 m深度的海域面積只占到全部海域面積的0.04%[22].所以z′/R?1 , 對泰勒級數展開式取一階近似可得:

(18)式與Etter[16]給出的修正方法一致.運用此方法進行地球曲率修正, 計算修正后的聲場, 誤差來源有兩部分.一部分是上述一階近似引入的誤差,但因為海深遠小于地球半徑, 此誤差可以忽略.另外一部分就是環境參數修正后, 運用聲場計算模型計算時產生的誤差.
修正后的深度z是在真實海深z′的基礎上, 加上一個修正項得到的.也就是說, 修正后的深度z要略大于真實深度z′, 修正后深度相較于修正前深度的差值會隨著深度的增加逐漸變大.修正后的聲速變化情況與深度類似.
采用深海中典型Munk剖面[23]來比較修正前后的聲速剖面的變化情況.Munk剖面的表達式為

其中η=2(h-h0)/B,h0為聲道軸的深度,B為波導的垂直寬度,c0為聲速最小值(聲道軸處),ε描述了偏移極小值的量級.采用Munk給出的典型數據,B=1300m,z0=1300m,c0=1500m/s,ε=0.00737.
由圖5可以看到, 地球曲率修正后的聲速剖面在同深度的聲速值要大于未修正的聲速剖面的聲速值.在較淺的深度, 修正前后的聲速剖面差異不明顯, 在深海, 修正前后的聲速剖面有較大的差異.Munk剖面地球曲率修正前后在深海的聲速梯度對比如圖6所示, 可以看到, 修正前后深海聲速梯度變化明顯.

圖5 地球曲率修正前后Munk聲速剖面對比Fig.5.Comparison of Munk sound speed profiles before and after the earth curvature correction.

圖6 地球曲率修正前后Munk聲速剖面在深海聲速梯度對比Fig.6.Comparison of sound velocity gradient of Munk sound velocity profile before and after earth curvature correction in deep sea.
會聚區是深海中良好的聲信道, 在實際的水聲探測、通信中有較廣泛的應用.以下通過仿真研究地球曲率對會聚區傳播損失的影響.仿真的聲速剖面是上文提到的Munk聲速剖面, 采用的海底模型包含兩層液態海底, 由一層沉積層和半無限基底組成.如圖7所示, 沉積層厚度為 1 00m , 聲速為1550m/s , 密度為 1.65g/cm3, 吸收系數為 0.1dB/λ.基底中聲速為 1 800m/s , 密度 1.9g/cm3, 吸收系數為 0.1dB/λ.

圖7 仿真采用的海底模型Fig.7.Seabed model used in simulation.
圖8 給出了拋物方程近似聲場模型(RAMPE)[24]計算的引入地球曲率修正前后不同深度上的聲傳播損失(transmission loss, TL)隨距離的變化情況對比.仿真設置聲源深度為 2 00m , 中心頻率 1 00Hz , 計算的頻帶帶寬為 2 3Hz , 頻點數為11.仿真的最遠距離為 1 000km.從圖8可以看到, 修正前后都出現了明顯的會聚區傳播現象.
不同深度下的地球曲率修正前后的傳播損失對比如圖9所示, 圖9(a)對應的是會聚區傳播的情況.可以看到, 相較于未修正前, 會聚區位置前移, 并且隨著距離的增大, 前移的距離越來越大.圖9(b)和圖9(c)代表中層深度的情況, 地球曲率修正前后的傳播損失曲線沒有明顯的變化規律.圖9(d)展示了5000 m深度上的會聚區傳播中下反轉點區域的地球曲率修正前后的傳播損失變化情況, 可以看到, 進行地球曲率修正后, 傳播損失曲線變化情況與200 m深度時類似.

圖8 二維傳播損失圖 (a)地球曲率修正前; (b)地球曲率修正后Fig.8.Numerical TL: (a) Before the earth curvature correction; (b) after the earth curvature correction.

圖9 不同接收深度下的地球曲率修正前后的傳播損失對比 (a) 200 m; (b) 1000 m; (c) 2000 m; (d) 5000 mFig.9.Comparison of TL before and after correction of earth curvature at the three different receiver depths: (a) 200 m;(b) 1000 m; (c) 2000 m; (d) 5000 m.
圖10 給出了地球曲率修正前后不同深度和不同距離傳播損失的差值, 可以看到清晰的類似會聚區分布的圖樣, 體現出地球曲率導致會聚區整體向前移動的趨勢, 不同的深度上, 影響程度不同.在500—3500 m左右的深度, 修正前后傳播損失相差不大.會聚區傳播的下反轉點深度位置處, 地球曲率修正前后, 傳播損失差異較大.
王彥磊等[25]對西北太平洋的聲速剖面進行了分類, 選取文中提到的四種典型的西太平洋聲速剖面, 分別命名為Ⅰ型聲速剖面(130°E, 15°N), Ⅱ型聲速剖面(150°E, 10°N), Ⅲ型聲速剖面(135°E,25°N)和Ⅳ型聲速剖面(155°E, 35°N), 位置分布如圖11所示.聲速剖面根據WOA(world ocean atlas)[26,27]數據庫得到.

圖10 地球曲率修正前后傳播損失的差值Fig.10.The difference of the TL before and after the earth curvature correction.

圖11 西北太平洋四類典型聲速剖面的分布位置Fig.11.Distribution locations of four typical sound speed profiles in the Northwest Pacific.
圖12 給出的是西北太平洋四類典型的聲速剖面圖, 可以看到, 四類聲速聲道均為深海完全聲道,具備實現會聚區傳播的條件.
仿真設置聲源深度為200 m, 中心頻率100 Hz,頻帶帶寬23 Hz, 計算的頻點數為11.仿真的最遠距離為1000 km.海底模型如圖6所示.圖13給出的是200 m深度的接收情況.可以看到, 四類不同的聲速剖面環境下, 均出現了明顯的會聚區傳播現象.地球曲率修正后, 會聚區向聲源方向前移, 前移的幅度隨著距離的增大而逐漸增大.這一結論與Munk剖面下的仿真結果相同.
定義地球曲率修正前的各個會聚區傳播損失峰值所在的位置為會聚區的位置.針對某個特定會聚區, 先大致劃分該會聚區所在的區間, 然后對該區間內的地球曲率修正前的傳播損失曲線和地球曲率修正后的傳播損失曲線做互相關處理, 根據互相關曲線的極大值出現的位置相較于中心點的偏移, 得到一個偏移量, 將這個偏移量定義為該會聚區在進行地球曲率修正后的會聚區偏移量.經過上述方法處理, 得到了各類聲速剖面下的各個會聚區的偏移值, 如圖14所示.

圖12 西北太平洋四類典型聲速剖面Fig.12.Four types of typical sound speed profiles in the Northwest Pacific.

圖13 西北太平洋四類典型聲速剖面的200 m深度的傳播損失 (a)Ⅰ型聲速剖面的傳播損失; (b) Ⅱ型聲速剖面的傳播損失; (c) Ⅲ型聲速剖面的傳播損失; (d) Ⅳ型聲速剖面的傳播損失Fig.13.TLs at a depth of 200 m in four types of typical sound speed profiles in the Northwest Pacific: (a) TLs of type Ⅰ sound speed profile; (b) TLs of type Ⅱ sound speed profile; (c) TLs of type Ⅲ sound speed profile; (d) TLs of type Ⅳ sound speed profile.

圖14 不同的聲速剖面下計算的會聚區移動情況Fig.14.Movement of the convergence zone calculated under different sound speed profiles.
由圖14可以看到, 不同聲速剖面下, 相同距離上, 地球曲率修正后的會聚區偏移程度相近, 傳播200 km距離會聚區偏移2 km左右, 傳播500 km距離會聚區偏移5 km左右, 傳播1000 km的距離會聚區偏移可達10 km.對某特定的聲速剖面, 會聚區偏移隨距離變化呈現近似的線性關系, 隨距離變大, 地球曲率會導致會聚區偏移程度越來越大, 傳播距離每增長100 km, 會聚區偏移增加約1 km.
深海聲道, 又被稱為SOFAR聲道, 在大洋層析、遠程水聲導航等領域中有重要應用.仿真設置的海洋環境為前文所述的Munk剖面, 聲源深度為1300 m, 仿真時采用的聲源脈沖S(t)[28]為

這是Hanning加權的四周期正弦波.聲源的中心頻率fc為200 Hz, 脈沖的時間窗長度為10 s.以下通過仿真分析了不同距離上的脈沖到達結構.
圖15給出了歸一化后1000 km距離處地球曲率修正前后的到達結構, 圖中不同顏色體現了到達脈沖的相對強弱.可以看到, 地球曲率對到達結構產生的主要影響有三點.1)地球曲率使得到達結構在時間軸上整體前移, 前移的幅度約為136 ms.這是因為地球曲率修正后, 聲速剖面的聲速略大于修正前的聲速導致的.2)受地球曲率影響, 地球曲率修正后得到的到達結構相較于修正前, 出現了整體的擴展, 包括深度方向上的延拓和時間方向上的展寬, 越靠近到達結構左側, 展寬程度越大.這是也是因為修正后聲速增大, 導致部分聲線路徑的傳播時間變短, 有些路徑的聲線會比之前先到達接收點.3)地球曲率使得各個深度上的多途結構的能量分布也發生了變化.

圖15 1000 km距離上的到達結構 (a)地球曲率修正前;(b)地球曲率修正后Fig.15.Arrival structure at 1000 km: (a) Before the earth curvature is corrected; (b) after the earth curvature is corrected.
圖16 比較了500, 1300和2500 m三個深度上歸一化處理后的地球曲率修正前后的接收時域波形.如圖16(b)所示, 地球曲率修正前后, 接收的時域波形明顯前移, 前移幅度約為140 ms.并且由于地球曲率修正后到達結構整體的展寬以及能量的變化, 導致部分深度到達的時域波形產生明顯的差異, 如圖16(a)和圖16(c)所示.
圖17和圖18比較500和2000 km距離上的地球曲率修正前后的到達結構.從不同距離的到達結構的比較可以看到, 1000 km距離上地球曲率對到達結構產生的三類主要影響在各個距離上都有體現, 并且影響程度都隨著傳播距離的增大逐漸增大.

圖16 地球曲率修正前后1000 km距離上的不同接收深度的時域到達波形比較 (a) 500 m; (b) 1300 m; (c) 2500 mFig.16.Comparison of time-domain arrival waveforms at 1000 km before and after the earth curvature correction at different receiver depths: (a) 500 m; (b) 1300 m;(c) 2500 m.
以上的分析都是在將地球近似為標準球體時進行的.而實際上, 地球是一個扁平的橢球體, 如圖19所示.國際上對地球橢球建立了不同的數學模型[18], 這些數學模型對應不同的地球橢球參數,比較著名的有貝塞爾橢球、克拉克橢球、克拉索夫斯基橢球、以及全球定位系統(GPS)所使用的WGS84橢球, WGS84橢球也是最常用的地球橢球.
WGS84橢球對應的地球橢球參數為[18]: 子午橢圓的長半軸a=6378.137km ; 子午橢圓的短半軸b=6356.752km ; 子午橢圓的扁率α=(a-b)/a=1/298.2572; 子午橢圓的第一偏心率子午橢圓的第一偏心率

圖17 500 km距離上的到達結構 (a)地球曲率修正前; (b)地球曲率修正后Fig.17.Arrival structure at 500 km: (a) Before the earth curvature is corrected; (b) after the earth curvature is corrected.

圖18 2000 km距離上的到達結構 (a)地球曲率修正前; (b)地球曲率修正后Fig.18.Arrival structure at 2000 km: (a) Before the earth curvature is corrected; (b) after the earth curvature is corrected.

圖19 地球橢球模型Fig.19.Earth ellipsoid model.
在WGS84橢球下, 有幾種常用的地球曲率半徑, 包括子午線曲率半徑、卯酉線曲率半徑、平均曲率半徑和任意方向上的曲率半徑.
子午圈曲率半徑定義為子午圈上一點的曲率半徑, 設為M,

卯酉圈是與子午面垂直的閉合圓圈, 卯酉圈的曲率半徑用N來表示:

平均曲率半徑Ra是指經過曲面任意一點所有可能方向上的法截線曲率半徑RA的算數平均值.推導得到曲面上任意一點的平均曲率半徑表示為

圖20給出了從赤道到極地的地球平均曲率半徑變化情況.赤道上地球平均曲率半徑最小,R0為6356.863 km, 隨著緯度的升高, 地球平均曲率半徑增大, 到極點上平均曲率半徑取最大值,R90為6399.699 km.

圖20 地球平均曲率半徑隨緯度的變化情況Fig.20.The variation of the average radius of curvature of the earth with latitude.
下面分析不同地球半徑取值對地球曲率影響下環境參數修正方法的影響.設ca為地球平均半徑下經過地球曲率修正后的Munk剖面的聲速值,c0為赤道位置上地球平均曲率半徑下經過地球曲率修正后的Munk剖面的聲速值,c90為極點位置上地球平均曲率半徑下經過地球曲率修正后的Munk剖面的聲速值,
圖21給出了由于地球半徑選取造成的修正后的聲速剖面的差異.可以看到, 三種聲速剖面幾乎完全重合.5000 m深度上, 三類聲速剖面之間的差值在 1 0-3m/s量級, 在實際計算中, 此差異造成的影響可以忽略.所以將地球近似為一標準球體進行地球曲率修正, 完全可以滿足修正后計算精度的需要.

圖21 選取不同地球半徑計算聲速剖面的比較Fig.21.Comparison of sound speed profiles calculated by selecting different earth radius.
本文提出了一種計算簡便的地球曲率影響下環境參數的修正方法, 并對地球曲率對聲傳播的影響進行了分析, 主要結論如下.
1)對于會聚區傳播, 地球曲率會造成會聚區位置相對于未修正前的計算位置前移, 偏移程度隨著距離的增大而增大, 1000 km距離上偏移可達10 km.通過對不同的聲速剖面仿真結果分析得到,會聚區移動隨距離變化近似為線性關系, 傳播距離每增大100 km, 會聚區的偏移增加約1 km.
2)對于深海聲道傳播, 地球曲率會導致到達結構整體向前移動, 前移的幅度隨著距離的增大而逐漸增大, 1000 km距離上整體前移的幅度可達136 ms.地球曲率還會造成整個到達結構在深度和時間方向上的擴展.在特定深度上接收到的時域波形在地球曲率修正前后有較大差異.
3)比較不同地球近似模型下的修正效果表明,將地球近似為標準球體即可滿足計算精度需要.
以上工作是以水平不變的環境為前提進行分析的, 以后將進一步研究水平變化海洋環境下的地球曲率對聲傳播以及遠程信息傳輸的影響.