基于最小二乘配置的陌生海域海底地形反演方法

范 雕,李?yuàn)檴?馮進(jìn)凱,黃 炎,范昊鵬,張金輝,李新星

1. 信息工程大學(xué),河南 鄭州 450001; 2. 軍事科學(xué)院國(guó)防科技創(chuàng)新研究院,北京 100071

海洋是人類(lèi)生存和發(fā)展的重要資源寶庫(kù),其面積占地球表面積71%左右[1]。文獻(xiàn)[2]研究表明目前依據(jù)船載聲學(xué)測(cè)量手段僅有15%左右的海域面積完成了空間分辨率約800 m的海底地形繪制,人類(lèi)對(duì)于海洋的認(rèn)識(shí)仍處于較低水平和層次[3]。文獻(xiàn)[4]認(rèn)為依托船舶為載體實(shí)施海底地形測(cè)量,若以空間網(wǎng)格大小為1′作為考慮對(duì)象,全球大約還有90%的海域?qū)儆谖幢豢碧綘顟B(tài),而且船舶測(cè)線(xiàn)間存在巨大的空間間隔,測(cè)量結(jié)果分布很不均勻,人類(lèi)對(duì)海底地形的描繪詳細(xì)程度遠(yuǎn)不及對(duì)月球等其他行星的了解深入[5]。

然而,海底地形地貌數(shù)據(jù)作為開(kāi)發(fā)和利用海洋的基礎(chǔ)性地理空間信息數(shù)據(jù),廣泛影響其他科學(xué)研究和應(yīng)用實(shí)踐。一方面,海底地形變化直接影響海洋環(huán)流和混合,從而間接影響地球氣候、生物多樣性和糧食產(chǎn)量。另一方面,海底地形和海底構(gòu)造的精細(xì)化表達(dá)有利于研究板塊構(gòu)造規(guī)律,有利于揭示海盆和海底地質(zhì)運(yùn)動(dòng)歷史,有利于預(yù)判發(fā)生地震、火山爆發(fā)和海嘯等區(qū)域位置。因此,過(guò)去100多年來(lái)國(guó)內(nèi)外學(xué)者持續(xù)不斷努力研究繪制全球范圍海底地形圖[6-7],囿于技術(shù)裝備發(fā)展水平制約,前期研究進(jìn)展較為緩慢。直到19世紀(jì)90年代中期,由于GEOSAT和ERS-1等執(zhí)行大地測(cè)量任務(wù)的測(cè)高衛(wèi)星數(shù)據(jù)不斷積累和應(yīng)用,全球海洋區(qū)域重力場(chǎng)精度和分辨率不斷提升,依托海洋重力信息反演海底地形研究不斷推進(jìn)[8-10]。至此,運(yùn)用空天探測(cè)技術(shù)獲取全球范圍的海底地形信息方才取得重要突破。文獻(xiàn)[11]利用衛(wèi)星測(cè)高重力異常數(shù)據(jù),依據(jù)回歸方法(又稱(chēng)S&S方法)建立了范圍涵蓋南北緯72°之間,分辨率約15 km的海深格網(wǎng)模型。隨后數(shù)十年時(shí)間,由于模型構(gòu)建數(shù)據(jù)量不斷豐富,數(shù)據(jù)精度不斷提升,學(xué)者們?cè)诟呔取⒏叻直媛始案鼜V覆蓋范圍的海底地形建模方面進(jìn)行了積極探索,比如:①豐富的船測(cè)海深數(shù)據(jù)可用于校準(zhǔn)和修正預(yù)測(cè)海深結(jié)果[12],以使構(gòu)建的海深模型網(wǎng)格密集程度更高;②波形重跟蹤技術(shù)的廣泛使用,在算法層面不斷提高測(cè)高結(jié)果精度,從而不斷提升海洋重力場(chǎng)質(zhì)量[13]。目前利用重力數(shù)據(jù)反演海底地形算法大致可分為統(tǒng)計(jì)法和解析法兩大類(lèi):統(tǒng)計(jì)法包括回歸分析方法[14]、模擬退火方法[15]、最小二乘配置[16]及神經(jīng)網(wǎng)絡(luò)[17]等;解析方法包括重力地質(zhì)方法[18-19]、導(dǎo)納函數(shù)方法[20]及顧及地形高次項(xiàng)影響的解析反演方法[21-22]等。然而,重力地質(zhì)方法、回歸分析方法、模擬退火方法在試驗(yàn)海區(qū)需要外部船測(cè)水深數(shù)據(jù)支撐,在無(wú)船測(cè)水深數(shù)據(jù)支撐的陌生海區(qū)該方法將不再適用,進(jìn)而極大削弱算法應(yīng)用范圍?；赑arker公式的解析反演方法不需要外部船測(cè)水深數(shù)據(jù)為補(bǔ)償,文獻(xiàn)[23]利用此方法開(kāi)展了無(wú)船測(cè)水深數(shù)據(jù)(利用先驗(yàn)海深模型代替)的海底地形建模,但是獲取該方法涉及的諸多地球物理參數(shù)并非易事,若全球海域均采用參數(shù)經(jīng)驗(yàn)值代替,毋庸置疑會(huì)造成最終反演精度損失,進(jìn)而導(dǎo)致反演結(jié)果可信度降低。

基于此,為擴(kuò)展無(wú)船測(cè)水深數(shù)據(jù)支持環(huán)境下陌生海域海底地形恢復(fù)方法,本文提出以衛(wèi)星測(cè)高重力數(shù)據(jù)和先驗(yàn)海底地形模型為輸入,運(yùn)用最小二乘配置方法(統(tǒng)計(jì)法),解決無(wú)船測(cè)數(shù)據(jù)支撐環(huán)境下高效構(gòu)建陌生海區(qū)海底地形模型難題的新思路?；静呗允?首先,在船測(cè)水深數(shù)據(jù)和重力數(shù)據(jù)已知海區(qū)統(tǒng)計(jì)解算出反映海區(qū)地形的相關(guān)協(xié)方差模型;然后,采用最小二乘配置方法,以海區(qū)重力數(shù)據(jù)為輸入,構(gòu)建相應(yīng)的海底地形模型并檢驗(yàn)算法適用性;最后,將已知海區(qū)數(shù)據(jù)統(tǒng)計(jì)結(jié)果(不同類(lèi)型參量的協(xié)方差函數(shù))遷移、推廣并應(yīng)用于與已知海區(qū)相似的船測(cè)數(shù)據(jù)陌生海區(qū),實(shí)現(xiàn)僅依靠重力數(shù)據(jù)的海底地形模型構(gòu)建。同時(shí)將本文反演結(jié)果與導(dǎo)納函數(shù)反演海底地形結(jié)果、DTU18、ETOPO1模型進(jìn)行比對(duì)分析,簡(jiǎn)要評(píng)估模型精度和效能。

1 原理與方法

1.1 最小二乘配置方法

利用解析方法(例如導(dǎo)納函數(shù))反演海底地形時(shí),算法模型涉及的地球參數(shù)求解并非易事。文獻(xiàn)[24]認(rèn)為海底地形和重力信息為兩類(lèi)隨機(jī)變量,從而依據(jù)重力異常數(shù)據(jù),采用統(tǒng)計(jì)方法恢復(fù)海底地形模型為

(1)

式中,b為海底地形;g為海面重力異常信息;μb和μg表示海底地形和海面重力異常的數(shù)學(xué)期望(可用相應(yīng)的平均值代替);Cbg為海底地形和重力異常的互協(xié)方差矩陣;Cgg為海面重力異常的自協(xié)方差。式(1)化簡(jiǎn)為

(2)

由式(2)可知,使用海面重力異常數(shù)據(jù),運(yùn)用最小二乘配置方法恢復(fù)海底地形信息關(guān)鍵在于獲取地形-重力異?；f(xié)方差函數(shù)和重力異常信息間自協(xié)方差函數(shù)。其中重力信息的自協(xié)方差函數(shù)求解在物理大地測(cè)量重力擬合推估研究較多[26],目前已有成熟的解決手段。然而海底地形和重力異常分屬于兩種不同類(lèi)型的參數(shù),因此如何精確求解地形-重力異?；f(xié)方差函數(shù)成了利用式(2)恢復(fù)海底地形的主要難題。

1.2 地形-重力協(xié)方差估計(jì)

將地形和重力異常作為一維隨機(jī)變量參考時(shí),地形-重力異?；f(xié)方差的無(wú)偏估計(jì)通過(guò)遍歷數(shù)據(jù)得到[27]

(3)

(4)

式(3)離散化的結(jié)果為

(5)

式中,l表示重力異常和地形間距離,通常取0,1,…,N-1。易知,式(5)重力-地形協(xié)方差估計(jì)是無(wú)偏估計(jì)

(6)

式中,E(·)表示期望算子。同理,二維重力-地形互協(xié)方差函數(shù)無(wú)偏估計(jì)表示為

(7)

借鑒重力異常擬合推估自協(xié)方差函數(shù)求解策略[26, 28],對(duì)于方形格網(wǎng)研究區(qū)域(N行和M列)重力-地形的互協(xié)方差可采用式(8)統(tǒng)計(jì)縱向和橫向互協(xié)方差

(8)

式中,q和qj分別表示橫向距離和縱向距離。最終將所有具有距離l的縱向協(xié)方差和橫向協(xié)方差取平均值作為協(xié)方差Cbg(l)的統(tǒng)計(jì)數(shù)值

(9)

協(xié)方差統(tǒng)計(jì)數(shù)值結(jié)果反映了海區(qū)相關(guān)數(shù)據(jù)變化的內(nèi)在規(guī)律,通常情況下需構(gòu)建強(qiáng)相關(guān)函數(shù)模擬海區(qū)數(shù)據(jù)協(xié)方差統(tǒng)計(jì)結(jié)果,進(jìn)而依托強(qiáng)相關(guān)函數(shù)擬合結(jié)果構(gòu)造相應(yīng)的(互)協(xié)方差矩陣,從而據(jù)最小二乘配置算法模型恢復(fù)海底地形信息。強(qiáng)相關(guān)函數(shù)如Inverse協(xié)方差函數(shù)、Gauss協(xié)方差函數(shù)、Hirvonen協(xié)方差函數(shù)、似高斯協(xié)方差函數(shù)、周期指數(shù)函數(shù)、Tscherning/Rapp協(xié)方差函數(shù)和Mortiz協(xié)方差函數(shù)等[29]。本文選擇較為簡(jiǎn)單的似高斯協(xié)方差函數(shù)作為海區(qū)協(xié)方差統(tǒng)計(jì)備選擬合函數(shù),表達(dá)式為

Cbg(l)=C0·e-al

(10)

式中,C0是兩點(diǎn)距離為零對(duì)應(yīng)的協(xié)方差函數(shù)值,表示協(xié)方差函數(shù)的變化幅度;l為兩點(diǎn)之間距離;a為常數(shù)。式(10)通過(guò)恒等變換后為

(11)

從而依據(jù)最小二乘方法可求解似高斯協(xié)方差函數(shù)待求參量a。借鑒以上分析結(jié)果,同理可求得重力-地形協(xié)方差Cgb和重力異常的自協(xié)方差Cgg,從而可獲取以海面重力信息為輸入,采用最小二乘配置方法(式(2))構(gòu)建相應(yīng)波段的海底地形模型。

2 數(shù)值分析試驗(yàn)

2.1 數(shù)據(jù)準(zhǔn)備與前期處理

選擇西太平洋某1°×1°(155.2°E—156.2°E,22.5°E—23.5°E)海區(qū)作為研究試算區(qū)域,海區(qū)輸入數(shù)據(jù)包括衛(wèi)星測(cè)高重力異常數(shù)據(jù)和稀疏船測(cè)數(shù)據(jù),數(shù)據(jù)來(lái)源和預(yù)處理如下。

(1) 海區(qū)重力異常來(lái)源于SIO(Scripps Institution of Oceanography)于2022年8月3日發(fā)布的S&S V32.1重力異常模型,該版本重力異常模型增加了Altika、Cryosat LRM、Cryosat SAR和Sentinel-3A/B等數(shù)據(jù)。研究海區(qū)重力異?？臻g分布如圖1所示,相應(yīng)的統(tǒng)計(jì)結(jié)果見(jiàn)表1。重力異常較大值集中在海區(qū)中部區(qū)域,推測(cè)該區(qū)域可能存在明顯的海山地形。

表1 研究海域重力異常(S&S V28.1)數(shù)據(jù)統(tǒng)計(jì)

圖1 研究海區(qū)重力異常Fig.1 Gravity anomaly in study sea area

(2) 稀疏船測(cè)海深數(shù)據(jù)來(lái)源于NGDC(The National Geophysical Data Center)發(fā)布的單波束和多波束船測(cè)海深數(shù)據(jù)。 共收集到1061個(gè)實(shí)測(cè)船測(cè)點(diǎn),數(shù)據(jù)頻數(shù)分布直方圖如圖2(a)所示。依圖2(a)可知,研究海區(qū)水深測(cè)量結(jié)果主要集中于5000 m以上。采用一定的處理手段對(duì)原始水深測(cè)量數(shù)據(jù)進(jìn)行質(zhì)量控制,本文處理策略是依據(jù)2019年12月發(fā)布的海洋空間格網(wǎng)大小為15″的SRTM15+V2.0作為參考模型[30],采用雙線(xiàn)性?xún)?nèi)插技術(shù)將SRTM15+V2.0模型海深值內(nèi)插到離散船測(cè)點(diǎn)處,并與船測(cè)點(diǎn)海深值作差得到海深殘差;借鑒3σ誤差處理方法,認(rèn)為海深殘差絕對(duì)值大于3倍海深殘差標(biāo)準(zhǔn)差的船測(cè)點(diǎn)是粗差點(diǎn),然后作剔除處理。最終利用3σ方法處理原始船測(cè)海深數(shù)據(jù)后,共剔除34個(gè)可疑水深點(diǎn),剩余1027個(gè)船測(cè)水深點(diǎn)。原始水深、處理后水深及水深殘差統(tǒng)計(jì)結(jié)果見(jiàn)表2,表中括號(hào)內(nèi)數(shù)值代表參與統(tǒng)計(jì)的水深點(diǎn)個(gè)數(shù)。由表2可知,以SRTM15+V2.0作為對(duì)比模型,海深殘差標(biāo)準(zhǔn)差為44.59 m,從而粗差閾值為133.77 m。原始船測(cè)水深數(shù)據(jù)剔除粗差后結(jié)果的相應(yīng)直方圖如圖2(b)所示,船測(cè)數(shù)據(jù)空間分布情況如圖3所示(圖中背景為S&S V19.1海深模型)。

表2 海深數(shù)據(jù)統(tǒng)計(jì)結(jié)果

圖2 船測(cè)數(shù)據(jù)分布直方圖Fig.2 Distribution histogram of shipborne bathymetry data

圖3 研究海區(qū)粗差剔除后船測(cè)點(diǎn)分布Fig.3 Distribution of survey points after gross error elimination

利用衛(wèi)星測(cè)高重力數(shù)據(jù)通常僅可恢復(fù)有限波段海底地形信息,剩余頻段地形信息常使用實(shí)際船測(cè)水深數(shù)據(jù)進(jìn)行補(bǔ)償; 另外為評(píng)估最終構(gòu)建的海深模型精度,通常選擇部分船測(cè)點(diǎn)作為外部檢核參考[21]。從而本文基于1027個(gè)粗差處理后的船測(cè)點(diǎn),隨機(jī)均勻選擇847個(gè)船測(cè)點(diǎn)(約占船測(cè)點(diǎn)總數(shù)的82.47%)作為控制點(diǎn)以補(bǔ)償非反演頻段地形信息;剩余180個(gè)船測(cè)點(diǎn)(約占船測(cè)點(diǎn)總數(shù)的17.53%)作為評(píng)估海深模型的外部檢核點(diǎn)。

另外,S&S V19.1海深模型中奇數(shù)海深值由實(shí)際海深約束生成,偶數(shù)海深值為重力反演預(yù)測(cè)水深值,從而可認(rèn)為這些奇數(shù)海深值是實(shí)際船測(cè)海深值。研究海區(qū)S&S V19.1海深模型中奇數(shù)(實(shí)際船測(cè)數(shù)據(jù)約束)和偶數(shù)(重力預(yù)測(cè)海深)海深值空間分布情況如圖4(a)所示,其中“·”代表重力預(yù)測(cè)海深的位置,“?”代表實(shí)際水深約束生成的海深值位置。S&S V19.1模型在試驗(yàn)海區(qū)共有616個(gè)船測(cè)實(shí)際值,相應(yīng)的統(tǒng)計(jì)結(jié)果見(jiàn)表2,海深平均值和標(biāo)準(zhǔn)差分別為4 861.00 m和816.09 m;船測(cè)數(shù)據(jù)頻數(shù)分布直方圖如圖4(b)所示。

圖4 S&S V19.1船測(cè)水深數(shù)據(jù)Fig.4 S&S V19.1 shipborne bathymetry data

對(duì)比S&S V19.1模型中616個(gè)船測(cè)點(diǎn)與筆者收集到的1027個(gè)船測(cè)點(diǎn)空間分布如圖5(a)所示,圖中“·”號(hào)表示筆者收集的船測(cè)點(diǎn)位置,“?”號(hào)代表S&S V19.1中船測(cè)點(diǎn)位置。圖5(a)船測(cè)點(diǎn)空間分布比較就總體而言,筆者收集的船測(cè)數(shù)據(jù)與S&S V19.1在大部分區(qū)域重合;部分海區(qū)S&S V19.1相比于本文增加了多條測(cè)線(xiàn)的數(shù)據(jù),如研究海區(qū)北部,S&S V19.1中包含了筆者不掌握的船測(cè)水深數(shù)據(jù)。需要說(shuō)明的是,因S&S V19.1空間網(wǎng)格大小為1′,從而基于S&S V19.1提取的船測(cè)點(diǎn)間距也是1′;而筆者從NGDC下載的聲學(xué)水深測(cè)量結(jié)果為不規(guī)則離散分布,故認(rèn)為S&S V19.1使用了本文以外的其他水深測(cè)量數(shù)據(jù),數(shù)據(jù)量更加豐富。從而為盡可能充分利用可用船測(cè)數(shù)據(jù),將S&S V19.1中包含的船測(cè)點(diǎn)(616個(gè))和本文收集的船測(cè)控制點(diǎn)(847個(gè))合并以補(bǔ)償非反演頻段海底地形信息(共1463個(gè))。最終,試驗(yàn)海區(qū)控制點(diǎn)和檢核點(diǎn)空間分布如圖5(b),其中三角形表示檢核點(diǎn)位置,白色六邊形表示控制點(diǎn)的空間位置。

圖5 船測(cè)水深點(diǎn)分布Fig.5 Distribution of shipborne bathymetry points

2.2 海底地形模型構(gòu)建

海底地形主要影響海面重力中長(zhǎng)波段信息,因而利用衛(wèi)星測(cè)高重力數(shù)據(jù)也主要恢復(fù)海底地形中長(zhǎng)波段信息。依照文獻(xiàn)[31]和筆者前期研究相關(guān)結(jié)論,本文選擇15 km和160 km作為截止波段,即依托衛(wèi)星測(cè)高重力異常數(shù)據(jù)主要恢復(fù)15～160 km波段海底地形信息,15 km以?xún)?nèi)和160 km以外的海底地形信息依托船測(cè)海深進(jìn)行補(bǔ)償。首先,利用GMT(generic mapping tools)中格網(wǎng)化操作模塊(surface)將試驗(yàn)海區(qū)離散船測(cè)控制點(diǎn)格網(wǎng)化為網(wǎng)格大小為1′的海深格網(wǎng)[32];然后,利用文獻(xiàn)[21]基于撓曲補(bǔ)償和信噪比設(shè)計(jì)的高通、低通和帶通濾波器對(duì)格網(wǎng)海深和S&S V32.1重力異常實(shí)施濾波操作,分別獲得長(zhǎng)波(波長(zhǎng)大于160 km)和短波(波長(zhǎng)小于15 km)海底地形,反演波段(15～160 km)重力異常,相應(yīng)的數(shù)據(jù)空間分布和統(tǒng)計(jì)結(jié)果如圖6和表3所示。需要說(shuō)明的是,依據(jù)海底地形解析反演算法[33]可知,重力異常經(jīng)帶通濾波并向下延拓后與海底地形表現(xiàn)出較強(qiáng)的相關(guān)關(guān)系,因而為充分利用重力場(chǎng)元與海底地形間相關(guān)關(guān)系進(jìn)行統(tǒng)計(jì)分析,圖6(c)表示的反演波段重力異常為經(jīng)帶通濾波并延拓到平均海深面的重力異常結(jié)果。

表3 數(shù)據(jù)統(tǒng)計(jì)結(jié)果

圖6 不同波段數(shù)據(jù)Fig.6 Different band data

對(duì)比圖6(a)和圖6(b)可發(fā)現(xiàn),長(zhǎng)波海底地形表現(xiàn)一種大尺度的低頻特征,而短波海底地形數(shù)據(jù)則呈現(xiàn)出高頻變化特征,從而可利用重力梯度和重力異常對(duì)不同頻段地形敏感程度進(jìn)行聯(lián)合反演[34]。依據(jù)式(2)恢復(fù)海底地形信息需解決的關(guān)鍵技術(shù)問(wèn)題是地形-重力協(xié)方差矩陣和重力協(xié)方差矩陣的獲取,其中重力異常協(xié)方差矩陣為對(duì)稱(chēng)矩陣,地形-重力協(xié)方差矩陣是非對(duì)稱(chēng)矩陣。重力異常協(xié)方差矩陣求解技術(shù)目前已十分成熟[35],而地形-重力協(xié)方差矩陣研究成果相對(duì)較少。分析地形-重力協(xié)方差矩陣元素結(jié)構(gòu),結(jié)合互協(xié)方差性質(zhì)易知,地形-重力協(xié)方差矩陣中上三角矩陣元素值可將對(duì)應(yīng)地形和重力位置關(guān)系施以地形-重力協(xié)方差函數(shù)(如本文選擇的似高斯協(xié)方差函數(shù))得到;地形-重力協(xié)方差矩陣下三角矩陣元素值則是對(duì)應(yīng)重力和地形位置關(guān)系施以重力-地形協(xié)方差函數(shù)獲取。

基于以上分析,首先以反演波段重力異常格網(wǎng)數(shù)據(jù)為輸入,統(tǒng)計(jì)分析試驗(yàn)海區(qū)重力異常分布規(guī)律,為使統(tǒng)計(jì)結(jié)果充分反映統(tǒng)計(jì)數(shù)據(jù)間相互關(guān)系,參與統(tǒng)計(jì)樣本應(yīng)盡可能豐富,從而本文實(shí)施數(shù)據(jù)統(tǒng)計(jì)分析時(shí),數(shù)據(jù)間橫向距離和縱向距離最大設(shè)置為40′。考慮到協(xié)方差函數(shù)非負(fù)性質(zhì),去掉協(xié)方差統(tǒng)計(jì)為負(fù)值結(jié)果,最終得到的重力異常協(xié)方差統(tǒng)計(jì)結(jié)果如圖7(a)所示,距離為零時(shí),協(xié)方差統(tǒng)計(jì)結(jié)果為2 721.248 2 mGal2。依據(jù)重力異常協(xié)方差統(tǒng)計(jì)結(jié)果,根據(jù)式(11),采用最小二乘方法擬合參數(shù)結(jié)果如圖7(b)所示,圖中方點(diǎn)為協(xié)方差待擬合數(shù)值分布,直線(xiàn)為擬合結(jié)果,相應(yīng)的參數(shù)計(jì)算結(jié)果為0.207 8。采用相同的統(tǒng)計(jì)策略,得到的重力-地形互協(xié)方差和地形-重力互協(xié)方差結(jié)果如圖8和圖9所示。其中距離為零時(shí),重力-地形互協(xié)方差和地形-重力互協(xié)方差統(tǒng)計(jì)值大小均為33.776 7 mGal·km,其中重力-地形互協(xié)方差似高斯函數(shù)參數(shù)擬合值為0.210 8,地形-重力互協(xié)方差似高斯函數(shù)參數(shù)擬合值為0.192 1。最終,基于似高斯協(xié)方差函數(shù)的(互)協(xié)方差函數(shù)表達(dá)式為

(12)

圖7 重力異常協(xié)方差結(jié)果Fig.7 Gravity anomaly covariance results

圖8 重力-地形協(xié)方差結(jié)果Fig.8 Cross covariance results between gravity and terrain

圖9 地形-重力協(xié)方差結(jié)果Fig.9 Cross covariance results between terrain and gravity

研究海區(qū)格網(wǎng)點(diǎn)間歐氏距離(距離單位為′)解算結(jié)果如圖10所示,對(duì)圖10距離矩陣施以式(12)重力異常似高斯協(xié)方差函數(shù)解算,可得相應(yīng)的重力異常協(xié)方差矩陣(如圖11(a))。依前面分析知,地形-重力協(xié)方差矩陣為非對(duì)稱(chēng)矩陣,從而將距離矩陣所有元素直接施以地形-重力協(xié)方差函數(shù)計(jì)算結(jié)果并非地形-重力協(xié)方差矩陣。對(duì)于地形-重力協(xié)方差矩陣求解,分為兩個(gè)步驟:①對(duì)圖10距離矩陣的上三角所有元素施以地形-重力似高斯協(xié)方差函數(shù),距離矩陣的下三角矩陣所有元素施以重力-地形似高斯協(xié)方差函數(shù);②將各自計(jì)算結(jié)果合并可得地形-重力協(xié)方差矩陣,結(jié)果如圖11(b)所示?；讷@取的試驗(yàn)海區(qū)地形-重力協(xié)方差矩陣和重力異常協(xié)方差矩陣,以反演波段重力異常為輸入,采用最小二乘配置方法可恢復(fù)相應(yīng)波段的海底地形,然后疊加依據(jù)稀疏船測(cè)數(shù)據(jù)獲取的非反演波段海底地形可獲得試驗(yàn)區(qū)域海深模型。試驗(yàn)海區(qū)海深模型結(jié)果如圖12所示,為描述方便,稱(chēng)為BAT_LSC_1模型。

圖10 距離矩陣Fig.10 Distance matrix

圖12 BAT_LSC_1海深模型Fig.12 BAT_LSC_1 bathymetry model

2.3 模型精度評(píng)價(jià)分析

引入ETOPO1和DTU18海深模型評(píng)估BAT_LSC_1海深模型效能;同時(shí)采用未考慮地殼均衡補(bǔ)償?shù)闹亓Ξ惓?dǎo)納函數(shù)理論[36]恢復(fù)目標(biāo)海區(qū)海深模型,建立的海深模型稱(chēng)為BAT_AF_1模型,研究海區(qū)各海深模型數(shù)值統(tǒng)計(jì)結(jié)果見(jiàn)表4。

表4 海深模型統(tǒng)計(jì)結(jié)果

表4中各海深模型統(tǒng)計(jì)結(jié)果顯示,研究海區(qū)海深最深約為5500 m,其中BAT_LSC_1最大海深處海深達(dá)到了6000 m;各模型海深平均值為4700 m左右;海深標(biāo)準(zhǔn)差均處于750 m左右,說(shuō)明本文利用重力異常數(shù)據(jù)反演海深在試驗(yàn)海區(qū)與目前國(guó)際發(fā)布的海深模型表現(xiàn)相似。進(jìn)一步評(píng)價(jià)本文海深模型質(zhì)量,以事先準(zhǔn)備的外部實(shí)測(cè)海深數(shù)據(jù)作為外部檢核參考,采用雙線(xiàn)性插值方法將海深模型值內(nèi)插到檢核點(diǎn)處,然后與檢核點(diǎn)處海深作差并統(tǒng)計(jì)差值結(jié)果,各海深模型檢核精度統(tǒng)計(jì)結(jié)果見(jiàn)表5,表中相對(duì)精度為檢核均方差與檢核點(diǎn)平均海深絕對(duì)值之比。

表5 海深模型檢核差值統(tǒng)計(jì)結(jié)果

表5檢核差值統(tǒng)計(jì)結(jié)果顯示,ETOPO1海深模型檢核差值最大值達(dá)到了900 m,其余海深模型檢核差值最值為500 m。就差值標(biāo)準(zhǔn)差和均方差統(tǒng)計(jì)指標(biāo)而言,DTU18在研究海區(qū)檢核精度最高,約為百米;BAT_LSC_1和BAT_AF_1檢核精度優(yōu)于ETOPO1模型,相較于ETOPO1,BAT_LSC_1檢核精度提高了38%左右;BAT_AF_1檢核精度提升約37%;BAT_LSC_1檢核精度與BAT_AF_1相當(dāng),檢核精度略高于BAT_AF_1。4種海深模型與檢核點(diǎn)相關(guān)系數(shù)均達(dá)到了0.94以上,其中BAT_LSC_1和BAT_AF_1與檢核點(diǎn)相關(guān)系數(shù)均高于0.97。ETOPO1檢核相對(duì)精度為5.73%,BAT_LSC_1和BAT_AF_1相對(duì)精度均優(yōu)于4%。綜合海深模型檢核差值統(tǒng)計(jì)分析結(jié)果,驗(yàn)證了本文方法構(gòu)建海深模型技術(shù)方案可行,即利用統(tǒng)計(jì)方法構(gòu)造地形-重力互協(xié)方差函數(shù)方法有效。進(jìn)一步比對(duì)各海深模型檢核效果,統(tǒng)計(jì)不同檢核差值(差值取絕對(duì)值)范圍內(nèi)檢核點(diǎn)所占比例,結(jié)果如圖13所示。

圖13 不同檢核差值比例Fig.13 The ratio of checkpoints in different ranges of difference

對(duì)比不同海深模型不同檢核差值范圍內(nèi)檢核點(diǎn)比例可以看出,DTU18海深模型明顯優(yōu)于其他3種海深模型。差值小于120 m左右范圍時(shí),ETOPO1模型檢核點(diǎn)數(shù)量多于BAT_LSC_1和BAT_AF_1;檢核差值大于120 m左右時(shí),相同差值情況下ETOPO1模型檢核點(diǎn)數(shù)量明顯減少,檢核點(diǎn)數(shù)量明顯小于BAT_LSC_1和BAT_AF_1;檢核差值范圍大概120～280 m時(shí),BAT_AF_1檢核點(diǎn)數(shù)量多于BAT_LSC_1,之后BAT_LSC_1檢核點(diǎn)數(shù)量略?xún)?yōu)于BAT_AF_1。經(jīng)以上分析可知,試驗(yàn)海區(qū)利用最小二乘配置方法與利用導(dǎo)納函數(shù)(未考慮地殼均衡補(bǔ)償)方法構(gòu)建的海深模型檢核結(jié)果大體相當(dāng),說(shuō)明了本文方法的實(shí)際可操作性。

2.4 陌生海區(qū)海底地形模型構(gòu)建

據(jù)以上研討分析結(jié)果可知,采用本文方法可有效構(gòu)建目標(biāo)海區(qū)海底地形數(shù)值模型,根據(jù)所建立的海底地形模型可較為清晰地反映目標(biāo)海區(qū)地形地貌特征。最小二乘配置技術(shù)恢復(fù)海底地形模型關(guān)鍵核心是構(gòu)造地形-重力協(xié)方差矩陣和重力異常自協(xié)方差矩陣。依統(tǒng)計(jì)分析理論,協(xié)方差矩陣反映了目標(biāo)海區(qū)相關(guān)數(shù)據(jù)信息分布的內(nèi)在規(guī)律,因此若某一陌生海區(qū)數(shù)據(jù)分布與已知海區(qū)數(shù)據(jù)分布具有相似特性,則可將反映已知海區(qū)數(shù)據(jù)分布內(nèi)在特性的協(xié)方差函數(shù)遷移、推廣、應(yīng)用于陌生海區(qū),進(jìn)而實(shí)現(xiàn)以重力異常為輸入的無(wú)船測(cè)數(shù)據(jù)支撐環(huán)境下,運(yùn)用最小二乘配置方法反演海底地形目的。

選擇西太平洋某5°×5°(157°E—162°E,19°N—24°N)海域作為陌生備選海區(qū),以上述1°×1°(155.2°E—156.2°E,22.5°E—23.5°E)試驗(yàn)海區(qū)為已知海區(qū)(以下稱(chēng)為區(qū)域1)。將陌生備選海區(qū)S&S V32.1重力異常作為已知數(shù)據(jù)(圖14),相應(yīng)的數(shù)據(jù)統(tǒng)計(jì)結(jié)果見(jiàn)表6,通過(guò)比較、評(píng)判已知海區(qū)和陌生海區(qū)重力異常數(shù)據(jù)相似性程度,搜尋出備選海區(qū)中與已知海區(qū)數(shù)據(jù)特征最為相近區(qū)域。本文借鑒影像紋理特征匹配相關(guān)知識(shí),首先將已知海區(qū)和備選海區(qū)重力異常作灰度化處理,得到像素值范圍在[0,255]的標(biāo)準(zhǔn)灰度影像,影像結(jié)果如圖15所示。

表6 備選海區(qū)重力異常

圖14 重力異常Fig.14 Gravity anomaly

圖15 標(biāo)準(zhǔn)灰度影像Fig.15 Standard grayscale image

為確定陌生備選海區(qū)中與已知海區(qū)(區(qū)域1)相似區(qū)域準(zhǔn)確位置,同時(shí)兼顧影像匹配效率,本文分粗匹配和精匹配兩個(gè)步驟環(huán)節(jié)在備選海區(qū)完成相似性區(qū)域搜索。首先,以已知海區(qū)重力異常灰度影像和陌生備選海區(qū)重力異?；叶扔跋駷檩斎?根據(jù)SAD粗匹配算法模型,在備選海區(qū)中搜尋與已知海區(qū)最為相似的5幅子圖,尋出的子圖位置如圖16中小圖的方框區(qū)域。然后,按照灰度共生矩陣(GLCM)和灰度-梯度共生矩陣(GGCM)生成原則,構(gòu)造已知影像和5幅子圖的GLCM和GGCM,考慮到計(jì)算消耗,將灰度值壓縮為16個(gè)灰度級(jí)(具體構(gòu)造過(guò)程見(jiàn)前文描述),分別選擇GLCM特征值中和平均、和方差、和熵和反差,GGCM特征值中灰度熵、灰度平均和灰度均方差,以及Hu矩(7個(gè)不變矩)一共14個(gè)特征值作為影像特征向量構(gòu)成元素;以設(shè)計(jì)的特征向量為參考,分別計(jì)算5幅子圖與已知影像間的歐式距離,與已知影像距離最小的子圖即為最相似海區(qū),精匹配結(jié)果如圖16(a)中藍(lán)色方框所示區(qū)域(以下稱(chēng)為區(qū)域2)。另外,采用以上粗匹配和精匹配方法獲得與已知海域相似度最低(與已知海域歐式距離最遠(yuǎn))區(qū)域如圖16(b)中藍(lán)色方框所示(稱(chēng)為區(qū)域3)。

圖16 影像匹配結(jié)果Fig.16 Image matching results

從S&S V32.1重力異常模型中提取出區(qū)域2和區(qū)域3衛(wèi)星測(cè)高重力異常數(shù)據(jù),相應(yīng)的統(tǒng)計(jì)結(jié)果見(jiàn)表7,區(qū)域2和區(qū)域3重力異常平均值分別為4.13 mGal和25.36 mGal,標(biāo)準(zhǔn)差f分別為43.20 mGal和58.96 mGal。筆者收集到區(qū)域2實(shí)測(cè)水深數(shù)據(jù)分布如圖17中紅色圓點(diǎn)所示,共計(jì)265個(gè);S&S V19.1海深模型中通過(guò)實(shí)測(cè)水深數(shù)據(jù)約束生成的模型點(diǎn)共309個(gè)(認(rèn)為S&S V19.1中實(shí)測(cè)水深值約束生成的模型點(diǎn)為實(shí)際水深測(cè)量值),分布情況如圖17所示,相應(yīng)的統(tǒng)計(jì)結(jié)果見(jiàn)表7?？梢钥闯?S&S V19.1中實(shí)測(cè)水深數(shù)據(jù)相較于筆者收集的船測(cè)數(shù)據(jù)更為豐富,因此將S&S V19.1實(shí)測(cè)水深數(shù)據(jù)作為評(píng)估海深模型精度的外部檢核參考。區(qū)域3 S&S V19.1實(shí)測(cè)水深數(shù)據(jù)統(tǒng)計(jì)結(jié)果見(jiàn)表7。

表7 區(qū)域2、3數(shù)據(jù)

圖17 S&S V19.1實(shí)際船測(cè)水深數(shù)據(jù)和NGDC船測(cè)水深數(shù)據(jù)分布比較Fig.17 Comparison of shipborne bathymetry data between S&S V19.1 and NGDC

據(jù)以上研究分析,因區(qū)域1和區(qū)域2海底地形地貌相似程度高,以區(qū)域2反演波段(15～160 km) S&S V32.1重力異常為輸入,將區(qū)域1事先解算的地形-重力協(xié)方差函數(shù)、重力-地形協(xié)方差函數(shù)和重力異常協(xié)方差函數(shù)遷移于區(qū)域2,應(yīng)用最小二乘配置技術(shù)可獲得區(qū)域2反演波段海底地形;而后以S&S V19.1作為區(qū)域2先驗(yàn)海深模型,并施以濾波操作,從而可補(bǔ)充非反演波段(波長(zhǎng)小于15 km和波長(zhǎng)大于160 km)海底地形信息;疊加區(qū)域2所有波段地形信息,進(jìn)而可實(shí)現(xiàn)陌生海域(區(qū)域2)海底地形模型構(gòu)建。采用以上海底地形構(gòu)建策略,建立的區(qū)域2海深模型(稱(chēng)為BAT_LSC_2)效果如圖18所示。

圖18 BAT_LSC_2模型Fig.18 BAT_LSC_2 model

將區(qū)域1相關(guān)協(xié)方差函數(shù)應(yīng)用于區(qū)域2,采用最小二乘配置方法恢復(fù)的區(qū)域2海底地形模型(BAT_LSC_2模型)可明顯反映區(qū)域2的海底地形地貌特征, 區(qū)域2海區(qū)海底分布著5座大型海底高山。區(qū)域2海底地形的呈現(xiàn)狀態(tài)定性表明:以海洋表面重力異常為相似性參考,將已知海區(qū)解算的協(xié)方差函數(shù)遷移、推廣并應(yīng)用于相似未知海區(qū),采用最小二乘配置方法恢復(fù)得到的海底地形可較清晰地表達(dá)未知海區(qū)海底地形特征。為進(jìn)一步評(píng)估該方法的適用性,依據(jù)重力異常導(dǎo)納函數(shù)反演理論,建立區(qū)域2海底地形模型(下文簡(jiǎn)稱(chēng)BAT_AF_2模型)。以事先設(shè)置的外部船測(cè)檢核點(diǎn)為參考,各海深模型的檢核統(tǒng)計(jì)結(jié)果見(jiàn)表8。

表8 海深模型檢核結(jié)果

表8統(tǒng)計(jì)結(jié)果表明,BAT_LSC_2檢核統(tǒng)計(jì)結(jié)果指標(biāo)優(yōu)于其余模型,如以檢核均方差為例,BAT_LSC_2檢核均方差為254.07 m,而B(niǎo)AT_AF_2、ETOPO1和DTU18分別為386.62、454.87和308.66 m。BAT_LSC_2與檢核點(diǎn)相關(guān)系數(shù)最高,達(dá)到了0.970 6。4種海深模型中,ETOPO1模型檢核相對(duì)精度超過(guò)10%,其余模型相對(duì)精度均在10%以?xún)?nèi),其中BAT_LSC_2模型相對(duì)精度為6.16%?；谝陨虾?jiǎn)要分析,進(jìn)一步定量闡明了本文將已知海區(qū)相關(guān)數(shù)據(jù)統(tǒng)計(jì)計(jì)算結(jié)果遷移、推廣并應(yīng)用于相似陌生海區(qū)方法可行,構(gòu)建的海底地形模型具有較強(qiáng)表現(xiàn)能力。比較各海深模型不同檢核差值范圍內(nèi)檢核點(diǎn)數(shù)量占檢核點(diǎn)總數(shù)情況,結(jié)果如圖19所示。

圖19 不同檢核差值范圍比例Fig.19 The ratio of checkpoints in different ranges of difference

圖19各海深模型不同檢核差值范圍檢核點(diǎn)數(shù)量占比曲線(xiàn)比對(duì)可以發(fā)現(xiàn),DTU18模型在差值小于400 m左右范圍表現(xiàn)優(yōu)于BAT_LSC_2,而B(niǎo)AT_LSC_2表現(xiàn)整體優(yōu)于BAT_AF_2和ETOPO1;BAT_AF_2和ETOPO1模型表現(xiàn)相當(dāng)。經(jīng)以上分析可知,就模型檢核統(tǒng)計(jì)各指標(biāo)和不同檢核差值范圍檢核點(diǎn)比例而言,利用最小二乘配置方法構(gòu)建的區(qū)域2海底地形模型明顯好于基于重力異常導(dǎo)納函數(shù)構(gòu)建的海底地形模型,進(jìn)一步說(shuō)明了本文解決無(wú)船測(cè)水深數(shù)據(jù)支撐環(huán)境下海底地形模型構(gòu)建具有可行性。

類(lèi)似地,將區(qū)域1解算的地形-重力協(xié)方差函數(shù)、重力-地形協(xié)方差函數(shù)和重力異常協(xié)方差函數(shù)遷移于區(qū)域3,應(yīng)用最小二乘配置技術(shù)亦可反演獲得陌生區(qū)域3海底地形信息,稱(chēng)為BAT_LSC_3模型。同理,采用導(dǎo)納函數(shù)方法獲得的區(qū)域3海底地形模型命名為BAT_AF_3模型。BAT_LSC_3模型和BAT_AF_3模型統(tǒng)計(jì)結(jié)果及其外部水深檢核統(tǒng)計(jì)結(jié)果分別見(jiàn)表9和表10。

表9 海深模型

表10 海深模型檢核結(jié)果

表9顯示BAT_LSC_3模型和BAT_AF_3模型整體表現(xiàn)基本一致。由表10可知,BAT_LSC_3模型和BAT_AF_3模型檢核均方差分別為347.03 m和359.56 m,BAT_LSC_3模型精度略?xún)?yōu)于BAT_AF_3模型。相較于BAT_LSC_2,BAT_LSC_3質(zhì)量明顯下降:①BAT_LSC_3系統(tǒng)誤差明顯增大,BAT_LSC_2檢核差值平均值為46.08 m,BAT_LSC_3檢核差值平均值達(dá)到了182.73 m;②BAT_LSC_2檢核精度明顯優(yōu)于導(dǎo)納函數(shù)結(jié)果,而B(niǎo)AT_LSC_3與導(dǎo)納函數(shù)海底地形結(jié)果比較接近。不同檢核差值范圍內(nèi)BAT_LSC_3模型和BAT_AF_3模型檢核點(diǎn)數(shù)量占檢核點(diǎn)總數(shù)情況如圖20所示。

圖20 不同檢核差值范圍比例Fig.20 The ratio of checkpoints in different ranges of difference

對(duì)比圖20和圖19可知,BAT_LSC_3模型與導(dǎo)納函數(shù)反演結(jié)果較為接近,檢核差值小于250 m左右,BAT_AF_3模型表現(xiàn)明顯優(yōu)于BAT_LSC_3;而B(niǎo)AT_AF_2模型則在整個(gè)檢核差值范圍表現(xiàn)均好于導(dǎo)納函數(shù)結(jié)果。筆者認(rèn)為原因是:使用最小二乘配置方法構(gòu)建區(qū)域2和區(qū)域3海底地形模型時(shí),均使用了區(qū)域1地形-重力協(xié)方差函數(shù)、重力-地形協(xié)方差函數(shù)和重力異常協(xié)方差函數(shù)進(jìn)行解算,而區(qū)域2和區(qū)域3分別為與區(qū)域1相似度最高和最低區(qū)域,從而區(qū)域2協(xié)方差函數(shù)與區(qū)域1更為接近,進(jìn)而區(qū)域2解算結(jié)果質(zhì)量更高。另外值得說(shuō)明的是,本文利用重力異常恢復(fù)海底地形的最小截止波段為15 km,即奈奎斯特(Nyquist)波長(zhǎng)為15 km;則依據(jù)頻率域采樣定理可知,恢復(fù)的海底地形分辨率最高為30 km,而本文構(gòu)建的海底地形模型在波長(zhǎng)小于15 km部分利用了S&S V32.1模型補(bǔ)充,因此構(gòu)建的海底地形模型最終分辨率依賴(lài)S&S V32.1模型,海底地形網(wǎng)格大小為1.8 km左右。

3 結(jié) 論

本文嘗試探索基于已知海區(qū)數(shù)據(jù)統(tǒng)計(jì)信息恢復(fù)未知陌生海區(qū)海底地形的手段途徑。首先,利用已知海區(qū)可用的稀疏船測(cè)水深數(shù)據(jù)和衛(wèi)星測(cè)高重力異常數(shù)據(jù),統(tǒng)計(jì)建立了數(shù)據(jù)已知海區(qū)(區(qū)域1)地形-重力互協(xié)方差函數(shù)和重力自協(xié)方差函數(shù),依據(jù)最小二乘配置方法恢復(fù)了區(qū)域1的海底地形;然后,以衛(wèi)星測(cè)高海面重力異常作為相似性參考,采用影像粗匹配和精匹配算法模型,尋出了與已知海區(qū)重力異常最為相似的無(wú)船測(cè)水深數(shù)據(jù)海區(qū)(區(qū)域2);最后,將區(qū)域1相關(guān)數(shù)據(jù)的協(xié)方差函數(shù)遷移、推廣、應(yīng)用于區(qū)域2,采用最小二乘配置方法構(gòu)建了相應(yīng)的海底地形模型,并評(píng)估了所構(gòu)建的海底地形模型效能。經(jīng)過(guò)定性和定量?jī)蓚€(gè)層面比對(duì)分析試驗(yàn)結(jié)果可以發(fā)現(xiàn):本文方法反演的海底地形模型可有效描述研究海區(qū)海底地形地貌特點(diǎn),反演的研究海區(qū)海深模型整體呈現(xiàn)效果好于根據(jù)重力異常導(dǎo)納函數(shù)構(gòu)建的海深模型,同時(shí)本文方法可巧妙回避導(dǎo)納函數(shù)涉及的地球物理參數(shù)獲取難題。試驗(yàn)區(qū)域1中,BAT_LSC_1模型檢核精度相較于ETOPO1提升了約38%;試驗(yàn)區(qū)域2中,無(wú)論檢核點(diǎn)差值統(tǒng)計(jì)結(jié)果,還是不同檢核差值范圍檢核點(diǎn)比例情況,BAT_LSC_2表現(xiàn)均明顯優(yōu)于BAT_AF_2和ETOPO1。