EOF重構(gòu)聲速剖面對(duì)深水多波束的聲速改正分析

張孝首周興華*唐秋華王盼龍丁繼勝王永康

(1.山東科技大學(xué) 測(cè)繪科學(xué)與工程學(xué)院,山東 青島266590;2.自然資源部 第一海洋研究所,山東 青島266061)

海洋聲速是海洋環(huán)境觀測(cè)的基本要素之一,主要受溫度、鹽度與壓力的影響,準(zhǔn)確的聲速剖面對(duì)各種海洋聲學(xué)測(cè)量均具有重要的意義。大洋聲速剖面大部分符合混合層、主躍層與深海等溫層的三層剖面結(jié)構(gòu),具有一定的規(guī)律性。目前聲速剖面測(cè)量主要有直接測(cè)量和間接測(cè)量?jī)煞N方法。直接測(cè)量法通常采用聲速剖面儀直接測(cè)量聲速,測(cè)量精度較高;間接測(cè)量通常利用CTD 測(cè)量得到的溫度、鹽度和壓力等數(shù)據(jù),根據(jù)聲速經(jīng)驗(yàn)公式計(jì)算得到聲速,該方法的精度主要取決于所選擇的聲速經(jīng)驗(yàn)公式與儀器精度[1-2]。在大洋科考中,無(wú)論是采用直接法還是間接法,全水深聲速剖面都需停船定點(diǎn)作業(yè),對(duì)于4 000 m 的全水深聲速剖面測(cè)量,大約耗時(shí)5～6 h,采集效率較低,成本較高。在深海多波束海底地形測(cè)量過(guò)程中,通常采用定點(diǎn)CTD 作業(yè)和走航XCTD 相結(jié)合等方式獲取聲速數(shù)據(jù),但XCTD 作業(yè)深度有限,無(wú)法獲取全海深聲速剖面,易導(dǎo)致因聲速數(shù)據(jù)不準(zhǔn)確產(chǎn)生的多波束測(cè)量誤差,從而影響海底地形成果精度。為有效改正此種原因?qū)е碌恼`差,可以利用聲速剖面結(jié)構(gòu)參數(shù)化模型實(shí)現(xiàn)對(duì)聲速的有效表示,模型主要分為兩類,一類為解析函數(shù)模型,另一類為經(jīng)驗(yàn)正交函數(shù)模型[3](Empirical Orthogonal Functions,EOF)。

解析函數(shù)模型利用一系列數(shù)學(xué)表達(dá)式及參數(shù)表示聲速隨深度的變化,此類模型的缺點(diǎn)是公式較復(fù)雜,應(yīng)用者需要具備較強(qiáng)的數(shù)學(xué)基礎(chǔ),不利于推廣應(yīng)用。EOF 函數(shù)模型也叫主成分分析法,是描述聲速剖面最有效的基函數(shù),其將具有相同特征的聲速剖面群進(jìn)行模態(tài)分解,得到基函數(shù),一般3～6個(gè)基函數(shù)就可以近似表示任意一個(gè)聲速剖面[4-6]。根據(jù)沈遠(yuǎn)海[6]、張旭等[7]、孫文川等[8]的研究表明,EOF 方法構(gòu)建淺水聲速剖面場(chǎng)具有良好的效果,在構(gòu)建深遠(yuǎn)海聲速剖面場(chǎng)時(shí)的效果研究較少[11]。本文基于東南印度洋海洋調(diào)查工作,對(duì)獲取的CTD 實(shí)測(cè)數(shù)據(jù)計(jì)算得到的全水深聲速剖面進(jìn)行計(jì)算,利用EOF分析方法建立整個(gè)作業(yè)區(qū)的聲速剖面場(chǎng)模型;在此基礎(chǔ)上,提取深水多波束測(cè)量所需聲速剖面,對(duì)多波束水深地形數(shù)據(jù)進(jìn)行聲速改正,實(shí)驗(yàn)結(jié)果表明,通過(guò)EOF方法表示的聲速剖面可以有效對(duì)多波束測(cè)量數(shù)據(jù)進(jìn)行聲速改正。

1 EOF分析法構(gòu)建聲速剖面場(chǎng)原理

設(shè)測(cè)區(qū)有N個(gè)實(shí)測(cè)聲速剖面,每個(gè)聲速剖面有M個(gè)均勻深度的聲速值,(xi,yi)為聲速剖面對(duì)應(yīng)的地理坐標(biāo),xi為緯度,yi為經(jīng)度,將N個(gè)聲速剖面表示為矩陣形式為[12]：

式中：每一列代表一個(gè)聲速剖面,每一行代表各個(gè)聲速剖面同一深度層的聲速。

對(duì)式(1)中的每一行取平均值,得到所有聲速剖面同一深度層的平均聲速值,最終得到平均聲速矩陣,將聲速剖面矩陣與平均聲速矩陣相減得到聲速擾動(dòng)矩陣ΔCM×N,進(jìn)而得到擾動(dòng)矩陣的協(xié)方差矩陣COVM×N：

計(jì)算協(xié)方差矩陣的特征值λ(λ1,λ2,…,λN)和特征向量V(V1,V2,...,VN)：

λN對(duì)應(yīng)的特征向量VN即為EOF的一個(gè)模態(tài),V即為EOF空間函數(shù)。

將EOF投影到聲速擾動(dòng)矩陣,得到所有空間特征向量對(duì)應(yīng)的時(shí)間系數(shù)(主成分)PCM×N：

式中,PCM×N中的每一列為對(duì)應(yīng)聲速擾動(dòng)剖面的時(shí)間系數(shù)。

聲速剖面矩陣可以表示為以下形式：

每一個(gè)特征向量對(duì)應(yīng)的特征值表示此特征向量的權(quán)重,根據(jù)張志偉等的研究,前六階經(jīng)驗(yàn)正交函數(shù)就可以有效表示聲速剖面,因此,本文用前六階經(jīng)驗(yàn)正交函數(shù)重構(gòu)聲速剖面,表達(dá)式如下：

運(yùn)用雙線性插值方法對(duì)得到的時(shí)間系數(shù)進(jìn)行插值,可得到任意一點(diǎn)的時(shí)間系數(shù)：

根據(jù)最小二乘法原理X=(ATA)-1(ATA),得到待求參數(shù)矩陣a4×6。

2 深海聲速剖面場(chǎng)建立

2.1 CTD數(shù)據(jù)預(yù)處理

本研究選取東印度洋海洋調(diào)查工作中獲取的17個(gè)CTD 站位數(shù)據(jù),其中全水深3個(gè)(8、12及13號(hào)站位),非全水深14個(gè)(測(cè)站位置如圖1所示)。首先將各站位CTD 測(cè)量得到的溫度、鹽度數(shù)據(jù)按1 m 深度間隔導(dǎo)出,溫度與鹽度變化如圖2所示,然后再將14個(gè)非全水深CTD 的溫度、鹽度數(shù)據(jù)擬合計(jì)算至4 000 m。

圖1 CTD 測(cè)站位置Fig.1 Location of CTD stations

圖2 溫度與鹽度變化Fig.2 Changes in temperature and salinity

通過(guò)3個(gè)全水深CTD 數(shù)據(jù)分析后發(fā)現(xiàn),該區(qū)域海洋溫度與鹽度在2 000 m 以下變化較小(圖2),其中溫度我們可將其變化近似為線性變化,可通過(guò)有理擬合公式進(jìn)行擬合計(jì)算：

式中：x為深度(m);y為溫度(℃);p1,p2,q1和q2為參數(shù)。

為驗(yàn)證溫度擬合公式的精度,隨機(jī)選取一個(gè)全水深測(cè)站所測(cè)溫度數(shù)據(jù),對(duì)1 000～2 000 m 溫度數(shù)據(jù)進(jìn)行擬合并預(yù)測(cè)2 000～4 000 m 的溫度,預(yù)測(cè)溫度誤差見圖3,可知預(yù)測(cè)最大誤差為0.07℃,引起的聲速誤差為0.3 m/s,表明運(yùn)用式(8)可以對(duì)溫度進(jìn)行有效的擬合預(yù)測(cè)。研究區(qū)內(nèi)距離較近的CTD 站位鹽度基本相同,可采用鄰近的全水深鹽度數(shù)據(jù)對(duì)非全水深鹽度數(shù)據(jù)進(jìn)行替代。通過(guò)以上計(jì)算,最終得到17組深度為4 000 m、間隔為1 m 的溫度、鹽度數(shù)據(jù)。

圖3 溫度擬合誤差Fig.3 The fitting error of temperature

圖4 17組聲速剖面Fig.4 17 groups of sound velocity profiles

2.2 CTD聲速剖面計(jì)算

基于CTD 獲取的溫、鹽、深數(shù)據(jù),通過(guò)聲速剖面經(jīng)驗(yàn)公式可以計(jì)算得到各測(cè)站聲速剖面,目前國(guó)內(nèi)外比較認(rèn)可的聲速剖面經(jīng)驗(yàn)公式有Chen-Millero[13],Wilson[14],Del Grosso[13-14]及Leroy[15]等,各模型的溫度、鹽度、壓力適用范圍有所差別。根據(jù)黃辰虎等的 研 究,Del Grosso、Chen-Millero、C.C.Leroy、Coppens四個(gè)聲速模型在全球海域具有較高精度及適用性,Del Grosso及Chen-Millero兩個(gè)模型具有最優(yōu)精度。在以Del Grosso作為最優(yōu)模型的假設(shè)下,C.C.Leroy公式與Del Grosso公式計(jì)算得到的聲速差值約為±0.1 m/s,隨著水深的增加,互差趨于0.04 m/s,且使用C.C.Leroy 模型計(jì)算聲速無(wú)需測(cè)量壓力[15-17],為消除對(duì)壓力數(shù)據(jù)延伸引起的聲速計(jì)算誤差,本文選取C.C.Leroy 聲速剖面經(jīng)驗(yàn)公式計(jì)算聲速,C.C.Leroy算法的表達(dá)式[15]：

式中：C為聲速;T為溫度(℃);S為鹽度;Z為水深(m);φ為測(cè)量處緯度。根據(jù)式(9)計(jì)算,共得到17組水深為4 000 m,間隔為1 m 的聲速剖面(圖4)。

2.3 聲速剖面的EOF表示

研究區(qū)17組聲速剖面包含3組全水深聲速剖面,14組非全水深聲速剖面,本文利用2組全水深聲速剖面,14組非全水深聲速剖面構(gòu)建兩個(gè)聲速剖面場(chǎng),最大水深分別為1 500和4 000 m,站號(hào)13的全水深聲速剖面數(shù)據(jù)作為驗(yàn)證。

2.3.1 1 500 m 水深聲速剖面場(chǎng)建立

選取1～1 500 m 聲速剖面,對(duì)所選數(shù)據(jù)組成的聲速矩陣進(jìn)行經(jīng)驗(yàn)正交分解,將特征向量按方差貢獻(xiàn)率從大到小的順序進(jìn)行排列,張旭等[3,7]認(rèn)為前3～6組特征向量累積方差貢獻(xiàn)率為89.4%～96.6%。本文取前6組EOF特征向量構(gòu)建聲速剖面場(chǎng),前6組特征向量累積方差貢獻(xiàn)率為95.01%。用預(yù)留全水深聲速剖面進(jìn)行驗(yàn)證,對(duì)重構(gòu)聲速剖面與實(shí)測(cè)聲速剖面進(jìn)行對(duì)比,結(jié)果如圖5。

由圖5可知,重構(gòu)聲速剖面與實(shí)測(cè)聲速剖面非常接近,無(wú)明顯差異點(diǎn),0～300 m 聲速誤差逐漸變大,300～700 m 聲速誤差逐漸變小,700～1 500 m 聲速誤差無(wú)顯著變化。聲速誤差最大值約為4 m/s,聲速差值較大處出現(xiàn)在混合層,深度約為300 m,且在混合層差值跳躍較大,這與溫度在海水中的變化基本吻合,為驗(yàn)證,對(duì)第13組CTD 溫度數(shù)據(jù)求導(dǎo),結(jié)果如圖6。

由圖6可知,在700 m 以內(nèi)導(dǎo)數(shù)絕對(duì)值較大,最大值約為200 m,700～1 000 m 導(dǎo)數(shù)絕對(duì)值相對(duì)較小,1 000～2 000 m 導(dǎo)數(shù)絕對(duì)值接近0,且基本無(wú)變化,2 000 m 以后基本為0,證明溫度在700 m 以內(nèi)變化較快,200 m 左右變化最為劇烈,700～1 000 m 變化逐漸趨于平緩,1 000～2 000 m 勻速下降,2 000 m 以后基本沒(méi)有變化?；旌蠈佑捎谑芄庹沼绊?一天之內(nèi)同一水深溫度變化較大,同一時(shí)間垂直方向溫度也存在較大差異,在影響聲速的3個(gè)主要因素中,又以水溫變化對(duì)聲速的影響最大,在深度與鹽度不變的情況下,溫度變化1 ℃,聲速變化約為4.5 m/s。溫度變化越快,聲速變化也越快,重構(gòu)聲速剖面場(chǎng)所需的特征向量也就越多。本文構(gòu)建聲速剖面場(chǎng)所用聲速剖面的聲速按照1 m 深度間隔均勻分布,同一深度構(gòu)建聲速剖面場(chǎng)所用聲速剖面數(shù)相同,在聲速剖面數(shù)相同的情況下,溫度變化越復(fù)雜,構(gòu)建的聲速剖面場(chǎng)精度越差,實(shí)測(cè)聲速與重構(gòu)聲速剖面差值也越大。

圖5 1 500 m 聲速誤差Fig.5 Sound velocity errors within a water depth of 1 500 m

經(jīng)計(jì)算,聲速誤差的標(biāo)準(zhǔn)差為1.383 4 m/s,其中大于3 m/s的數(shù)量為63 個(gè),占總數(shù)比例為4.2%;2～3 m/s為118個(gè),占總數(shù)比例為7.9%;1～2 m/s為527個(gè),占總數(shù)比例為35.1%;聲速差小于1 m/s的為792個(gè),占總數(shù)比例為52.8%,證明利用EOF 方法構(gòu)建聲速剖面場(chǎng)效果較好。

2.3.2 4 000 m 水深聲速剖面場(chǎng)建立

選取1～4 000 m 聲速剖面,對(duì)照構(gòu)建1 500 m 聲速剖面場(chǎng)的方法建立4 000 m 聲速剖面場(chǎng),前6組特征向量累積方差貢獻(xiàn)率為94.34%,用預(yù)留全水深聲速剖面進(jìn)行驗(yàn)證,對(duì)重構(gòu)聲速剖面與實(shí)測(cè)聲速剖面進(jìn)行對(duì)比,結(jié)果如圖7。

聲速誤差的標(biāo)準(zhǔn)差為0.880 2 m/s,其中大于3 m/s的數(shù)量為64個(gè),占總數(shù)比例為1.6%;2～3 m/s為117個(gè),占總數(shù)比例為2.925%;1～2 m/s為643個(gè),占總數(shù)比例為16.075%;聲速差小于1 m/s的為3 176個(gè),占總數(shù)比例為79.4%。與構(gòu)建的1 500 m 聲速剖面場(chǎng)進(jìn)行對(duì)比可知,將聲速延伸至4 000 m 并不會(huì)降低聲速剖面場(chǎng)的精度。

圖7 4 000 m 聲速誤差Fig.7 Sound velocity errors within a water depth of 4 000 m

3 重構(gòu)聲速剖面對(duì)深水多波束水深數(shù)據(jù)進(jìn)行改正

為驗(yàn)證EOF方法構(gòu)建的聲速剖面場(chǎng)對(duì)深水多波束進(jìn)行聲速改正的效果,利用構(gòu)建的4 000 m 聲速剖面場(chǎng),選取同一區(qū)域EM122深水多波束測(cè)量的6條測(cè)線進(jìn)行聲速改正,測(cè)線位置如圖8所示。其中測(cè)線1和測(cè)線2水深約為3 000 m、測(cè)線3和測(cè)線4水深約為4 000 m、測(cè)線5和測(cè)線6水深約為5 000 m。

圖8 多波束測(cè)線位置Fig.8 Locations of multi-beam survey lines

測(cè)線1實(shí)測(cè)聲速剖面深度為1 500 m,測(cè)線2、測(cè)線3和測(cè)線4實(shí)測(cè)聲速剖面深度為2 000 m,測(cè)線5和測(cè)線6實(shí)測(cè)聲速剖面深度為4 500 m。為驗(yàn)證本方法的有效性,采用Caris10.0軟件進(jìn)行多波束數(shù)據(jù)處理,對(duì)選取的6條測(cè)線分別運(yùn)用實(shí)測(cè)聲速剖面和EOF方法重構(gòu)聲速剖面進(jìn)行聲速改正,將改正后水深點(diǎn)輸出并進(jìn)行比對(duì),2種聲速剖面改正比對(duì)結(jié)果如圖9～11,圖中顏色表示差值大小,圖中箭頭表示測(cè)量船的航向。對(duì)比對(duì)結(jié)果進(jìn)行最大值、最小值、平均值及中誤差計(jì)算,結(jié)果如表1。

圖9 3 000 m 水深差值Fig.9 Differences within a water depth range of 3 000 m

圖10 4 000 m 水深差值Fig.10 Differences within a water depth range of 4 000 m

圖11 5 000 m 水深差值Fig.11 Differences within a water depth range of 5 000 m

表1 兩種聲速剖面改正結(jié)果誤差統(tǒng)計(jì)(m)Table 1 Error statistics of the results corrected by using two types of sound velocity profiles(m)

由圖11和表1中可知水深約為5 000 m 的測(cè)線5和測(cè)線6中誤差最小,約為0.7 m,相對(duì)中誤差為0.014%,誤差絕對(duì)值的最大值約為5 m,相對(duì)誤差為0.1%,且誤差范圍為最小,約為2.5 m,測(cè)線5和測(cè)線6聲速改正所用重構(gòu)聲速剖面深度接近實(shí)測(cè)聲速剖面深度,可以認(rèn)為2種聲速剖面深度不同對(duì)聲速改正結(jié)果帶來(lái)的誤差可以忽略不計(jì),誤差的主要來(lái)源為同一深度層的聲速差。測(cè)線5和測(cè)線6對(duì)比的結(jié)果表明,通過(guò)EOF方法重構(gòu)聲速剖面對(duì)多波束數(shù)據(jù)進(jìn)行聲速改正的結(jié)果接近實(shí)測(cè)聲速剖面對(duì)多波束數(shù)據(jù)進(jìn)行聲速改正的結(jié)果。

從圖10和表1中可以看出測(cè)線3和測(cè)線4的中誤差最大,約為17 m,相對(duì)中誤差為0.425%,接近《海道測(cè)量規(guī)范》[18]規(guī)定的1%的中誤差的一半,誤差絕對(duì)值的最大值約為58 m,相對(duì)誤差為1.45%,超過(guò)了《海道測(cè)量規(guī)范》規(guī)定的1%的中誤差,其誤差范圍也為最大,約為85 m,這是因?yàn)闇y(cè)線3和測(cè)線4進(jìn)行聲速改正所用實(shí)測(cè)聲速剖面深度為2 000 m,遠(yuǎn)小于本文重構(gòu)聲速剖面的4 000 m 深度,用2 000 m 處聲速表示2 000 m 以下聲速會(huì)對(duì)水深點(diǎn)最終結(jié)果帶入較大誤差,應(yīng)用本文方法構(gòu)建的聲速剖面進(jìn)行聲速改正可以有效解決此類問(wèn)題,提高測(cè)量成果精度,對(duì)邊緣波束的改正效果最為明顯;從圖9和表1中可以看出水深約為3 000 m的測(cè)線1和測(cè)線2中誤差約為5 m,相對(duì)中誤差為0.17%,誤差絕對(duì)值的最大值約為20 m,相對(duì)誤差為0.67%,誤差范圍大小約為26 m,大于測(cè)線5和測(cè)線6,小于測(cè)線3和測(cè)線4,這是因?yàn)闇y(cè)線1和測(cè)線2的水深小于測(cè)線3和測(cè)線4的水深,用聲速剖面最大深度點(diǎn)處的聲速表示最深點(diǎn)以下聲速對(duì)水深點(diǎn)結(jié)果帶來(lái)的誤差小于測(cè)線3和測(cè)線4。

綜上所述,在5 000 m 范圍內(nèi),利用EOF方法構(gòu)建的聲速剖面與利用實(shí)測(cè)聲速剖面對(duì)深水多波束數(shù)據(jù)進(jìn)行聲速改正最終得到的水深相差很小,可以滿足深遠(yuǎn)海水深測(cè)量的要求。

4 結(jié) 論

本文對(duì)印度洋實(shí)測(cè)CTD 數(shù)據(jù)運(yùn)用C.C.Leroy聲速經(jīng)驗(yàn)公式計(jì)算得到各測(cè)站聲速剖面,應(yīng)用有理擬合公式對(duì)非全水深聲速剖面進(jìn)行延伸,在此基礎(chǔ)上利用EOF方法建立研究區(qū)聲速剖面場(chǎng),通過(guò)分析處理并與實(shí)測(cè)聲速剖面對(duì)比,建立的聲速剖面場(chǎng)與實(shí)測(cè)聲速剖面具有良好的符合性,得到結(jié)論：

1)利用EOF方法重構(gòu)聲速剖面,對(duì)深水多波束水深地形數(shù)據(jù)進(jìn)行聲速改正。在5 000 m 水深范圍內(nèi),應(yīng)用重構(gòu)聲速剖面與實(shí)測(cè)全水深聲速剖面進(jìn)行聲速改正的水深點(diǎn)最大誤差為5 m,滿足深遠(yuǎn)海水深測(cè)量的要求。

2)在無(wú)法獲取調(diào)查區(qū)CTD 站位的全水深聲速剖面的情況下,對(duì)非全水深聲速剖面應(yīng)用簡(jiǎn)單有效的有理擬合公式進(jìn)行延伸,不會(huì)影響構(gòu)建聲速剖面場(chǎng)的精度,并可有效改善邊緣波束的精度。

3)在深遠(yuǎn)海多波束地形測(cè)量無(wú)法獲取實(shí)時(shí)、全水深聲速剖面的情況下,EOF 反演聲速剖面方法可以作為一種實(shí)測(cè)聲速剖面的有效補(bǔ)充。