基于分層EOF 的深海聲速剖面時變特征建模

張林虎，劉 楊*，劉焱雄，陳冠旭，李夢昊

（1.自然資源部 第一海洋研究所,山東 青島 266061；2.自然資源部 海洋測繪重點實驗室,山東 青島 266061）

海洋占地球表面積71%，平均深度3 700 m[1]，如今正成為世界各國爭奪資源的主戰(zhàn)場。一切海洋活動都離不開精密的位置服務(wù)，而深海聲速隨海洋動態(tài)環(huán)境在1 400～1 600 m/s 的某個區(qū)間內(nèi)變化[2]，成為制約水聲精密定位精度的主要因素，因此需對深海聲速的變化規(guī)律展開研究。高精度的深海聲速剖面主要依靠聲速儀、溫鹽深儀進(jìn)行點位觀測獲取，受觀測手段制約，尚無法對整個目標(biāo)觀測海域的聲速場進(jìn)行全面、連續(xù)的觀測，因而需要根據(jù)實測數(shù)據(jù)構(gòu)建觀測區(qū)域內(nèi)的聲速剖面模型，以獲得特定時間和地點的聲速剖面。經(jīng)驗正交函數(shù)（Empirical Orthogonal Function,EOF）模型法是聲速剖面的主要構(gòu)建方法，它可以通過分析目標(biāo)海域內(nèi)歷史觀測數(shù)據(jù)的矩陣特征值與特征向量有效提取聲速剖面的主要特征[3]。

Leblanc 和Middleton[3]于1980 年首次提出了采用經(jīng)驗正交函數(shù)（EOF）分析海洋聲速剖面的方法，其后相關(guān)研究陸續(xù)開展。Davis 試驗驗證了EOF 作為基函數(shù)描述聲速剖面的方法具有最小均方根最優(yōu)性[4]；丁繼勝等[5]基于實測聲速剖面簇，利用經(jīng)驗正交函數(shù)計算了目標(biāo)海域任意位置處的聲速剖面，并依此對多波束邊緣數(shù)據(jù)進(jìn)行聲線折射改正，取得了較好的效果；張旭等[6]基于Argo 數(shù)據(jù)利用EOF 方法分析了臺灣以東海域聲速垂直結(jié)構(gòu)的時空變化特征；這些研究表明采用經(jīng)驗正交函數(shù)模型分析聲速剖面數(shù)據(jù)是可行的。不同階數(shù)的EOF 表達(dá)聲速剖面的精度和效率存在差異。Park 等研究認(rèn)為利用聲速剖面簇的前5 階EOF 就能較為準(zhǔn)確地描述任意聲速剖面[7]。周士弘和張茂有利用少于5 階的EOF 描述臺灣東北部海域聲速場，并嘗試采用EOF 方法預(yù)報聲速剖面[8]。國內(nèi)外相關(guān)工作集中于對淺海聲速剖面的研究，而對深海聲速剖面的研究則主要基于Argo 浮標(biāo)觀測數(shù)據(jù)開展。由于Argo 浮標(biāo)聲速剖面觀測數(shù)據(jù)分布稀疏（空間分辨率約為1°）且采樣時間間隔較大（10～30 d），不適用于局域較小時間尺度的深海聲速剖面特征分析。

由于不同深度處的深海聲速剖面存在不同的聲速變化規(guī)律和特征，以及聲速觀測設(shè)備的觀測深度范圍不同，所以需要采用分層的方式分析聲速剖面。劉楊范等[9]針對聲速在淺層水域的復(fù)雜性與深層水域的平穩(wěn)性特點，提出了雙層經(jīng)驗正交函數(shù)建模法，并采用自適應(yīng)方法確定出合理的階次。由于經(jīng)典深海聲速剖面分為混合層、深海聲道層和深海等溫層三層[10]，因而采用雙層正交函數(shù)建模法不適用于深海聲速剖面。聲速時空變化特征與海洋環(huán)境相關(guān)，陳小宇采用EOF 將聲速剖面分為背景場和擾動場，分析了聲速變化與潮汐現(xiàn)象之間的關(guān)系[11]；孫文舟等提出了一種簡化模型解析聲速剖面EOF 第一模態(tài)時間系數(shù)和空間函數(shù)變化規(guī)律，并認(rèn)為溫度是影響時空變化規(guī)律的主要因素[12]。

因此，本文基于較小時間尺度（時間間隔為6 h）的深海聲速剖面實測數(shù)據(jù)，首先分析深海聲速觀測數(shù)據(jù)的統(tǒng)計特征，結(jié)合經(jīng)典深海聲速剖面分層，將目標(biāo)海域聲速剖面分為上、中、下三層；然后，采用分層EOF 方法分析目標(biāo)海域不同深度層第一模態(tài)時間系數(shù)的變化規(guī)律，采用2 種時域擬合模型描述時變規(guī)律并比較其擬合精度，并利用等效平均聲速驗證擬合模型的合理性；最后，探究聲速剖面周期性變化特性的影響因素。

1 聲速剖面分層EOF 建模

1.1 深海聲速剖面分層深度選取

不同深度層聲速剖面隨時間和空間變化的特征不同，需要根據(jù)聲速時空變化特征對剖面進(jìn)行分層，以減小不同分層中聲速剖面時空變化特征對構(gòu)建模型的影響。本文利用不同聲速剖面簇的同深度標(biāo)準(zhǔn)差、聲速標(biāo)準(zhǔn)差梯度、平均聲速剖面垂直梯度三種聲速剖面統(tǒng)計特征對深海聲速剖面進(jìn)行分層。其中，不同聲速剖面簇的同深度標(biāo)準(zhǔn)差表示聲速隨時間變化的劇烈程度，聲速標(biāo)準(zhǔn)差梯度表示同深度標(biāo)準(zhǔn)差隨深度變化的程度，平均聲速剖面垂直梯度表示聲速剖面隨深度變化的程度。

假設(shè)觀測區(qū)域內(nèi)獲得M個N層的標(biāo)準(zhǔn)聲速剖面數(shù)據(jù)集。第m層的同深度聲速標(biāo)準(zhǔn)差 σm定 義為：

式中：M為第m層聲速觀測個數(shù)；Si為第m層第i次聲速觀測值；Sˉ為第m層聲速觀測均值。

第m層 的聲速標(biāo)準(zhǔn)差梯度?σm定義為：

式中：σm+1為第m+1層同深度聲速標(biāo)準(zhǔn)差；Δzm為第m層 的水層厚度。

第m層的聲速垂直梯度gm定義為：

式中：Sm+1為第m+1層聲速；Sm為第m層 聲速。

1.2 分層EOF 方法

假設(shè)觀測區(qū)域內(nèi)的聲速剖面有M個，經(jīng)過數(shù)據(jù)預(yù)處理，最終獲得M個N層的標(biāo)準(zhǔn)聲速剖面數(shù)據(jù)集S。根據(jù)經(jīng)典深海聲速剖面變化特點并結(jié)合目標(biāo)海域的聲速變化統(tǒng)計特征，可以將標(biāo)準(zhǔn)聲速剖面數(shù)據(jù)集S劃分為上、中、下三個深度區(qū)間層：

進(jìn)而，對協(xié)方差矩陣求解特征值和特征向量：

式中：λu、λb、λd分別為Ru、Rb、Rd的特征值矩陣，F(xiàn)u、Fb、Fd分別為Ru、Rb、Rd與λu、λb、λd對應(yīng)的特征向量。

將特征值按照從大到小的順序重新排列:

式中：λu、λb、λd分別為Ru、Rb、bd的排序后的特征值矩陣，F(xiàn)u、Fb、Fd別為Ru、Rb、Rd與λu、λb、λd對應(yīng)的特征向量，D1、D2-D1、N-D2分別為上、中、下深度區(qū)間層解算的特征值個數(shù)。

每一個深度層內(nèi)，模態(tài)系數(shù)矩陣A定義[14]為：

式中：aN為第N階模態(tài)系數(shù)。

基于EOF 模態(tài)矩陣的方差貢獻(xiàn)率[15]，分別在上、中、下三個深度區(qū)間層內(nèi)取前i、j、k個特征值，及其對應(yīng)的特征向量和模態(tài)系數(shù)矩陣。以下層為例，特征值矩陣 λdk、特征向量矩陣Fdk及模態(tài)系數(shù)矩陣Adk分別表示為：

同理可得到λui、Fui、Aui、λbj、Fbj、Abj。對模態(tài)系數(shù)矩陣進(jìn)行時變分析，可得到聲速隨時間的變化規(guī)律。對模態(tài)系數(shù)進(jìn)行擬合便可獲得測量期間任意時刻的模態(tài)系數(shù)并解算得到任意時刻的聲速剖面：

2 結(jié)果與分析

2.1 深海聲速剖面數(shù)據(jù)獲取及分層

實驗數(shù)據(jù)來自2021 年7 月在南海海域利用聲速儀（Sound Velocity Profile,SVP）及溫鹽深儀（Conductivity,Temperature,Depth,CTD）觀測的深海聲速剖面數(shù)據(jù)，其觀測間隔為6 h，時間跨度為7 d。采用聲速儀及溫鹽深儀觀測獲得數(shù)據(jù)后，需要預(yù)處理以獲得標(biāo)準(zhǔn)分層的聲速剖面數(shù)據(jù)，處理流程見圖1。其中，溫鹽深儀測量海水中的溫度、深度、鹽度信息后，根據(jù)聲速經(jīng)驗公式計算海水聲速值。陳長安等[16]對比了10 種聲速經(jīng)驗公式在南海海域的精度，建議使用Del Grosso 聲速經(jīng)驗公式計算聲速，本文采用Del Grosso 聲速經(jīng)驗公式計算聲速值[17-18]。聲速觀測過程中，受投放儀器設(shè)備故障及海流影響，觀測剖面深度不是順序增加，需基于深度對數(shù)據(jù)進(jìn)行重排序，并剔除粗差；對于偏離聲速剖面主廓線過大的聲速數(shù)據(jù)予以剔除。基于不同分層的數(shù)據(jù)分布情況進(jìn)行不同處理：觀測深度以內(nèi)的缺失數(shù)據(jù)進(jìn)行插值；觀測深度以外數(shù)據(jù)進(jìn)行數(shù)據(jù)延拓；對相同深度的觀測值采用調(diào)和平均方法獲得該深度處聲速值，對于標(biāo)準(zhǔn)間隔深度處的聲速則取前后一定深度區(qū)間內(nèi)聲速的加權(quán)平均值；最后對分層內(nèi)的聲速數(shù)據(jù)取平均以獲得深度間隔為0.5 m 的標(biāo)準(zhǔn)化聲速剖面。

圖1 聲速數(shù)據(jù)預(yù)處理流程Fig.1 Preprocessing flow of the sound velocity data

經(jīng)過數(shù)據(jù)預(yù)處理，SVP 測量獲得12 條有效數(shù)據(jù)，CTD獲得30 條有效數(shù)據(jù)（圖2）。聲速剖面具有典型的深海聲速剖面分層特征[10]，隨深度變化可劃分為表面層（混合層）、躍變層和深海等溫層，其中躍變層又分為季節(jié)性躍變層和主躍變層。聲速剖面在水深0～2 m 處，受儀器自身溫度的影響，聲速變化劇烈且差異較大；在水深2～30 m 左右處，考慮到風(fēng)浪的攪拌作用導(dǎo)致該水深區(qū)間水溫相近，聲速變化幅度較小，具有明顯的表面層特征；在水深50～1 100 m 處，聲速值隨深度增加而減小，為深海垂直結(jié)構(gòu)中的躍變層；在水深1 000 m 以深的聲速剖面結(jié)構(gòu)穩(wěn)定，層內(nèi)水溫與鹽度變化趨于一致，聲速主要受壓力變化影響，所以，聲速隨深度增加而增大，具有明顯的深海等溫層特征。

圖2 測區(qū)聲速剖面示意圖Fig.2 Sound velocity profiles in the study area

基于實測的聲速剖面數(shù)據(jù)，本文分別統(tǒng)計了SVP 和CTD 測得的聲速剖面簇的同深度標(biāo)準(zhǔn)差、同深度聲速標(biāo)準(zhǔn)差梯度以及平均聲速剖面的聲速垂直梯度，見圖3 和圖4。

由圖3 可以看出，不同深度聲速剖面標(biāo)準(zhǔn)差為0.02～3.58 m/s，其最大值出現(xiàn)在109 m 水深處，最小值出現(xiàn)在3 130 m 水深處。由圖4 可以看出，不同深度聲速剖面標(biāo)準(zhǔn)差為0.02～2.86 m/s，其最大值出現(xiàn)在87 m 水深處，最小值出現(xiàn)在2 650 m 水深處。綜合圖3 和圖4 可以發(fā)現(xiàn)：在水深0～500 m 左右，聲速標(biāo)準(zhǔn)差梯度數(shù)值和聲速剖面標(biāo)準(zhǔn)差隨深度變化劇烈且數(shù)值較大，聲速垂直梯度為負(fù)值且梯度絕對值由大到小，并趨近于0；在水深500～1 500 m 左右，同深度聲速標(biāo)準(zhǔn)差梯度數(shù)值變化趨緩且呈現(xiàn)變小趨勢，不同深度聲速剖面標(biāo)準(zhǔn)差變化劇烈程度變小且具有隨深度變化緩慢降低的趨勢，聲速垂直梯度緩慢趨近于0；在水深1 500～3 200 m 左右，不同深度聲速剖面標(biāo)準(zhǔn)差、同深度聲速標(biāo)準(zhǔn)差梯度和聲速垂直梯度基本穩(wěn)定，數(shù)值均趨近于0。

圖3 SVP 測量聲速剖面數(shù)據(jù)統(tǒng)計特征Fig.3 Statistical characteristics of the sound velocity profile data measured with SVP

圖4 CTD 測量數(shù)據(jù)聲速剖面統(tǒng)計特征Fig.4 Statistical characteristics of the sound velocity profile data measured with CTD

針對不同深度區(qū)間聲速變化差異并顧及分層點，應(yīng)盡量選擇在聲速變化較緩慢的深度[19]，故本文將實測數(shù)據(jù)聲速剖面分為0～550 m、450～1 550 m、1 450～3 200 m 三層，為避免出現(xiàn)邊界聲速不匹配的情況，不同分層之間設(shè)置了100 m 深度區(qū)域的重疊部分。

2.2 分層EOF 第一模態(tài)系數(shù)時變分析

對聲速剖面簇進(jìn)行分層EOF 分析時，在不同分層，采用適當(dāng)?shù)腅OF 階數(shù)能提高運算速度，降低計算數(shù)據(jù)量，因此，需要選取合適階數(shù)的特征向量和特征值。通常，EOF 各階次方差貢獻(xiàn)率及累計方差貢獻(xiàn)率[15]能夠反映模態(tài)分析的有效性，并為分層EOF 階數(shù)的選擇提供依據(jù)。前m階EOF 累積貢獻(xiàn)率Ekm定義為：

式中：k為分層內(nèi)協(xié)方差矩陣的行數(shù)；λi，λj分別為第i和第j個特征值。

對測量數(shù)據(jù)進(jìn)行分層EOF 分解，SVP、CTD、SVP+CTD 觀測數(shù)據(jù)的結(jié)果獲得上、中、下三層各階次方差貢獻(xiàn)率與階數(shù)的關(guān)系，上層聲速剖面利用前3～4 階獲得的反演聲速剖面，EOF 方差貢獻(xiàn)率即可達(dá)到99%；中、下層聲速剖面利用前2 階獲得的反演聲速剖面，EOF 方差貢獻(xiàn)率即可達(dá)到99%，如圖5 所示。整體而言，不同聲速剖面簇的第一階次分層EOF 方差貢獻(xiàn)率均在79%以上，分層EOF 第一模態(tài)系數(shù)能夠反映聲速剖面隨時間變化的主要特征。

圖5 EOF 累計方差貢獻(xiàn)率與階數(shù)的關(guān)系Fig.5 The relation between the accumulated variance contribution rate and the order of the EOF

本文通過分析分層EOF 第一模態(tài)系數(shù)研究聲速剖面的時域變化規(guī)律，假設(shè)聲速剖面時變規(guī)律CT符合正弦變化，則：

式中：TD為聲速剖面觀測的以天為周期的時間；A1、A21、A22為常數(shù)項和周期項系數(shù)，其中A21=A2·cos(φ)，A22=A2·sin(φ)。

基于SVP 和CTD 聲速剖面觀測數(shù)據(jù)集，以及2 種數(shù)據(jù)的所有觀測結(jié)果組成的聲速剖面觀測數(shù)據(jù)集（簡稱SVP+CTD 數(shù)據(jù)集）對深海聲速剖面進(jìn)行分層EOF 分析，采用一階傅里葉函數(shù)擬合分層EOF 第一模態(tài)系數(shù)的日周期變化規(guī)律，擬合結(jié)果（圖6）表明，上、中層EOF 第一模態(tài)系數(shù)的日變化規(guī)律均較明顯，其中中層較上層日變化規(guī)律更加明顯，其數(shù)據(jù)擬合誤差見表1。

表1 第一模態(tài)系數(shù)擬合殘差標(biāo)準(zhǔn)差Table 1 The standard deviation of fitted residuals of the first mode coefficients of EOF

圖6 不考慮長期線性變化項的第一模態(tài)系數(shù)時間變化Fig.6 The temporal variations of the first mode coefficient of EOF without considering the long-term linear variation term

聲速剖面的時域變化具有日周期、季節(jié)性和年際變化特征，聲速剖面測量時間跨度約為一周，其時域變化會受到季節(jié)性和年際長期項的影響，同時長期項引起的聲速剖面時域變化在觀測期間近似可視為線性變化，故在一階傅里葉擬合基礎(chǔ)上添加線性項并進(jìn)行對比分析。擬合日周期的分層EOF 第一模態(tài)系數(shù)變化的方法為：

式中：MN1(TD,H)為第一模態(tài)系數(shù)的時變規(guī)律項；TD為聲速剖面觀測的以天為周期的時間；H為聲速剖面觀測與首次聲速剖面觀測的時間延遲；a1、a2、a3、a4為擬合系數(shù)，由最小二乘法計算獲得。

對深海聲速剖面進(jìn)行分層EOF 分析，采用添加線性項的一階傅里葉函數(shù)擬合分層EOF 第一模態(tài)系數(shù)的日周期變化特征，擬合結(jié)果如圖7 所示。由圖7 可知：不同分層的EOF 第一模態(tài)系數(shù)時變規(guī)律特征不同；中層比上層第一模態(tài)系數(shù)日周期變化特征更明顯。對比2 種模型的擬合殘差標(biāo)準(zhǔn)差（表1）發(fā)現(xiàn)，所有分層及不同種類數(shù)據(jù)在添加線性項之后，擬合效果基本有所改進(jìn)，且中層改進(jìn)幅度大于其他層。

圖7 考慮長期線性變化項的第一模態(tài)系數(shù)時間變化Fig.7 The temporal variations of the first mode coefficient of EOF with the consideration of the long-term linear variation term

2.3 聲速剖面時變特征分析

在EOF 分析過程中，對分層EOF 模態(tài)系數(shù)進(jìn)行擬合可獲得測量期間任意時刻的模態(tài)系數(shù)，并解算得到任意時刻的聲速剖面，但模態(tài)系數(shù)不具有明確的物理意義。因此，本文采用添加線性項的一階傅里葉函數(shù)擬合分層等效平均聲速的日周期變化特征，基于等效平均聲速的日周期變化特征進(jìn)一步反映聲速剖面的時變特征。擬合日周期的分層等效平均聲速變化的方法為：

式中：ES(TD,H)為等效平均聲速的時間變化項；TD為聲速剖面觀測的以天為周期的時間；H為聲速剖面觀測與首次聲速剖面觀測的時間延遲；a1、a2、a3、a4為擬合系數(shù)，采用最小二乘法計算獲得。通過采用添加線性項的一階傅里葉函數(shù)擬合分層等效平均聲速結(jié)果發(fā)現(xiàn)，不同分層聲速剖面均具有日周期變化特性，中層日周期變化特征較上層日周期變化特征更明顯（表2 和圖8）。聲速剖面的等效平均聲速與第一模態(tài)系數(shù)變化規(guī)律一致，說明第一模態(tài)系數(shù)的時域變化特征反映聲速剖面的時域變化特征。

圖8 考慮長期線性變化項的等效平均聲速時間變化Fig.8 The temporal variations of the equivalent average sound velocity with the considerationof the long-term linear variation term

表2 等效平均聲速擬合殘差標(biāo)準(zhǔn)差（m·s-1）Table 2 The standard deviation (m·s-1) of fitted residuals of the equivalent average sound velocity

2.4 影響因素分析

EOF 第一模態(tài)系數(shù)變化的影響因素有溫度日變化和年變化引起的海水聲速周期性變化，且影響隨深度增大而減小[11]。本文首先將第一模態(tài)系數(shù)與溫度做相關(guān)性分析，結(jié)果見表3。從表3 中可以看出，分層EOF第一模態(tài)系數(shù)與溫度數(shù)據(jù)具有顯著相關(guān)性?？紤]到太陽輻射對海水溫度的影響主要集中于表層至500 m 以淺的區(qū)間[11]，推測中層聲速或溫度變化主要由潮流所致。為此，本文分析聲速與潮汐之間的關(guān)系。

表3 EOF 第一模態(tài)系數(shù)與等效平均溫度相關(guān)系數(shù)Table 3 The correlation coefficients between the first mode coefficient of EOF and equivalent average temperature

基于實驗海區(qū)附近東沙島驗潮站2021 年4 月1 日至5 月30 日的潮汐逐時預(yù)報數(shù)據(jù)（http://global-tide.nmdis.org.cn/），利用T_TIDE 潮汐處理軟件包[20]對潮汐數(shù)據(jù)進(jìn)行調(diào)和分析，獲得主要分潮的調(diào)和常數(shù)，結(jié)果見表4。由表4 可以看出，聲速實驗區(qū)域內(nèi)聲速剖面觀測時間段，潮汐的主要影響分潮為 O1分潮和 K1分潮。

表4 主要分潮調(diào)和常數(shù)Table 4 Tidal constants of the major tidal constituents

根據(jù)獲得的分潮振幅，即可獲得區(qū)域的潮汐類型[22]，潮汐類型指標(biāo)F計算方法為：

式中:HO1、HK1、HM2為O1、K1、M2分潮振幅。當(dāng)F≤0.5時，潮型為半日潮；當(dāng)0.5 ＜F≤2.0時，潮型為不規(guī)則半日潮；當(dāng)2.0 ＜F≤4.0時，潮型為不規(guī)則全日潮；當(dāng)F＞4.0時，潮型為全日潮[23]。

在聲速剖面觀測時間段，實驗海域潮汐的F值為2.47，潮汐類型為不正規(guī)全日潮，與余慕耕[24]對南海該區(qū)域潮汐分布特征的描述一致，與本文提取的聲速剖面時變特征也基本吻合。

3 結(jié)論

針對南海海域?qū)崪y的全海深聲速剖面，本文提出基于分三層的經(jīng)驗正交函數(shù)法分析聲速剖面的時域變化特征，并構(gòu)建局域小時間尺度聲速剖面模型對聲速時域變化特征進(jìn)行擬合。根據(jù)典型深海聲速剖面結(jié)構(gòu)和實測數(shù)據(jù)不同分層的統(tǒng)計特征，將實測數(shù)據(jù)深度分為上（0～550 m）、中（450～1 550 m）、下（1 450～3 200 m）三層。對每層聲速剖面進(jìn)行經(jīng)驗正交函數(shù)分解，獲得模態(tài)系數(shù)矩陣，并分析第一模態(tài)系數(shù)的時域變化，獲得其日變化特征，進(jìn)而采用2 種模型對聲速剖面第一模態(tài)系數(shù)進(jìn)行擬合，得到結(jié)論如下。

1）實驗海區(qū)聲速剖面分層EOF 第一模態(tài)系數(shù)及等效平均聲速具有日周期變化特征，上層聲速日周期變化特征不明顯，中層聲速日周期變化特征較明顯，下層聲速變化較小但仍具有日周期變化特征。

2）對于局部海域小時間尺度聲速剖面變化擬合過程，應(yīng)考慮聲速剖面長周期變化項的影響。

3）實驗海區(qū)聲速剖面EOF 第一模態(tài)系數(shù)變化與溫度顯著相關(guān)，聲速剖面時變特征與海區(qū)潮汐周期特征相似。

聲速剖面的時域變化特征影響因素的分析局限于對相關(guān)關(guān)系的研究。下一步本研究工作將對影響因素進(jìn)行全面分析并分析其因果關(guān)系，同時采用不同聲速觀測值進(jìn)行模型計算，由于觀測儀器精度和系統(tǒng)誤差存在差異，需要在聲速剖面數(shù)據(jù)集解算過程中對不同聲速剖面數(shù)據(jù)進(jìn)行定權(quán)處理。