徐 煥，于錦海，安 邦，萬曉云

1. 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100049; 2. 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京 100049; 3. 中國地質(zhì)大學(xué)(北京)土地科學(xué)與技術(shù)學(xué)院，北京 100083

海底地形和海底基巖密度可以反映出海底構(gòu)造演化活動(dòng)[1]，可幫助探明海底礦產(chǎn)資源[2]，此外，海底地形對于海事搜救等活動(dòng)也是非常重要的，因此研究海底地形不僅有理論意義，而且有實(shí)用價(jià)值。目前，利用多波束回聲測量裝置直接測量海深是最好的方法，但是由于成本高、耗時(shí)長，因此，測得的海深數(shù)據(jù)比較稀疏，難以覆蓋整個(gè)海洋[3]。近年來，利用海洋重力數(shù)據(jù)反演海底地形的工作日益受到重視，其原因是成本低、覆蓋范圍廣。從目前已有的研究來看，反演海底地形主要利用重力異常和垂直重力梯度異常。

重力地質(zhì)法(gravity-geologic method，GGM)[4]，其原理是將海面上重力異常中的短波部分與海底地形存在著線性關(guān)系，然后在實(shí)際使用時(shí)，需要結(jié)合船測的海深數(shù)據(jù)，擬合出該線性關(guān)系，從而得到海底地形。GGM方法目前已經(jīng)被廣泛地應(yīng)用于海底地形的反演[5-8]。事實(shí)上，GGM方法在海底較為平坦的區(qū)域是非常有效的，但對起伏較大的海底來說，反演海底地形的精度較差[9-10]，因此GGM方法在應(yīng)用上有一定的局限性。

頻域方法源自Parker將引力位表達(dá)式進(jìn)行傅里葉展開[11]，給出了海面重力數(shù)據(jù)與海深之間的頻率域表達(dá)式，略去海深的高階小量后，即得頻域內(nèi)的線性關(guān)系，這為利用重力異常直接反演海底地形提供了理論依據(jù)。當(dāng)考慮地殼的均衡補(bǔ)償，Parker理論得到進(jìn)一步完善，基于該理論，文獻(xiàn)[12]反演皇帝海山鏈的海深，由于地形變化劇烈，海山鏈沿線附近精度不甚理想。文獻(xiàn)[13]聯(lián)合多源重力數(shù)據(jù)反演菲律賓海域的海底地形，精度提高30%左右。此外，文獻(xiàn)[14]證實(shí)了在50～300 km波段范圍內(nèi)，海底地形與衛(wèi)星測高數(shù)據(jù)之間存在高度相關(guān)性。在文獻(xiàn)[12,14]的研究基礎(chǔ)上，文獻(xiàn)[15]說明了在15～160 km的波段內(nèi)，海底地形和重力異常間具有線性相關(guān)的特征。因?yàn)樵趹?yīng)用Parker理論反演海底地形時(shí)都會(huì)舍去海深二階以上的量，所以精確度存在不確定性[16]。例如，文獻(xiàn)[17]對比了GGM和Parker理論反演海底地形的結(jié)果后指出，GGM方法在可靠性和精度方面具有明顯的優(yōu)勢。

作為重力數(shù)據(jù)的延伸，垂直重力梯度(vertical gravity gradient,VGG)對海底地形高頻部分的敏感度超過重力本身[18-19]。文獻(xiàn)[20]推導(dǎo)出長方體產(chǎn)生的VGG公式，該公式已在地球物理勘探方面得到了許多應(yīng)用。例如，文獻(xiàn)[21]利用重力梯度全張量數(shù)據(jù)，采用共軛梯度法反演地質(zhì)體與密度，但其僅僅限于模擬計(jì)算。文獻(xiàn)[22]采用聚焦反演方法來預(yù)測異常地質(zhì)體的埋深信息。此外，文獻(xiàn)[23]也利用該公式，計(jì)算了不規(guī)則地質(zhì)異常體產(chǎn)生的垂直重力梯度，并將其用于評估西墨西哥Jasper海山的密度。在將VGG應(yīng)用于預(yù)測海深方面，文獻(xiàn)[24]利用SIO(scrioos institude of oceangraphy)的VGG數(shù)據(jù)，并結(jié)合船測數(shù)據(jù)來反演海底地形，例如，在南印度洋海域的反演結(jié)果與船測深度的均方根誤差約為188 m。文獻(xiàn)[25]利用VGG異常數(shù)據(jù)，顧及了Park理論中的二階量影響和采用模擬退火法反演了海底地形，并給出其精度提高了22%的結(jié)論。

綜合上述各種方法，其本質(zhì)上還是利用經(jīng)驗(yàn)公式來建立起重力數(shù)據(jù)與海底地形之間的線性關(guān)系，而這樣的線性關(guān)系是隨著海域的不同而變化的。為避開使用經(jīng)驗(yàn)公式帶來的影響，本文建立了一種解析的利用VGG數(shù)據(jù)來反演海底地形的方法，其原理是使用矩形分割將研究區(qū)域進(jìn)行格網(wǎng)化，然后利用各個(gè)子矩形柱體產(chǎn)生的VGG表達(dá)式建立起各個(gè)分割子域的海深與該海域上VGG異常之間的函數(shù)關(guān)系，即導(dǎo)出了關(guān)于各個(gè)分割子域海深的觀測方程組。但是在求解該觀測方程組時(shí)，研究區(qū)域邊界附近海深的求解易于受到區(qū)域外VGG異常的影響，即會(huì)產(chǎn)生邊界效應(yīng)，因此，如何評估邊界效應(yīng)的影響以及分離該影響的方法也就成了十分重要的研究課題。

為了減少邊界效應(yīng)的影響，本文將研究海域適當(dāng)進(jìn)行擴(kuò)充，然后將擴(kuò)充海域內(nèi)的海深同時(shí)進(jìn)行求解。為避免觀測方程組出現(xiàn)奇異性，在求解觀測方程組時(shí)引入了正則化參數(shù)，通過仿真計(jì)算后，確認(rèn)正則化參數(shù)選取為10-3時(shí)，反演的海域下方海深具有較好的精度。在對實(shí)際海域的海底地形進(jìn)行反演計(jì)算時(shí)，除了研究海域下方海底地形會(huì)對VGG產(chǎn)生影響外，其他的海底地形也會(huì)對VGG異常產(chǎn)生影響，本文將擴(kuò)充海域外的海底地形對VGG的影響統(tǒng)稱為遠(yuǎn)區(qū)影響。為了消除遠(yuǎn)區(qū)影響，本文將其視為干擾常數(shù)ε，并在求解觀測方程組時(shí)對其同步進(jìn)行求解。

1 理論與方法

1.1 長方體產(chǎn)生的垂直重力梯度

假設(shè)有一個(gè)海底長方體海山Ω(圖1)，其長、寬、高分別為2a、2b、H-h，其中H為Ω底部至海平面的高度，h為Ω頂部至海平面的高度(即Ω的海深)。記R是Ω在海面對應(yīng)的海域， 顯然確定Ω的形狀僅需求解h。假設(shè)P是在海平面上任意一點(diǎn)，下面計(jì)算Ω在P處產(chǎn)生的VGG。

圖1 長方體海山 Fig.1 Schematic diagram of synthetic cylindrical seamount

建立圖1所示的坐標(biāo)系，則Ω在P點(diǎn)產(chǎn)生的重力位為

(1)

式中，G是萬有引力常數(shù)(6.672×10-11Nm2kg-2)；ρ=2700 kg/m3是海山Ω的密度；(x1,y1,z1)為積分變量。因此，海山Ω在P點(diǎn)產(chǎn)生的VGG為

Gρ[IR(H,x,y)-IR(h,x,y)]

(2)

式中

(3)

積分區(qū)間R={(x1,y1);-a≤x1≤a,-b≤y1≤b}為矩形區(qū)間，進(jìn)一步計(jì)算后可得IR(H,x,y)的表達(dá)式為

(4)

對于一般的矩形積分區(qū)間R={(x1,y1);-a≤|x1-x0|≤a,-b≤|y1-y0|≤b}，結(jié)合式(4)，IR(H,x,y)可寫為

IR(H,x,y)=I(a,b,H,x-x0,y-y0)

(5)

同理，式(2)中的IR(h,x,y)與IR(H,x,y)有相同的表達(dá)形式。因此，得到了矩形海山Ω在海平面上任何一點(diǎn)P產(chǎn)生垂直重力梯度的表達(dá)式。

1.2 觀測方程與可解性分析

為了討論方便，下面恒設(shè)R={(x,y);-a≤x,y≤a}是海面上一個(gè)正方形區(qū)域，Ω是R下方對應(yīng)的海山，其密度為常數(shù)ρ，形狀為Ω={(x,y,z);(x,y)∈R,H-h(x,y)≤z≤H}，其中h(x,y)是Ω的頂部至海面的高，因此，要確定R下方的海底地形就是確定函數(shù)h(x,y)。本節(jié)將在R上VGG異常僅是由Ω產(chǎn)生的前提下來研究如何求解h(x,y)。

選擇一個(gè)步長t，將區(qū)間[-a,a]進(jìn)行2N等分，即a=N·t，這樣就能將區(qū)域R進(jìn)行格網(wǎng)化，其中典型的子格網(wǎng)區(qū)域是Rij=[xi,xi+1]×[yj,yj+1]，xi=it,yj=jt，i,j=-N,…,0,…N-1。若步長t足夠小，在Rij中h(x,y)可以被看成是常數(shù)hij，于是，利用式(4)、式(5)和疊加原理可知，由Ω生成的VGG在海面P(x,y)處的值Γ(x,y)為

Γ(x,y)=

(6)

在式(6)中去掉海水的影響，則得Ω在海面上P處產(chǎn)生的VGG異常為

δΓ(x,y)=

(7)

式中，Δρ=ρ-ρ0，而ρ0表示海水密度，在后面的計(jì)算中總是取ρ0=1030 kg/m3。

如果事先能夠給出海域R上的VGG異常的觀測值，例如在R的所有格網(wǎng)點(diǎn)上VGG異常是已知的，可得

δΓ(xp,yq)=G·Δρ[IR(H,xp,yq)-

(8)

-δΓ(xp,yq)+GΔρIR(H,xp,yq)-

(9)

式(9)就是線性化后的海底深度與VGG異常之間的觀測方程組。

圖2 模擬20 km×20 km海底地形及對應(yīng)格網(wǎng)Fig.2 The simulated 20 km×20 km seafloor topography and corresponding grid diagram

1.3 觀測方程組抗誤差干擾分析

圖3 抗隨機(jī)誤差和系統(tǒng)誤差的曲線圖Fig.3 The graph of the RMS error of the random error and the system error

2 邊界效應(yīng)的處理

上文建立的關(guān)于海深的觀測方程式(8)或式(9)都是在海底僅有一個(gè)異常海山的前提下導(dǎo)出的，當(dāng)處理實(shí)際海域時(shí)，顯然海域下方的情況要復(fù)雜得多，特別是還存在著其他的異常海山。例如，圖4(a)模擬了40 km×40 km海域D下方的海山，如果僅需要求解中心20 km×20 km區(qū)域(稱為內(nèi)區(qū)R)下方的海底地形，顯然在內(nèi)區(qū)外的海山會(huì)對內(nèi)區(qū)R上的VGG異常δΓ有貢獻(xiàn)，因此，式(8)或式(9)需要進(jìn)行調(diào)整。本節(jié)主要以圖4(a)模擬的海山為例來介紹求解方法和進(jìn)行結(jié)果分析。

以后稱40 km×40 km區(qū)域D為內(nèi)區(qū)R的擴(kuò)展區(qū)域，選取步長t=2 km對D進(jìn)行格網(wǎng)化(圖4(b))。稱D-R為內(nèi)區(qū)R的邊界區(qū)域，其下方海山對內(nèi)區(qū)R上VGG異常的貢獻(xiàn)稱為R的邊界效應(yīng)。假設(shè)在內(nèi)區(qū)R上給出了VGG異常值(包含了邊界效應(yīng))，本文的做法是：使用內(nèi)區(qū)R上給出的VGG異常值來求解出整個(gè)擴(kuò)展海域D下方的海底地形，然后剔除D-R下方的海山，僅取R下方海山作為最終的反演結(jié)果。注意：采用的已知數(shù)據(jù)僅是內(nèi)區(qū)R上給出了VGG異常值。

圖4 模擬40 km×40 km海底地形及對應(yīng)格網(wǎng)化Fig.4 The simulated 40 km×40 km seafloor topography and corresponding grid diagram

事實(shí)上，在式(7)中將區(qū)域取為D，而將P(x,y)點(diǎn)僅限于R中，則有

δΓ(x,y)=G·Δρ[ID(H,x,y)-

(10)

(11)

(ATA+αE)h=ATb

(12)

這里E是單位矩陣，α>0是正則化參數(shù)(其單位為10-24m-2s-4)。

因?yàn)棣?0，所以形如式(12)的方程關(guān)于h總是存在唯一解。因此通過迭代求解式(12)的方程組后，便能從式(11)解得所有的hij(i,j=-20,-19，…，18,19)。從中選擇出hij(i,j=-10,-9，…，8,9)，便可得R的子分割Rij的平均海深hij。關(guān)于正則化參數(shù)α的選擇，本文以得到的Rij下方平均海深hij的RMS誤差作為標(biāo)準(zhǔn)，繪制了相應(yīng)的曲線(圖5)。

圖5 不同正則化參數(shù)α對應(yīng)的反演結(jié)果的RMS誤差Fig.5 The RMS error of the results corresponding to different regularization coefficients α

由圖5可知，當(dāng)正則化參數(shù)α=10-3時(shí)海深hij的RMS誤差達(dá)到最小，此時(shí)在內(nèi)區(qū)R上恢復(fù)的海深的誤差分布由圖6給出。由圖6可知，最大誤差為0.64 m，RMS誤差為0.48 m。因此通過將內(nèi)區(qū)R擴(kuò)充到D，并引入正則化參數(shù)，就可以有效地處理D-R下方海山引起的邊界效應(yīng)。

圖6 考慮邊界效應(yīng)后，反演的內(nèi)區(qū)下方海底地形的誤差分布Fig.6 The error distribution of the inverted seabed beneath R after introducing boundary effect

3 遠(yuǎn)區(qū)海山的影響

如果選定了內(nèi)區(qū)R，通過適當(dāng)擴(kuò)充R得到擴(kuò)展區(qū)域D，就可以求解R下方海底地形的邊界效應(yīng)。事實(shí)上，在擴(kuò)展區(qū)域D外部的下方仍有海山會(huì)對R上的VGG值產(chǎn)生影響，本文稱這樣的影響為遠(yuǎn)區(qū)影響。盡管遠(yuǎn)區(qū)影響總體來看在R上是低頻的，但其對R下方海底地形的反演到底影響到何種程度是需要探討的。本節(jié)將在上節(jié)基礎(chǔ)上，通過引入遠(yuǎn)區(qū)海山來進(jìn)一步地分析遠(yuǎn)區(qū)影響及處理方法。

圖7 沿y=0，從D的左邊起算的Ω生成的VGG異常的圖像Fig.7 Along y=0, the graph of VGG anomalies generated by Ω from the left of D

如果將遠(yuǎn)區(qū)影響看成是R上VGG異常的系統(tǒng)誤差，則可直接使用本文處理邊界效應(yīng)的方法來反演R下方海底地形。增加遠(yuǎn)區(qū)影響Ω后，其反演結(jié)果的誤差分布由圖8(a)給出。從圖8(a)可知，受Ω影響，反演的內(nèi)區(qū)R下方的海底地形的最大與RMS誤差分別是71.8 m與39.4 m。

在處理實(shí)際海域問題時(shí)，一旦選定了內(nèi)區(qū)，其遠(yuǎn)區(qū)影響是非常大的，這是因?yàn)閷?shí)際的遠(yuǎn)區(qū)影響是由多種因素產(chǎn)生的。例如：除了本文描述的擴(kuò)展區(qū)域外海山的影響外，還存在著地球深部質(zhì)量分布不均勻等因素的影響。如果相對于內(nèi)區(qū)R來說，遠(yuǎn)區(qū)影響太大，則觀測方程組(11)中會(huì)出現(xiàn)很大的誤差項(xiàng)，這將會(huì)導(dǎo)致反演結(jié)果完全不正確。因此，必須選擇某種途徑去處理遠(yuǎn)區(qū)影響。

盡管遠(yuǎn)區(qū)影響的量級會(huì)很大，但是在R上VGG異常會(huì)表現(xiàn)出低頻的特征，其理由是遠(yuǎn)區(qū)影響距內(nèi)區(qū)R相對較遠(yuǎn)，而VGG自身在離開異常體后會(huì)表現(xiàn)出衰減很快。因此，本文將引入一個(gè)常數(shù)ε來表示R的遠(yuǎn)區(qū)影響，將式(11)改寫為

(13)

然后在式(13)中同時(shí)對hij和ε進(jìn)行求解。至于如何求解式(13)，仍然采用對hij進(jìn)行線性化，再使用正則化參數(shù)的方法來迭代解算hij和ε。

按照上面的方法，對Ω的遠(yuǎn)區(qū)影響問題進(jìn)行了計(jì)算，反演得到了內(nèi)區(qū)R下方的海底地形，其誤差分布由圖8(b)給出，其中最大與均方根誤差分別是58.7 m與32.5 m。相較于圖8(a)的結(jié)果，最大誤差和均方根誤差分別降低18.2%和17.5%。這說明加入遠(yuǎn)區(qū)影響項(xiàng)ε后，確實(shí)是可以消除部分遠(yuǎn)區(qū)影響的。

圖8 考慮遠(yuǎn)區(qū)影響Ω后的反演結(jié)果誤差分布Fig.8 The error distribution of the inverted seabed topography with the disturbance seamount Ω

4 實(shí)際算例及分析

前面各節(jié)分別討論了單個(gè)海山、周邊海山(邊界效應(yīng))、遠(yuǎn)區(qū)海山(遠(yuǎn)區(qū)影響)等情況下如何反演某研究海域R下方海底地形的反演方法。本節(jié)將選取某個(gè)實(shí)際海域作為研究對象，來反演其下方的海底地形。這里選取南中國海72 km×72 km的海域(112.6°E～113.25°E,11.0°N～11.65°N)作為研究對象，密度差Δρ=ρ-ρ0選取為1670 kg/m3。此外，根據(jù)已有的海深數(shù)據(jù)，選擇H=4.4 km。

首先在72 km×72 km的海域D中選擇內(nèi)部52 km×52 km海域R作為內(nèi)區(qū)，其上的VGG異常數(shù)據(jù)(圖9(a))來自GFZ。然后選擇步長t=2 km，由此便可以建立方程式(13)，對hij進(jìn)行線性化，引入正則化參數(shù)α=10-3，迭代后便能解出R的子分割Rij的平均海深hij。最后將海深hij置于Rij的中心，便可生成R下方的海底地形如圖10所示。

圖9 VGG異常分布及船測數(shù)據(jù)點(diǎn)分布Fig.9 The distribution of VGG anomalies and the distribution of bathymetric points

圖10 利用VGG異常反演海底地形的最終結(jié)果Fig.10 The final result of seafloor topography inversion using VGG anomalies

為了檢驗(yàn)反演出的R下方的海底地形的精度，本文選取了該海域全部289個(gè)來自NGDC(National Geophysical Data Center)船測點(diǎn)的海深數(shù)據(jù)，其分布由圖9(b)給出。將本文反演結(jié)果與NGDC數(shù)據(jù)進(jìn)行比較后可知：其最大、最小和均方根誤差分別為339、-157和109 m。統(tǒng)計(jì)反演結(jié)果的平均深度為4173 m，均方根誤差相當(dāng)于平均深度的2.6%(109 m/4173 m)。另外，289個(gè)點(diǎn)上的誤差絕對值統(tǒng)計(jì)見表1，結(jié)果表明：誤差主要集中在150 m以下，占比為84.1%，誤差大于150 m的點(diǎn)，占比為15.9%，且這部分點(diǎn)主要分布于R的邊緣處。如若將最終的反演結(jié)果再向內(nèi)收縮一個(gè)格網(wǎng)，即2 km，反演結(jié)果的最大誤差和均方根誤差將分別下降為207 m和99 m，同比分別下降38.9%和9.2%。

表1 反演結(jié)果誤差統(tǒng)計(jì)Tab.1 Inversion results errors

5 總結(jié)與討論

目前，研究重力場與海底地形的關(guān)系主要還是使用經(jīng)驗(yàn)公式，即，通過部分實(shí)測的海深數(shù)據(jù)來擬合出海底地形與測高、重力和重力梯度等數(shù)據(jù)之間的線性關(guān)系，以此來進(jìn)一步地推算海域的海底地形。顯然，這樣的算法是隨著海域的不同而變化的。本文工作則是依靠建立起的海底地形與VGG異常之間的嚴(yán)格函數(shù)關(guān)系(或觀測方程組)而展開的，理論上看，本文工作是解析的方法。本文：①通過在觀測方程組中加入隨機(jī)誤差和系統(tǒng)誤差后，進(jìn)行模擬計(jì)算，結(jié)果表明觀測方程組具有極強(qiáng)的抗誤差干擾特征；②為了處理邊界效應(yīng)，引入了擴(kuò)充區(qū)域，并在求解線性化方程組時(shí)，引入正則化參數(shù)，確保了線性化方程組的可解性，同時(shí)經(jīng)模擬計(jì)算表明由此解算的海底地形具有很高的精度，統(tǒng)計(jì)出的均方根誤差為0.48 m；③分析了遠(yuǎn)區(qū)影響對反演結(jié)果的影響，使用某個(gè)常數(shù)來吸收遠(yuǎn)區(qū)影響，這樣就能確保遠(yuǎn)區(qū)影響較大時(shí)，觀測方程組依然有很好的可解性；④將本文的理論方法用于實(shí)際海域中，并得到統(tǒng)計(jì)的均方根誤差為109 m的海底地形反演結(jié)果，從而驗(yàn)證了本文方法的有效性。

本文對海域格網(wǎng)化時(shí)，采用的步長為2 km，即在2 km×2 km范圍內(nèi)，僅使用平均海深來表示海山，這意味著在2 km范圍內(nèi)海底地形的起伏是難以分辨出來的。如果將步長取得更小，理論上講對海山的解算精度也會(huì)增加，但是重力場數(shù)據(jù)的分辨率會(huì)制約步長的選擇。例如，EGM2008重力場模型的階次是2160[26]，其空間分辨率約為8 km，如果步長小于4 km，那么來自EGM2008模型的重力數(shù)據(jù)是不能反映出4 km×4 km內(nèi)的地形起伏的。因此，處理實(shí)際問題時(shí)，步長的選擇與重力數(shù)據(jù)的分辨率是相關(guān)的。基于此，若要提高海底地形的反演精度，不僅要提高重力數(shù)據(jù)的精度，而且還需要提高數(shù)據(jù)的分辨率。

本文主要的研究內(nèi)容還是側(cè)重于理論與方法上，即針對建立的從垂直重力梯度數(shù)據(jù)反演海底深度的觀測方程組，進(jìn)行抗誤差性分析以及處理邊界效應(yīng)等分析、驗(yàn)證。因此，在對實(shí)際海底地形進(jìn)行預(yù)測時(shí)，都將海山和海水密度取成了常用的推薦值。事實(shí)上，在實(shí)際預(yù)測海底地形時(shí)，應(yīng)該盡可能地參考當(dāng)?shù)氐牡刭|(zhì)與海洋資料來選擇海山與海水的密度分布。

致謝：感謝GFZ和NGDC提供的計(jì)算數(shù)據(jù)。