郭金運，魏志杰，祝程程，吳云龍，紀 兵

(1.山東科技大學(xué) 測繪與空間信息學(xué)院，山東 青島 266590； 2.中國地震局地震研究所 中國地震局地震大地測量重點實驗室，湖北 武漢 430071； 3.海軍工程大學(xué) 導(dǎo)航工程系，湖北 武漢 430033)

在海洋科學(xué)中，海底地形數(shù)據(jù)分析可為海洋地質(zhì)學(xué)、海洋地球物理學(xué)、海洋生態(tài)保護等提供可靠的研究資料。自20世紀70年代以來，衛(wèi)星測高技術(shù)經(jīng)過50多年的發(fā)展，其技術(shù)性能愈發(fā)成熟，測高精度和分辨率有了很大提升，海面高精度不斷提高[1]，海洋重力場研究不斷發(fā)展。1978年，Rummel等[2]利用逆Stokes公式通過大地水準面高確定了空間點的重力異常。最小二乘配置法(least squares collocation, LSC)作為大地測量學(xué)的經(jīng)典方法[3]，得到不斷改進和完善，現(xiàn)已成為利用衛(wèi)星測高數(shù)據(jù)恢復(fù)重力場模型的重要方法。此外，逆Vening-Meinesz公式法[4]也是重力場反演的一種算法。隨著發(fā)射衛(wèi)星的增多以及衛(wèi)星數(shù)據(jù)質(zhì)量的提高，聯(lián)合多衛(wèi)星數(shù)據(jù)資料提高重力場反演精度[5-7]逐漸成為一種普遍的方法。

在重力場精度逐漸提高過程中，重力場反演海底地形的理論方法不斷發(fā)展。目前常見的反演海底地形的方法有重力地質(zhì)法[8]、最小二乘配置法[9]和導(dǎo)納函數(shù)法[10-12]等，導(dǎo)納函數(shù)法相對于前兩種方法的優(yōu)點是利用重力場和海底地形的導(dǎo)納函數(shù)關(guān)系，可以很好地應(yīng)用于有觀測重力數(shù)據(jù)的海域，體現(xiàn)了衛(wèi)星測高數(shù)據(jù)反演重力場的優(yōu)勢，大大提升了反演海底地形的效率。導(dǎo)納函數(shù)不僅包含地形起伏信息，還包括巖石圈的均衡狀態(tài)信息，為了提高海底地形反演的精度，需考慮巖石圈的均衡補償效應(yīng)[13]。地形反演過程對重力異常數(shù)據(jù)質(zhì)量要求很高，一般需要將海面上的重力異常延拓到平均水深面?？焖俑道锶~變換(fast Fourier transform, FFT)是一種常用方法，向上延拓是穩(wěn)定的過程，但向下延拓是誤差放大的過程，結(jié)果不穩(wěn)定。因迭代延拓法具有延拓深度范圍大、結(jié)果穩(wěn)定性強的特點，迭代法向下延拓開始得到應(yīng)用，但對重力異常的迭代延拓法在地形反演研究的應(yīng)用很少提及。

本研究以南海海域中沙群島南部海域(114°E～118°E，12°N～15°N)為試驗區(qū)，利用最小二乘配置法得到HY-2A重力異常，通過迭代延拓法得到延拓面上的重力異常，確定重力異常和海底地形之間相干性大于0.5的波長范圍為20～115 km。采用CRUST 1.0全球地殼模型計算得到平均地殼厚度、地幔密度和有效彈性厚度等參數(shù)，進而求得試驗區(qū)域均衡補償模型。利用重力異常數(shù)據(jù)反演得到未加補償?shù)膶?dǎo)納水深模型和均衡補償導(dǎo)納水深模型，并將兩個水深模型與ETOPO1水深模型和DTU10水深模型進行精度對比分析，討論引起誤差的原因。

1 原理與方法

1.1 迭代延拓法原理

重力異常向上延拓是收斂的過程，結(jié)果具有唯一性，而向下延拓是發(fā)散過程，結(jié)果穩(wěn)定性差，誤差較大。圖1為海平面和延拓面的位置。

圖1 海平面和延拓面的位置

假設(shè)海平面為S1=h，延拓面為S0=0，延拓高度為h=S1-S0。在S0面上重力異常的信息為p0(x,y,0),通過FFT得到P0(kx,ky,0)=F[P0(x,y,0)]，kx、ky是各自方向的波數(shù)。然后通過逆FFT延拓公式得到對應(yīng)S1面的重力異常信息，公式為：

(1)

其中F-1為逆傅里葉變換。當(dāng)h>0時，向上延拓；當(dāng)h<0時，向下延拓。采用迭代法通過向上延拓計算得到向下延拓面的重力異常[14]，步驟如下：

1) 將海平面S1上已知點P的重力異常Pg0作為對應(yīng)延拓面S0上重力異常的初始值；

2) 通過公式(1)計算海平面上的重力異常值Pgi，根據(jù)Pg0與Pgi的差值對Pg0進行改正；

3) 反復(fù)計算迭代，當(dāng)計算值與初始值在誤差允許范圍時，即可得到延拓面S0上的重力異常。

1.2 導(dǎo)納函數(shù)理論

“導(dǎo)納”描述了海洋重力異常和海底地形在頻率域上的頻譜關(guān)系。Parker[15]在1973年將FFT引入位場理論，提高了密度界面反演的計算速度。在頻率域內(nèi)，重力異常和海底地形之間的關(guān)系為[16]:

(2)

一般地，巖石圈對海底地形載荷的均衡響應(yīng)是通過區(qū)域補償機制來實現(xiàn)的，巖石圈強度決定了撓曲量的大小。當(dāng)海底地形載荷的波長小于巖石圈彈性厚度時，公式(2)等式右邊第一項起主要作用，則公式(2)簡化為：

G(k)=Z(k)H(k)，

(3)

其中：G(k)為重力異常的傅里葉變換；H(k)為海深的傅里葉變換；導(dǎo)納函數(shù)為

Z(k)=2πG(ρc-ρw)exp(-2πkd)。

(4)

公式(4)為未考慮未補償模型的導(dǎo)納函數(shù)式。

1.3 撓曲均衡補償導(dǎo)納理論

實際問題中，Z(k)不僅包含海底地形信息，還包含巖石圈的均衡狀態(tài)信息，只有確定巖石圈撓曲產(chǎn)生的重力異常，才能確定Z(k)。由于無法直接觀測到巖石圈發(fā)生的撓曲情況，不能直接計算巖石圈的撓曲所產(chǎn)生的重力異常。為研究方便，可通過撓曲均衡模型來研究巖石圈的內(nèi)部信息。當(dāng)海底地形載荷的波長與巖石圈撓曲波長相近或者更長時，巖石圈在海底地形載荷作用下發(fā)生撓曲。根據(jù)撓曲均衡模型[13]，均衡補償物質(zhì)產(chǎn)生的重力異常與海底地形之間的關(guān)系為：

G(k)=2πG(ρc-ρw)Φ′e(k)exp(-2πk(d+t))H(k)，

(5)

(6)

其中：v為巖石圈的泊松比；E為彈性模量。結(jié)合公式(3)、(5)，整理得到在撓曲均衡模型下地形和補償物質(zhì)產(chǎn)生的重力異常為：

G(k)flex=2πG(ρc-ρw)e-2πkd(1-Φ′e(k)e-2πkt)H(k)=Z(k)flexH(k) ，

(7)

其中，重力異常的均衡補償導(dǎo)納函數(shù)

(8)

由公式(8)知，撓曲均衡補償導(dǎo)納函數(shù)與地殼平均厚度t、海底洋殼密度ρc、巖石圈撓曲剛度D以及彈性厚度Te等有關(guān)。

1.4 相干性分析

(9)

其中：*表示復(fù)數(shù)共軛；〈〉表示括號內(nèi)功率譜密度乘積的實部在圓環(huán)區(qū)域內(nèi)(ki-Δk≤k≤ki+Δk)的平均值。

2 試驗數(shù)據(jù)及結(jié)果分析

2.1 研究區(qū)域

南海盆地的中央?yún)^(qū)域為深海平原，平均水深約為4 100 m，將114°E～118°E、12°N～15°N作為試驗區(qū)域。在頻率域上為消除計算過程中邊緣效應(yīng)的影響，實際計算時經(jīng)緯度分別向外擴展2°，然后對計算結(jié)果進行處理，得到研究區(qū)域的海底地形。

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

2.2.1 重力異常數(shù)據(jù)

各鎮(zhèn)、各單位要利用多種形式，積極宣傳農(nóng)田林網(wǎng)建設(shè)與改善生態(tài)環(huán)境、提高人民群眾生活質(zhì)量的關(guān)系，與抵御自然災(zāi)害、改善農(nóng)業(yè)生產(chǎn)條件的內(nèi)在聯(lián)系。積極宣傳綠化先進典型，曝光毀林案件，通過一系列宣傳教育，不斷提高廣大干部群眾的綠化、美化意識，形成“植綠、愛綠、護綠、興綠”和“保護生態(tài)、綠化家園”的文明行為，進一步增強農(nóng)民群眾建設(shè)農(nóng)田林網(wǎng)的熱情，形成全民動員、全民參與、全民動手，推進農(nóng)田林網(wǎng)建設(shè)的濃烈氛圍。

根據(jù)2011年8月發(fā)布的HY-2A衛(wèi)星大地測量任務(wù)(geomatic mission, GM)的測高數(shù)據(jù)，獲得南海海域(105°E～122°E, 2°N～26°N)1′×1′網(wǎng)格的重力異常。2016年3月30日進行變軌后開始GM，GM的周期從14天增加到168 天，工作范圍為81°S～81°N。選取時間2016年3月30日—2018年8月22日，將采樣頻率1 Hz的HY-2A/GM 衛(wèi)星高度計數(shù)據(jù)Level 2 Plus (L2P)產(chǎn)品[17]作為研究數(shù)據(jù)。首先，對HY-2A/GM 的海面高度進行校正和粗差消除預(yù)處理；然后，利用預(yù)處理的海面高來計算沿軌的大地水準面梯度，通過移去EGM2008的大地水準面梯度獲得殘余大地水準面梯度；最后，通過計算窗口半徑為0.5°的最小二乘配置法，從殘余大地水準面梯度求得1′×1′網(wǎng)格上的殘余重力異常，通過恢復(fù)EGM2008的重力異常，從殘余重力異常中獲得最終的重力異常模型[18]，試驗區(qū)的重力異常如圖2(a)所示。

2.2.2 ETOPO1水深模型和DTU10水深模型

ETOPO1水深模型是由美國國家海洋與大氣管理局(National Oceanic and Atmospheric Administration, NOAA)的國家環(huán)境信息中心(National Centers for Environmental Information, NCEI)提供的分辨率為1′的全球地形數(shù)據(jù)模型[19]，集成了來自多個國際組織的船測水深數(shù)據(jù)、測高重力數(shù)據(jù)、海岸線數(shù)據(jù)和陸地高程數(shù)據(jù)的評估和匯總。ETOPO1模型是海洋地球物理研究中使用最廣泛的數(shù)據(jù)之一。

DTU10水深模型來自丹麥科技大學(xué)(Technical University of Denmark, DTU)發(fā)布的分辨率為1′的全球地形數(shù)據(jù)模型[20]，通過ERS-1和GEOSAT大地測量任務(wù)的高度計數(shù)據(jù)計算得到全球地形模型數(shù)據(jù)。

2.2.3 船測水深數(shù)據(jù)

船載水深數(shù)據(jù)從NCEI官網(wǎng)下載得到[21]，時間跨度為1960—2016年。NCEI控制數(shù)據(jù)質(zhì)量并將合格的數(shù)據(jù)收集到數(shù)據(jù)庫中。盡管總體數(shù)據(jù)質(zhì)量準確可靠，但一些早期數(shù)據(jù)存在較大的測量誤差，因此有必要剔除這些誤差數(shù)據(jù)。將船測水深數(shù)據(jù)減去ETOPO1數(shù)據(jù)，差值大于船測水深數(shù)據(jù)的點將被刪除。圖2(b)為試驗區(qū)的船測軌跡，黑色點表示試驗控制點(13 934控制點)，白色點表示水深反演結(jié)果的檢核點(9 464個)，背景為ETOPO1海深模型。

圖2 重力異常圖和船測水深軌跡分布圖

2.2.4 CRUST 1.0全球地殼模型

CRUST 1.0為全球三維地殼模型，其中的海水深度和地形數(shù)據(jù)是從ETOPO1地形模型中提取的，并通過平均處理得到間隔1°的網(wǎng)格數(shù)據(jù)。官方下載網(wǎng)址為：http://igppweb.ucsd.edu/～gabi/crust1.html[22]。

2.3 迭代延拓法延拓結(jié)果

將HY-2A重力異常數(shù)據(jù)向下延拓到平均水深面得到迭代次數(shù)和誤差的關(guān)系，如圖3所示。

圖3 迭代次數(shù)和延拓誤差的關(guān)系

2.4 確定海底地形反演波段范圍

根據(jù)公式(9)對重力數(shù)據(jù)和海底地形在頻率域進行相干性分析。信噪比(signal noise ratio, SNR)描述信號與噪聲的比例關(guān)系：

(10)

其中：當(dāng)相干coh=0.5時，SNR為1；當(dāng)coh>0.5時，信噪比大于1。

試驗區(qū)內(nèi)重力異常數(shù)據(jù)和海底地形的相干性如圖4所示。由圖4看出，在波長范圍20～115 km內(nèi)，重力異常與海底地形的相干性>0.5，選取該波段進行海底地形反演。

圖4 重力異常和海深的相干性

2.5 撓曲均衡補償模型參數(shù)的確定

在試驗區(qū)域內(nèi)通過均衡補償模型計算得到地球物理參數(shù)如表1所示。

表1 撓曲均衡補償模型所需要的參數(shù)

通過CRUST 1.0全球地殼模型獲取海水密度ρw，海底洋殼密度ρc，地幔密度ρm，平均地殼厚度t，巖石圈泊松比v和彈性模量E。平均海水深度由NCEI發(fā)布的船測水深數(shù)據(jù)獲得。

在實際計算時，通過分析觀測導(dǎo)納和理論導(dǎo)納曲線的走勢來確定巖石圈的有效彈性厚度Te，并通過重力異常和海深分析二者的“觀測導(dǎo)納”關(guān)系，公式如下：

(11)

當(dāng)Te分別取2、6、10 km時，導(dǎo)納函數(shù)Z(k)flex的曲線變化如圖5所示。低頻部分曲線主要受有效彈性厚度Te的影響，反映了深層地殼物質(zhì)對重力場的影響；當(dāng)Te取6 km時，低頻部分“觀測導(dǎo)納”和“理論導(dǎo)納”曲線走勢最接近，因此確定研究區(qū)域的彈性厚度為6 km。

圖5 不同有效彈性厚度下“理論導(dǎo)納”和“觀測導(dǎo)納”的曲線

2.6 海底地形反演

采用移去恢復(fù)法構(gòu)建海底地形，在未補償條件和均衡補償條件下，通過導(dǎo)納函數(shù)法反演海底地形的步驟如下：

1) 將試驗區(qū)的船測水深數(shù)據(jù)進行網(wǎng)格化，得到網(wǎng)格化水深模型；

2) 將網(wǎng)格化的水深進行二維傅里葉變換，通過濾波處理得到短波段(<20 km)海深模型和長波段(>115 km)海深模型；

3) 重力異常數(shù)據(jù)經(jīng)帶通濾波處理后，在未補償條件和均衡補償條件下，分別利用迭代延拓法向下延拓得到延拓面的重力異常，在波段20～115 km波長范圍內(nèi)反演海深；

4) 將各個波段的海深相加，得到未補償條件下的海深模型和均衡補償條件下的海深模型。將未補償導(dǎo)納水深模型記為Model1，均衡補償導(dǎo)納水深模型記為Model2。

2.7 試驗結(jié)果

圖6顯示了試驗區(qū)域反演的兩個海深模型Model1、Model2與ETOPO1水深模型和DTU10水深模型的結(jié)果，可以看出試驗區(qū)域整體的地形分布，4個模型都描繪了試驗區(qū)域的主要地形。表2為格網(wǎng)化的NCEI船測水深數(shù)據(jù)、Model1、Model2、ETOPO1水深模型和DTU10水深模型的數(shù)據(jù)統(tǒng)計。NCEI船測水深數(shù)據(jù)、Model1、Model2和ETOPO1水深模型的平均水深相差不大，Model2相對于Model1精度有所提升，而DTU10水深模型的海深平均值相對于其他數(shù)據(jù)海深平均值偏差較大，誤差超過100 m。

(a)Model1；(b)Model2；(c)ETOPO1水深模型；(d)DTU10水深模型圖6 海底地形模型

表2 NCEI、Model1、Model2、ETOPO1和DTU10水深數(shù)據(jù)統(tǒng)計表Tab. 2 Data statistics of NCEI, Model1, Model2, ETOPO1 and DTU10 bathymetry m

2.8 精度分析

為進一步分析各模型精度，將4個海深模型在NCEI船測水深檢核點上進行模型插值計算得到相應(yīng)的海深值，通過與NCEI船測水深數(shù)據(jù)作差比較各自的模型精度。表3給出了各海深模型在檢核點的差異統(tǒng)計結(jié)果。本研究反演的兩個水深模型差異的平均值分別為2.5、-1.4 m，ETOPO1水深模型與反演的兩個水深模型平均海深差異接近，但DTU10水深模型差異的平均值則相對較大。Model1、Model2和ETOPO1水深模型的相對誤差接近，而DTU10水深模型相對誤差達4.48%。按照標(biāo)準差質(zhì)量統(tǒng)計來看，水深模型精度從高到低依次為Model2、Model1、ETOPO1水深模型和DTU10水深模型(標(biāo)準差分別為150.8、 159.7、164.3、196.1 m)。由表3可知，通過迭代延拓法處理得到的重力異常，經(jīng)反演計算得到的兩個水深模型精度均優(yōu)于ETOPO1水深模型和DTU10水深模型。

表3 Model1、Model2、ETOPO1和DTU10數(shù)據(jù)在檢核點上誤差的統(tǒng)計結(jié)果Tab. 3 Data statistical results of error of Model1, Model2, ETOPO1 and DTU10 at check points

圖7分別顯示了Model1、Model2水深差分布直方圖，海深差在50 m范圍內(nèi)分別占51.7%和62.3%，海深差在200 m范圍內(nèi)分別占90.6%和91.9%。圖8分別顯示了兩個反演海深模型與NCEI船測水深數(shù)據(jù)的相對誤差和海深差分布關(guān)系，結(jié)果沒有明顯的差異特征，但在大于-3 000 m的深海區(qū)域反演海深精度相對較高。圖9顯示了Model1、Model2與NCEI船測水深數(shù)據(jù)相對誤差的空間位置分布。

圖7 海深差統(tǒng)計直方圖

圖8 相對誤差和海深差隨海深變化的情況

圖9 相對誤差的空間位置分布情況

通過上述統(tǒng)計結(jié)果，結(jié)合圖6和圖9海底地形模型分析可知，在中南暗沙附近(114.75°E～115.75°E，13.25°N～14.5°N)海底地形起伏較大，反演的水深精度相對較差，最大海深差可達2 000 m以上，ETOPO1水深模型和DTU10水深模型也存在相同的情況，這是地形變化劇烈影響的結(jié)果，是目前海底地形模型存在的問題。Model1在整體精度上優(yōu)于Model2和ETOPO1模型，但Model1和Model2水深模型相對于ETOPO1水深模型在水深刻畫上存在不足。

通過上述分析，認為影響水深反演的可能原因為：①船測水深數(shù)據(jù)時間跨度大，由于時效性問題前期測量的數(shù)據(jù)精度有待提高，粗差很難完全剔除，可能導(dǎo)致在分析海底地形與重力異常的導(dǎo)納性上出現(xiàn)誤差；②理論上海底地形與重力數(shù)據(jù)在頻率域上存在良好的線性關(guān)系，但實際上海底地形復(fù)雜，并且在頻率域上求二者相干性時用到的圓環(huán)內(nèi)的平均值也會產(chǎn)生誤差，導(dǎo)致反演水深時產(chǎn)生誤差；③將船測水深數(shù)據(jù)和離散數(shù)據(jù)網(wǎng)格化導(dǎo)致格網(wǎng)數(shù)據(jù)精度不高，得到的水深模型也會受到影響；數(shù)據(jù)在濾波過程中有信號損失，是導(dǎo)致精度損失的原因之一。

3 結(jié)論

通過重力導(dǎo)納理論，利用HY-2A重力異常延拓數(shù)據(jù)反演得到未補償導(dǎo)納水深模型Model1和均衡補償導(dǎo)納水深模型Model2，并將反演結(jié)果與目前公認的ETOPO1水深模型和DTU10水深模型進行外部檢核比較，得到如下結(jié)論：

1) 利用迭代延拓法將HY-2A重力數(shù)據(jù)向下延拓，得到平均水深面上的重力異常，其精度誤差在0.1 mGal以內(nèi)。

2) 通過實際計算得到的“觀測導(dǎo)納”和均衡補償理論計算的“理論導(dǎo)納”進行比較，可以看出低頻部分主要受地殼有效彈性厚度的影響，最后確定試驗區(qū)域內(nèi)的有效彈性厚度為6 km。

3) 雖然整體上Model2精度略高于Model1，標(biāo)準差提升約9 m，但精度提升效果不明顯，說明有效彈性厚度為6 km時，均衡補償模型對試驗區(qū)地形的精度改善效果有限。

4) 利用迭代延拓法求得延拓面上的重力異常，經(jīng)反演計算得到的兩個水深模型精度均優(yōu)于ETOPO1水深模型和DTU10水深模型，并且Model2優(yōu)于Model1。通過誤差分析，在地形變化劇烈的區(qū)域Model1、Model2、ETOPO1水深模型和DTU10水深模型相對于船測水深數(shù)據(jù)精度較差，這也是目前水深模型存在的問題。