王子蘅， 王振杰， ， 聶志喜， ， 張遠帆

(1. 中國石油大學(xué)(華東) 海洋與空間信息學(xué)院， 山東 青島， 266580； 2. 青島海洋科學(xué)與技術(shù)試點國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室， 山東 青島， 266071)

聲波在水中的良好傳播特性使其成為獲取和傳遞水下信息最為有效的手段[1-3]。聲速是影響水下多波束系統(tǒng)作業(yè)精度的重要外部影響因素， 通過影響聲線跟蹤的精度， 最終影響到測深精度。多波束測深通過換能器實時接收其發(fā)射出的各波束經(jīng)海底反射和散射后返回的到達角和旅行時[3]， 利用聲速剖面數(shù)據(jù)，由公式計算得到不同波束點對應(yīng)的水深值。聲速剖面是聲速的垂直結(jié)構(gòu)分布， 海水的介質(zhì)特性導(dǎo)致聲波的傳播軌跡發(fā)生彎曲， 要獲取波束腳印的確切位置， 需要沿著波束的傳播路徑追蹤聲線， 計算波束腳印的水平位移和深度， 即聲線跟蹤[3]。因此， 聲速剖面的正確與否直接影響多波束測深結(jié)果的精度和可靠性[4]。

海洋中的聲速剖面通常由以下兩種方法獲得：直接測量法和間接測量法[4]。直接測量法通過測量聲波在海水中已知的固定距離內(nèi)往返所需的時間或相位來計算聲速值。間接測量法利用溫鹽深測量系統(tǒng)(conductivity， temperature， depth， CTD) 測量海水的溫度、鹽度及深度等物理量， 再通過聲速經(jīng)驗公式計算作業(yè)水域內(nèi)各采樣深度的聲速值， 進而構(gòu)成聲速剖面。然而在深水大洋進行實際作業(yè)時， 受到水深和洋流的影響， 停船投放CTD采集聲速剖面數(shù)據(jù)時間成本較高。采用走航式的拋棄式聲速剖面計雖然可以快速有效的獲得聲速剖面， 但其作業(yè)成本隨之大大增加。相較于以上方法， 聲速剖面重構(gòu)只需少量的實測聲速數(shù)據(jù)即可實現(xiàn)全海深聲速剖面的重構(gòu)， 提高效率的同時還有效的降低了成本。

由于海洋環(huán)境的復(fù)雜性， 聲速剖面難以用一個簡單的函數(shù)表示。Davis等[5]證明EOF是在最小均方差意義下反映聲速剖面最有效的基函數(shù)； 孫文川等[6]論證了EOF重構(gòu)的聲速剖面具有較高的內(nèi)符合精度， 能較好地描述實際聲速剖面； 張志偉等[7]研究表明取前6階及以上的EOF表示聲速剖面可滿足多波束測深的精度要求； 張鎮(zhèn)邁等[8]采用部分深度的實測剖面， 并基于測區(qū)海域內(nèi)已有的數(shù)據(jù)， 對全海深聲速剖面實現(xiàn)了重構(gòu)；張孝首等[9]利用實測CTD數(shù)據(jù)通過有理擬合延伸和EOF方法實現(xiàn)了深海5 000 m聲速剖面場的重構(gòu)； 李洪超等[10]對不同深度位置的部分聲速剖面的重構(gòu)效果進行了探討； 張維等[11]在對殘缺樣本聲速合理外延的基礎(chǔ)上， 選取聲速剖面變化較劇烈深度的3個聲速值實現(xiàn)了剖面重構(gòu)。以上研究表明， 基于EOF方法利用部分實測聲速數(shù)據(jù)重構(gòu)全海深剖面是可行有效的， 但未對參與重構(gòu)的實測數(shù)據(jù)的采樣深度給出明確依據(jù)。

本文將測區(qū)Argo浮標的溫鹽深數(shù)據(jù)通過聲速經(jīng)驗公式計算得到的聲速剖面作為測區(qū)歷史聲速剖面資料， 基于EOF空間函數(shù)的方差貢獻率選取K(EOF階次)個最大梯度對應(yīng)深度處的聲速值進行聲速剖面重構(gòu)， 在保證多波束測深精度的前提下實現(xiàn)了對全海深聲速剖面的重構(gòu)。

1 聲速剖面EOF表示原理

1.1 聲速剖面EOF分解

經(jīng)驗正交函數(shù)分析方法也稱特征向量分析或主成分分析， 是分析矩陣數(shù)據(jù)的特性、提取主要數(shù)據(jù)特征量的一種方法[12]。

在某海域采集N個聲速剖面：C= [c(1)，c(2)，…，c(N)]， 內(nèi)插成為M個垂直標準層， 得到聲速矩陣CM×N：

式中， 每一列為一個聲速剖面的標準層插值， 每一行為所有剖面在同一深度的聲速。將N條聲速剖面

擾動矩陣的協(xié)方差矩陣為：

其中每個元素rij為：

將RM×M特征分解， 有：

式中，DM×M為特征值矩陣， 令特征值λi按從大到小順序排列：FM×N為特征值對應(yīng)的特征向量矩陣， 也就是EOF空間函數(shù)， 可表示為：

式中，f(k)(k= 1，2，…，N)為M×1的列向量， 即為確定的第k階EOF。特征值λi對應(yīng)的分量f(i)即為第i階EOF。特征值的大小代表了其對應(yīng)EOF空間函數(shù)對線性空間的影響權(quán)重[7]， 即越大的特征值對應(yīng)的特征向量中包含的聲速重構(gòu)信息越豐富。

1.2 聲速剖面EOF重構(gòu)

測區(qū)內(nèi)的聲速剖面完成EOF分解后， 利用前幾階EOF即可完成測區(qū)內(nèi)任一聲速剖面的重構(gòu)：

式中，c(z)為重構(gòu)的聲速剖面；z為各層海水深度；為平均聲速；K為EOF階次；αi為重構(gòu)系數(shù)；fi(z)為EOF空間函數(shù)。移項得：

則重構(gòu)系數(shù)矩陣AK×N為：

式中，ai= [ai(1)，ai(2)，…，ai(K)]T為每個聲速剖面對應(yīng)的重構(gòu)系數(shù)。得到重構(gòu)系數(shù)后， 即可按式(9)對測區(qū)內(nèi)任一聲速剖面進行EOF重構(gòu)。

2 聲速剖面EOF重構(gòu)的數(shù)據(jù)采樣深度選取方法

本文給出了用于聲速剖面EOF重構(gòu)的數(shù)據(jù)采樣深度選取的方法步驟。算法流程圖見圖1， 詳細步驟和說明如下：

(1) 從中國Argo實時資料中心(http：//www.argo.org.cn/)獲取某海域溫鹽壓數(shù)據(jù)， 預(yù)處理后利用聲速經(jīng)驗公式計算得到采樣位置處的聲速剖面。目前國內(nèi)外較為認可的聲速經(jīng)驗公式有Chen-Millero，Wilson， Del Grosso及Leroy等， 各模型的溫度、鹽度、壓力的適用范圍有所差異。根據(jù)文獻[13]的研究， Del Grosso， Chen-Millero， C.C.Leroy、Coppens四個經(jīng)驗聲速模型在全球海域具有較高的精度及適用性， 文獻[14]指出， 在350～1 000 m的深度上， Chen-Millero公式的誤差最小， 因此本文選用Chen-Millero 聲速經(jīng)驗公式計算的聲速剖面作為樣本聲速剖面， 利用Akima插值[15]將聲速剖面數(shù)據(jù)內(nèi)插到垂直標準層。

(2) 按照1.1節(jié)中的方法進行EOF分解， 提取EOF函數(shù)， 計算各階EOF的方差貢獻率和累計方差貢獻率。方差貢獻率[16]指的是每一階EOF包含聲速場信息的百分比。其計算公式為：

式中，σi表示第i階方差貢獻率，λi表示協(xié)方差矩陣R經(jīng)降序排列后的第i個特征值。

前K階EOF對區(qū)域的累計貢獻率可表示為：

一般地， 如果w(K)的值大于等于95%， 則認為利用前K階EOF能較好的表示該區(qū)域的主要信息[5]。

(3) 計算測區(qū)內(nèi)每一聲速剖面每一深度層的聲速梯度， 然后計算平均分層聲速梯度：

式中，ci，n和ci+1，n分別表示第n條聲速剖面第i層和第i+1層的聲速值，zi+1和zi+1，n分別表示第n條聲速剖面第i層和第i+1層的深度值。按照平均分層聲速梯度由大到小， 對聲速層進行降序排列， 形成新的聲速層矩陣。

(4) 從新構(gòu)建的聲速層矩陣中選取K個深度處的聲速值(即選取K個聲速變化最劇烈處的聲速值)，結(jié)合K階EOF計算得到重構(gòu)系數(shù)矩陣AK×N， 重構(gòu)全海深聲速剖面。

(5) 基于實測聲速剖面， 預(yù)設(shè)水深值H和初始入射角， 根據(jù)分層常梯度聲線跟蹤法[2]， 以H作為約束迭代計算聲波的傳播時間T：

式中，Y表示實測聲速剖面的層數(shù)，θi和θi+1分別表示實測剖面第i層和第i+1層的入射角，gi表示實測剖面第i層的梯度，ci和ci+1分別表示實測剖面第i層和第i+1層的聲速， Δzi表示第i層的深度差。p為Snell系數(shù)：

(6) 用重構(gòu)的聲速剖面和T反算波束水深值[6]H′。聲波在重構(gòu)的聲速剖面第i層的傳播時間ti為：

式中，iθ′和θi1+′分別表示重構(gòu)聲速剖面第i層和第i+1層的入射角，gi′表示重構(gòu)剖面第i層的梯度，ic′和ci1+′分別表示重構(gòu)剖面第i層和第i+1層的聲速。

式中，J表示使等式成立的ti的個數(shù)。所以：

(7) 計算測深誤差σ：

定義有效波束比為滿足深度限差的波束占總波束的比率[17]， 首先統(tǒng)計滿足0.25%水深限差要求的波束數(shù)數(shù)， 然后判斷多波束測深的有效波束比是否達到100%， 如果是， 則輸出測深偏差， 實驗結(jié)束，反之， 則將選取階數(shù)K加1， 重新完成重構(gòu)系數(shù)矩陣的計算以及之后的步驟， 直至符合滿足測深精度的要求。

3 實驗分析

實驗數(shù)據(jù)源于中國Argo資料實時中心發(fā)布的《全球Argo浮標剖面觀測資料質(zhì)量再控制數(shù)據(jù)集》[18]， 選取2006—2009年每年1月份的Argo浮標數(shù)據(jù)。將153.5°W—155W°， 27.5°N—29°N的海區(qū)作為區(qū)域1；153.5°W—155.5W°， 22.5°N—24.5°N的海區(qū)作為區(qū)域2。去除不合格的剖面后， 區(qū)域1選取了18個剖面， 區(qū)域2選取了15個剖面。

典型深海聲速剖面結(jié)構(gòu)如圖2所示[2]， 主要劃分為三部分： 混合層(由表面層和季節(jié)躍變層構(gòu)成)、主躍層和深海等溫層。其中， 聲速在混合層和主躍層的變化情況復(fù)雜， 而在深海等溫層中， 聲速隨深度的增加呈趨勢穩(wěn)定的正梯度增加， 接近線性變化。因此本文重點對混合層和主躍層這兩部分聲速變化情況更為復(fù)雜的水層的聲速剖面進行研究， 統(tǒng)一截取海面至海深1 000 m的深度之間作為本次實驗的剖面范圍。

圖2 深海典型聲速剖面圖Fig. 2 Typical deep-sea sound speed profile

隨機選取區(qū)域內(nèi)1條剖面用作檢核重構(gòu)精度，其余剖面作為已知剖面數(shù)據(jù)。實驗區(qū)域內(nèi)的Argo剖面浮標分布如圖3(a)和(b)所示， 圖中“×”點代表已知的樣本剖面，“▲”代表用于檢核的剖面。

圖3 實驗海域內(nèi)浮標剖面分布圖Fig. 3 Distribution map of the buoy profiles in the experimental sea area

3.1 數(shù)據(jù)預(yù)處理

(1) Argo剖面浮標測量的是海水的溫度、鹽度和壓力數(shù)據(jù)， 采用 Chen-Millero聲速經(jīng)驗公式[19]計算得到相應(yīng)深度處的聲速值：

式中，

其中，T為溫度， ℃；S為鹽度；P為壓力， bar。

(2) 將計算得到的聲速數(shù)據(jù)利用Akima插值函數(shù)內(nèi)插成1 m的等間距標準層。樣本聲速剖面結(jié)構(gòu)如圖4。

圖4 聲速剖面結(jié)構(gòu)示意圖Fig. 4 Schematic diagrams of SSP

3.2 確定參與重構(gòu)的數(shù)據(jù)深度及聲速值

分別對兩個區(qū)域的樣本聲速剖面進行EOF分解，得到測區(qū)EOF空間函數(shù)的方差貢獻率和累計方差貢獻率， 如表1所示?？梢钥闯?， 在區(qū)域1中， 前5階EOF的累計方差貢獻率可達97.985%， 在區(qū)域2中，前5階EOF的累計方差貢獻率達95.813%， 均已大于95%， 因此利用前5階EOF能夠較好的表示當前區(qū)域的主要信息。

表1 聲速剖面不同階次EOF方差貢獻率Tab. 1 Variance contribution rate of the different order EOF of the SSP

圖5 給出的是由兩個實驗海域內(nèi)的樣本聲速剖面提取的1～5階EOF， 可以明顯看出， 第1階EOF空間函數(shù)的變化較小， 整體呈趨勢穩(wěn)定， 這表明前幾階EOF中包含了實驗海域內(nèi)的聲速場的主要信息， 而2～5階EOF空間函數(shù)的變化明顯， 變化集中在淺層區(qū)域， 尤其是在越靠后的階次中， 這反映了聲速場相對于平均SSP的細節(jié)變化[7]。

圖5 1～5階EOF空間函數(shù)Fig. 5 1～5 order EOFs

計算每個樣本聲速剖面每層的梯度變化情況，得到該區(qū)域的平均分層聲速梯度及對應(yīng)的深度位置， 按照平均分層聲速梯度由大到小對聲速層降序排列， 自上而下選取5個深度處的聲速值， 選取結(jié)果如表2所示。

表2 選取的深度及聲速Tab. 2 Selected depths and sound speed values

3.3 聲速剖面重構(gòu)

選取表2提取的聲速值計算重構(gòu)系數(shù)矩陣， 進而重構(gòu)全海深聲速剖面。兩個實驗海域內(nèi)聲速剖面的重構(gòu)結(jié)果如圖6所示， 重構(gòu)誤差如圖7所示， 區(qū)域1的均方根誤差為0.514 m/s， 區(qū)域2的均方根誤差為0.609 m/s?？梢钥闯觯?重構(gòu)聲速剖面與實測聲速剖面整體符合程度較好； 在兩個區(qū)域的淺層部分， 尤其是區(qū)域1的100 m深度附近和區(qū)域2的350 m深度附近，重構(gòu)聲速剖面與實測聲速剖面有明顯出入， 這可能是由于只選取了前幾階的EOF空間函數(shù)進行重構(gòu)， 部分階次的EOF空間函數(shù)并不能包含測區(qū)聲速場的全部信息， 造成了部分信息的損失， 因而出現(xiàn)了一定的偏差。此外， 以500 m水深作為分界， 深層部分的聲速剖面重構(gòu)結(jié)果較為穩(wěn)定， 整體的重構(gòu)誤差較??； 淺層部分的聲速剖面重構(gòu)結(jié)果相比較差， 造成該現(xiàn)象的原因可能是淺層部分的海水包含混合層和主躍層， 海洋要素在其中的變化情況復(fù)雜， 而在深層海水中已趨于穩(wěn)定。

圖6 重構(gòu)的聲速剖面Fig. 6 Reconstructed SSP

圖7 聲速剖面重構(gòu)誤差Fig. 7 Error of reconstructed SSP

3.4 多波束測深檢驗

在多波束測深中， 利用EOF重構(gòu)的聲速剖面最終都要應(yīng)用于聲線改正。因此對于利用重構(gòu)的聲速剖面進行聲線改正能否滿足實際應(yīng)用的要求， 還需要作進一步的檢驗。美國國家海洋與大氣局(national oceanic and atmospheric administration， NOAA)對聲速剖面引起的水深誤差規(guī)定[20]： “若使用的聲速剖面與實際聲速剖面對測深改正造成的互差超過0.25%水深則視為超限”。本文實驗以NOAA的測深限差作為檢驗聲速剖面重構(gòu)精度的標準， 并結(jié)合有效波束比來判斷利用EOF重構(gòu)聲速剖面的有效性。設(shè)置水深和波束的初始入射角， 基于實測聲速剖面， 分層常梯度聲線跟蹤波束的往返時間， 利用該時間與重構(gòu)的聲速剖面再通過分層常梯度聲線跟蹤法反算水深， 并與真實水深值進行對比， 比較兩者差值是否滿足0.25%水深限差， 最后統(tǒng)計有效波束比。

表3 為使用兩則區(qū)域重構(gòu)的聲速剖面進行多波束測深的結(jié)果， 可以看出， 使用重構(gòu)的聲速剖面進行聲速改正得到的各波束點水深均符合限差要求，有效波束比均達到了100%。

表3 重構(gòu)聲速剖面的最大偏差和有效波束比Tab. 3 Maximum deviations and effective beam ratios of the reconstructed SSPs

4 結(jié)論

本文對重構(gòu)全海深聲速剖面所需實測數(shù)據(jù)的采樣深度進行了研究， 實驗證明在對歷史聲速剖面資料進行分析后， 按照K階EOF函數(shù)方差累計貢獻率大于95%， 選取K個最大平均分層聲速梯度處的聲速值作為已知值， 重構(gòu)的聲速剖面能夠滿足NOAA對多波束測量聲速剖面造成的誤差要求， 為實際測量作業(yè)中聲剖數(shù)據(jù)的采樣深度提供了參考， 提高了工作效率。由于海洋聲速易受各種環(huán)境復(fù)雜性的影響， 因此其他海區(qū)的選取階次及深度會有所不同， 但本文介紹的方法和所得規(guī)律性結(jié)論具有一定的普適性。