李林洋, 徐天河, 王君婷, 黃威, 高凡, 舒建旭

(1.長安大學 地質(zhì)工程與測繪學院, 陜西 西安 710054; 2.山東大學 空間科學與物理學院, 山東 威海 264209; 3.中國海洋大學 信息科學與工程學部, 山東 青島 266100)

21世紀海洋經(jīng)濟的重心下移,近年來全球重大油氣發(fā)現(xiàn)有70%來自水深超過1 000 m的水域,并且呈現(xiàn)逐年升高的趨勢[1]。對于固定資源的深海采集工作高精度的水下導(dǎo)航、定位必不可少,而聲速誤差作為高精度水下導(dǎo)航、定位的主要誤差源之一,構(gòu)建高精度的聲速場有著重要意義。水下聲速剖面可以由聲速剖面儀(sound velocity profiler, SVP)直接測量得到,或是通過溫鹽深剖面儀(temperature, conductivity, depth profiler, CTD)以及拋棄式溫鹽深剖面儀(expendable CTD, XCTD)經(jīng)過聲速經(jīng)驗公式獲得。直接測量法獲得聲速剖面的優(yōu)勢是結(jié)果準確,但是采集效率低、成本高,無法進行大規(guī)模的聲速剖面測量。

利用歷史數(shù)據(jù)反演得到的聲速剖面是海區(qū)歷史或者平均聲速剖面,不能據(jù)此準確地預(yù)測實際的聲速剖面[2-3]。據(jù)此,周士弘等[4]有提出了至少利用5階經(jīng)驗正交函數(shù)描述聲速場,進而建立模型描述海區(qū)聲速剖面,探討了海面溫度衛(wèi)星遙感或航空遙感技術(shù)、海水溫躍層水溫預(yù)報技術(shù)同聲速場經(jīng)驗正交函數(shù)描述方法結(jié)合在一起來預(yù)報海表以下深度的聲速剖面。艾銳峰[5]利用后向傳播(back propagation, BP)神經(jīng)網(wǎng)絡(luò)的自適應(yīng)能力和自學習能力,以歷史數(shù)據(jù)為樣本對模型進行訓(xùn)練,得到反演預(yù)測模型進行聲速剖面反演預(yù)測。禹小康[2]利用遙感數(shù)據(jù)提供的表面溫度和神經(jīng)網(wǎng)絡(luò)進行聲速剖面反演預(yù)測。但是使用神經(jīng)網(wǎng)絡(luò)模型構(gòu)建需要大量的歷史經(jīng)驗聲速剖面作為參考樣本,而對于局部海域經(jīng)常面臨參考樣本稀少,容易導(dǎo)致模型過擬合而精度性能降低。

海洋聲層析最初是由Munk和Wunsch提出的,在射線理論的基礎(chǔ)上,通過匹配聲場觀測估計海洋環(huán)境參數(shù)[6-8]。Tolstoy[9]提出用匹配場方法(match field process,MFP)反演聲速剖面。為了加快匹配場的反演過程,引入啟發(fā)式算法,如粒子群算法(particle swarm optimization,PSO)[10-11]、模擬退火算法(simulated annealing, SA)[12]、遺傳算法(genetic algorithm,GA)[3,13]。匹配場處理方法回避了建立聲場信息到聲速分布的反向映射關(guān)系的難題,為聲速場構(gòu)建提供了有效實現(xiàn)途徑[14]。

針對上述方法所存在的問題,本文提出了聯(lián)合匹配場處理和神經(jīng)網(wǎng)絡(luò)的聲速時間場構(gòu)建方法。母船走航或浮標采集的通信數(shù)據(jù)相較于直接測量聲速剖面是便捷、經(jīng)濟的,而且其采集的通信數(shù)據(jù)較為完整且連續(xù)。首先該算法通過通信數(shù)據(jù)依照匹配場反演算法進行聲速剖面仿真,加密參考樣本。然后通過對非線性函數(shù)有較好擬合效果的BP神經(jīng)網(wǎng)絡(luò)進行EOF重構(gòu)系數(shù)的擬合,形成時間(本文所使用的時間是距離2023年3月的秒計時)到EOF重構(gòu)系數(shù)的非線性映射關(guān)系。最后根據(jù)BP神經(jīng)網(wǎng)絡(luò)構(gòu)建出的模型輸入XCTD時間生成聲速剖面與實測XCTD結(jié)果進行對比。

1 MFPBP聲速時間場構(gòu)建

1.1 實驗原理

1)EOF分解。

由于經(jīng)驗正交函數(shù)分解(orthogonal empirical function, EOF)能夠顯著降低反演參數(shù)數(shù)量,這對于匹配場處理有極其重要的作用,使用CTD與SVP的14條全水深聲速剖面當作參考樣本進行EOF分解,將分解得到的特征向量當作這片海區(qū)的基函數(shù)。通過反演EOF重構(gòu)系數(shù)來生成聲速剖面,經(jīng)驗正交函數(shù)求解過程如下。

將實測CTD與SVP的14條全水深聲速剖面組合成矩陣的形式:

(1)

(2)

計算協(xié)方差矩陣R的特征向量VN×N和特征根λ(1,…,n),兩者滿足:

RN×N×VN×N=VN×N×ΛN×N

(3)

其中ΛN×N為:

(4)

對角線元素為特征根,將特征根從大到小排列,每一個特征根對應(yīng)一個特征向量,也稱EOF函數(shù)。取前k階的正交函數(shù),則任意聲速剖面可以表示為:

(5)

式中:ai稱為EOF重構(gòu)系數(shù);fi(z)為經(jīng)驗正交函數(shù)。將EOF投影到原始聲速剖面數(shù)據(jù)上,即可利用最小二乘原理計算各聲速剖面對應(yīng)的EOF重構(gòu)系數(shù):

A=(VTV)-1V·ΔC

(6)

式中:A為EOF重構(gòu)系數(shù);ΔC為實際與均值的偏差矩陣。

2)匹配場處理。

由于經(jīng)過EOF分解,任意聲速剖面可以表示為:

匹配場算法實質(zhì)是利用啟發(fā)式算法擬合EOF重構(gòu)系數(shù),利用射線追蹤原理模擬計算理論聲場信息,通過理論聲場信息與實測聲場信息進行匹配用于指導(dǎo)啟發(fā)式算法。

射線追蹤理論在分層線性聲速剖面近似模型下能夠建立聲速剖面分布到聲場信息的映射關(guān)系,對于某一給定的聲速剖面c(z),文獻[15]給出了線性梯度聲速層中信號傳播時間、水平傳播距離關(guān)于掠射角的計算公式:

(7)

(8)

式中:Δtj,i表示由信號源向第j個接收節(jié)點發(fā)送的信號在第i個深度層的傳播時間;Δhj,i是相應(yīng)水平傳播距離;gi是第i個深度層的聲速梯度;θj,i是由信號源向第j個接收節(jié)點發(fā)送的信號在第i個聲速層的信號掠射角。當i=0時,θj,0是信號的初始掠射角,而c0是聲速剖面的初始聲速值。分層線性聲速剖面近似模型中相鄰聲速剖面采樣點之間聲速變化是線性的,即:

ci=ci-1+giΔdi,i=1,2,…,I

式中Δdi是第i個聲速采樣點和第i+1個聲速采樣點之間的深度差。

上式信號傳播時間和水平傳播距離可以通過斯涅爾折射定律推導(dǎo)成與初始掠射角的關(guān)系。掠射角斯涅爾定律為:

c0cosθj,i=cicosθj,0,i=1,2,…,I

(9)

式中:ci為聲速;θj,i為第j組收發(fā)節(jié)點的時延信息。在第i個深度層的掠射角,通過三角函數(shù)關(guān)系可得:

(10)

(11)

(12)

當使用啟發(fā)式算法生成聲速剖面時,可以通過對初始掠射角設(shè)定45°通過二分法由式(11)計算直達聲所對應(yīng)的h,再通過計算h和實際測量設(shè)置閾值來求得初始掠射角。通過初始掠射角求得時延t,通過式(12)計算損失函數(shù)來指導(dǎo)啟發(fā)式算法,式中tcal為通信數(shù)據(jù)的測量值,即時延信息;tj為匹配場生成的時延;N為單個時段內(nèi)通信數(shù)據(jù)的個數(shù)。

本文匹配場算法使用的啟發(fā)式算法為遺傳算法(genetic algorithm,GA)。遺傳算法是借鑒生物的自然選擇和遺傳進化機制而開發(fā)出的一種全局優(yōu)化自適應(yīng)概率搜索算法[16]。能夠提高匹配場收斂速度,提高計算效率。

其運算過程如下。

①初始化。采用二進制編碼將EOF重構(gòu)系數(shù)編碼為基因型X,個體的表現(xiàn)型x和基因型X之間可以通過編碼解碼相互轉(zhuǎn)化,設(shè)置進化代數(shù)計數(shù)器t←0;設(shè)置最大進化代數(shù)T;隨機生成M個個體作為初始群體P(0)。

②個體評價。計算群體P(t)中各個個體的適應(yīng)度,以式(12)作為適應(yīng)度函數(shù)。

③選擇運算。先計算出群體中所有個體的適應(yīng)度綜合∑fi;其次計算出每個個體的相對適應(yīng)度大小fi/∑fi;最后產(chǎn)生0到1之間的隨機數(shù),通過隨機數(shù)出現(xiàn)在哪一個概率區(qū)域來確定各個個體被選中的次數(shù)。

④交叉運算。先對群體進行隨機配對;其次隨機設(shè)置交叉點位置;最后相互交換配對染色體的部分基因。

⑤變異運算。首先確定各個個體的基因編譯位置;然后依照某一概率將變異點原有基因值取反。群體P(t)經(jīng)過選擇、交叉、變異運算后得到下一代群體P(t+1)。

⑥終止條件判斷。若t≤T,則:t←t+1,轉(zhuǎn)到②;若t>T,則以進化過程中所得到的足有最大適應(yīng)度的個體作為最優(yōu)解輸出,終止計算。

3)神經(jīng)網(wǎng)絡(luò)。

BP神經(jīng)網(wǎng)絡(luò)具有高度非線性映射能力,對于海洋復(fù)雜多變的環(huán)境引起的聲速剖面的變化具有很強的適應(yīng)性。圖1為一個3層BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意。

圖1 神經(jīng)網(wǎng)絡(luò)示意Fig.1 Schematic of neural network

其中,輸入層向量為X=[x1x2…xI]T,共計I個輸入分量。隱藏層包括P個神經(jīng)元。輸出層向量為Y=[y1y2…yO]T,共計O個輸出分量。輸入層到隱藏層的權(quán)重向量為V=[vip],i=1,2,…,I,p=1,2,…,P,閥值為A=[a1a2…ap]T,其中ap=x0vop;激活函數(shù)為f1。隱藏層到輸出層的權(quán)重向量為W=[wpo],p=1,2,…,P,o=1,2,…,O,閥值為B=[b1b2…bO]T,其中bo=h0w0o;激活函數(shù)為f2。于是:

(13)

(14)

1.2 MFPBP聲速時間場構(gòu)建確定步驟

本文給出了MFPBP聲速時間場構(gòu)建確定方法。算法流程見圖2,詳細流程如下。

圖2 MFPBP聲速時間場構(gòu)建確定流程Fig.2 MFPBP sound velocity time field construction determination process

根據(jù)實測的數(shù)據(jù)進行聲速剖面預(yù)處理,得到參考樣本CTD與SVP的聲速剖面14條和驗證樣本XCTD的聲速剖面9條。

在通信數(shù)據(jù)預(yù)處理環(huán)節(jié),通過母船換能器與海底應(yīng)答器的距離除以大致聲速1 500 m/s與實際測量的時延數(shù)據(jù)作對比,進行通信數(shù)據(jù)的選取。

對參考樣本進行EOF分解得到參考樣本的特征向量和EOF重構(gòu)系數(shù),取前5階特征向量當作這片海區(qū)的基函數(shù)。

本文所使用的匹配場算法主要包含4個部分,通過EOF分解獲得參考樣本的特征向量;搜索主成分系數(shù)(EOF重構(gòu)系數(shù))并生成模擬聲速剖面;根據(jù)射線追蹤理論模擬計算理論聲場信息;理論聲場信息與實測聲場信息匹配并將匹配時的模擬聲速剖面確立為聲速剖面反演結(jié)果。

匹配場反演出的EOF重構(gòu)系數(shù)和參考樣本的EOF重構(gòu)系數(shù)作為神經(jīng)網(wǎng)絡(luò)的訓(xùn)練樣本。以反演聲速剖面所處時間(本文所使用的時間是距離2023年3月的秒計時)作為輸入,EOF重構(gòu)系數(shù)作為輸出。隱含層節(jié)點數(shù)采用經(jīng)驗公式與試湊法確定:

式中:P為隱藏層節(jié)點個數(shù);I為輸入層節(jié)點個數(shù);O為輸出層節(jié)點個數(shù)a=[0,10]。

根據(jù)匹配場與參考樣本的EOF重構(gòu)系數(shù)的統(tǒng)計特性劃定區(qū)間范圍,以及相鄰2個深度層之間聲速差異是否滿足閾值,進行神經(jīng)網(wǎng)絡(luò)訓(xùn)練。

將XCTD的測量時間換算距離2023年3月的秒計時輸入到訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)中反演出XCTD的聲速剖面,與實際測量的XCTD聲速剖面求解其RMSE值作為精度指標。

由于神經(jīng)網(wǎng)絡(luò)擬合具有隨機性,這里使用擬合100次的EOF重構(gòu)系數(shù)求平均來作為最終反演的聲速剖面結(jié)果。

2 實驗結(jié)果與分析

2.1 數(shù)據(jù)來源

本文所使用的數(shù)據(jù)是2023年南海海試數(shù)據(jù),經(jīng)處理后包括7條全水深(3 500 m)CTD數(shù)據(jù)、7條全水深(3 500 m)SVP數(shù)據(jù)和9條部分水深(2 200 m)XCTD數(shù)據(jù),如圖3所示。用于匹配場算法的數(shù)據(jù)有14條全水深(3 500 m)CTD與SVP數(shù)據(jù)、母船走航時的時延數(shù)據(jù),以及在4個應(yīng)答器上方和其中心靜止時的時延數(shù)據(jù)。

圖3 預(yù)處理后聲速剖面數(shù)據(jù)Fig.3 Preprocessed sound velocity profile data

圖4為實測聲速剖面與其均值偏差剖面圖,可以看到,其在1 500 m以后趨于穩(wěn)定,并且相同儀器在深海與自身平均值的偏差趨向于0。圖3末尾的分叉是由于SVP是直接測量,而CTD是通過經(jīng)驗公式間接測量,測量機理不同導(dǎo)致的系統(tǒng)偏差和數(shù)據(jù)預(yù)處理對深度不夠的聲速剖面進行線性延拓導(dǎo)致。

2.2 實驗分析

聲速剖面采用的是W.D.Wilson簡化公式:

C=1 449.2+4.6T-0.055T2+0.000 29T3+

(1.34-0.01T)(S-35)+0.016D

經(jīng)過聲速剖面預(yù)處理插值為1 m為間隔的標準聲速剖面層。對于某些聲速剖面深度不夠的采用線性延拓。

使用匹配場處理需要參考聲速剖面樣本的特征向量,以及應(yīng)答器、換能器的位置和時延信息。圖5為母船走航時的位置信息構(gòu)成的軌跡圖。采用每10 min為一組進行聲速剖面反演,其代表時間為中間時刻。

圖5 匹配場(MFP)使用的時延數(shù)據(jù)軌跡Fig.5 Time-delayed data traces used by matched field (MFP)

圖6為通信數(shù)據(jù)預(yù)處理,由于通信數(shù)據(jù)的質(zhì)量影響匹配場反演,需要對通信數(shù)據(jù)進行質(zhì)量控制。采用母船換能器與海底應(yīng)答器的坐標距離除以1 500 m/s簡略計算時延,用計算的時延與真實測量的時延信息進行對比,觀測其變化趨勢。設(shè)定閾值從而進行數(shù)據(jù)的質(zhì)量控制。

圖6 時延數(shù)據(jù)預(yù)處理Fig.6 Preprocessing of time delayed data

本文使用平均聲速剖面(AVG)、僅使用實測數(shù)據(jù)的BP神經(jīng)網(wǎng)絡(luò)和參考通信數(shù)據(jù)的MFPBP。3種算法進行對比。AVG計算公式為:

式中:CAVG,j為平均聲速剖面(AVG)第j層的聲速值;Cref,j為參考聲速剖面第j層的聲速值,本文指代CTD與SVP的14條聲速剖面;n為參考聲速剖面的個數(shù)。

表1 EOF重構(gòu)前5個階次的方差貢獻率

表2 結(jié)果精度RMSETable 2 Results accuracy RMSE

圖7 AVG、BP、MFPBP反演聲速剖面與XCTD實測聲速剖面RMSEFig.7 AVG, BP, MFPBP inverted sound velocity profiles and XCTD measured sound velocity profile RMSE

圖8 BP反演EOF第一模態(tài)系數(shù)與XCTD實測聲速剖面RMSEFig.8 BP inversion EOF first mode coefficient and XCTD measured sound velocity profile RMSE

圖9 MFPBP反演EOF第一模態(tài)系數(shù)與XCTD實測聲速剖面RMSEFig.9 MFPBP inversion EOF first mode coefficient and XCTD measured sound velocity profile RMSE

從圖10擬合聲速剖面圖可以得到,MFPBP反演的聲速剖面與實測XCTD最為接近。BP反演的聲速剖面在前5條驗證數(shù)據(jù)在1 000～1 500 m處有較大突起。

圖10 AVG、BP、MFPBP擬合聲速剖面和XCTD實測聲速剖面示意Fig.10 Schematic diagram of AVG, BP, MFPBP fitted sound velocity profiles and XCTD measured sound velocity profile

圖11為通信數(shù)據(jù)與實測聲速剖面數(shù)據(jù)的分布,其x軸為走航時的時間,是距離2023年3月的秒計時,y軸為由通信數(shù)據(jù)經(jīng)過匹配場反演得到的聲速剖面?zhèn)€數(shù)累加。從表2可以看出AVG的精度在7、8、9的精度較高,這是由于參考樣本的區(qū)間范圍在其附近,其平均值更能代表這附近的聲速剖面。而使用了母船走航時的時延數(shù)據(jù)從圖11可以看出其分布比較均勻和連續(xù),其中出現(xiàn)時延數(shù)據(jù)間斷可能是因為天氣等其他原因避險所致。相較于AVG使用了匹配場反演的聲速剖面MFPBP方法前6個驗證數(shù)據(jù)精度提升幅度分別為0.693、0.569、0.789、0.560、0.348、0.477,這是由于前6個驗證數(shù)據(jù)范圍內(nèi)沒有實測參考數(shù)據(jù),AVG算法精度不高,而MFPBP算法能夠有效利用通信數(shù)據(jù)增加參考樣本。第7個驗證數(shù)據(jù)其精度相較于前6個提升有明顯下降,提升幅度分別為0.187。而第8條驗證數(shù)據(jù)和AVG精度相當。這可能是由于CTD的測量時段在其附近CTD平均值更能體現(xiàn)這片時間范圍內(nèi)的聲速剖面。第9條驗證數(shù)據(jù)從圖9單次反演可以看出第一模態(tài)波動范圍-20～50,波動幅度達到70,其精度RMSE也在0.5～1.9范圍波動,可以看出單次反演其精度和第一模態(tài)波動相較于有匹配場反演數(shù)據(jù)的前8條驗證數(shù)據(jù)波動巨大,穩(wěn)定性差。從整體來看傳統(tǒng)AVG的整體RMSE為1.079 m/s,而MFPBP反演的整體RMSE能夠達到0.665 m/s,提高了38.4%。可以看出對于使用通信信息利用MFP反演的聲速剖面對于參考樣本過少時有著一定的可用性。

圖11 通信數(shù)據(jù)與聲速剖面數(shù)據(jù)的時間信息Fig.11 Time information for communication data and sound velocity profile data

3 結(jié)論

1)聲速剖面反演結(jié)果表明,單純用神經(jīng)網(wǎng)絡(luò)反演RMSE為1.315 m/s,如果添加匹配場反演的結(jié)果RMSE能夠達到0.665 m/s,其精度提升了49.4%且穩(wěn)定性更好。相較于傳統(tǒng)AVG的RMSE為1.079 m/s,提高了38.4%。經(jīng)過MFP進行加密后的BP神經(jīng)網(wǎng)絡(luò)擬合出的效果更好且更穩(wěn)定。

2)通過對比前8條有匹配場反演數(shù)據(jù)的驗證數(shù)據(jù)和第9條沒有匹配場反演數(shù)據(jù)的驗證數(shù)據(jù)從單次反演EOF第一模態(tài)的波動情況,可以看出匹配場對BP神經(jīng)網(wǎng)絡(luò)的穩(wěn)定性有較大影響。

存在的問題:使用匹配場加密樣本量時間復(fù)雜度過高,需要對其進行優(yōu)化。使用MFPBP構(gòu)建聲速時間場精度的好壞還需從定位、導(dǎo)航方面進行分析。