一種用于重力場分離的約束插值切割方法

孟慶奎 舒 晴 高 維* 周堅(jiān)鑫 張文志 李 勇

(①自然資源部航空地球物理與遙感地質(zhì)重點(diǎn)實(shí)驗(yàn)室，北京 100083； ②中國自然資源航空物探遙感中心，北京100083； ③中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100083； ④中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所，河北廊坊 065000； ⑤自然資源部地球物理電磁法探測技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北廊坊 065000)

0 引言

重力場分離是重力數(shù)據(jù)反演、地質(zhì)解釋與推斷的前提[1-3]，是業(yè)內(nèi)公認(rèn)的重點(diǎn)研究課題和技術(shù)難題。早在1952年，Dobrin等[4]就提出重力解釋中兩個(gè)最重要的問題：一個(gè)是區(qū)域場與局部場的分離； 另一個(gè)是基于分離得到的異常獲取目標(biāo)地質(zhì)體的空間位置和密度信息。程方道等[5]指出，與磁異常相比，重力異常普遍存在強(qiáng)大的區(qū)域背景場，往往會(huì)造成局部異常發(fā)生畸變，甚至難于辨識(shí)，要對(duì)重力資料進(jìn)行有效的地質(zhì)解釋，首先必須分離重力區(qū)域場與局部場； 管志寧等[6]強(qiáng)調(diào)，正確劃分不同規(guī)模的區(qū)域場與局部場是面積性測量資料解釋能否獲得良好地質(zhì)效果的一個(gè)重要因素； 曾華霖[7]認(rèn)為，在密度界面反演中，只要平均深度取值合理，界面之間的密度差異比較接近實(shí)際情況，無論采用何種反演方法，都能得到達(dá)到一定的精度及相似的結(jié)果，但前提是用于反演的異常必須單純由目標(biāo)體引起； 徐世浙等[8]認(rèn)為重力場三維反演是個(gè)難度極大的課題，提出異于常規(guī)解決方案的開拓性思路，主要包括三個(gè)方面：對(duì)平面上的重力場進(jìn)行不同深度的分離； 用大深度位場向下延拓的迭代法將各個(gè)切割層在地面的重力場下延至相應(yīng)深度； 將頂面的重力場反演為視密度?，F(xiàn)今，位場大深度的向下延拓技術(shù)和層源異常的反演技術(shù)已經(jīng)得到了很好發(fā)展，位場分離方法是目前影響反演方法效果的關(guān)鍵技術(shù)，繼續(xù)加強(qiáng)對(duì)異常分離方法的研究具有重要意義。

重力場分離方法大致可分為空間域方法和頻率域(波數(shù)域)方法。頻率域方法主要包括維納濾波法[9]、匹配濾波法[10-12]、小波多尺度分析法[13-17]、優(yōu)選濾波法[18-19]和優(yōu)化濾波法[20]等，這些方法各有特點(diǎn)和優(yōu)勢，但同時(shí)又有不同的局限性。維納濾波法是維納濾波器在信號(hào)與干擾不相關(guān)條件下的一個(gè)特例，適用于白噪聲的去除，不適合于場的分離，目前相關(guān)研究與應(yīng)用不多； 匹配濾波法對(duì)于埋深和尺寸相差較大的場源分離效果較好，但是若異常的功率譜不具明顯的線性特征時(shí)，如何更加有效地確定深、淺源場的頻段分布，是該方法在實(shí)際資料處理過程中的一個(gè)難點(diǎn)； 小波多尺度分析法是一種新興的位場分離方法，在區(qū)域性深部地殼結(jié)構(gòu)研究方面取得了一些成果，但是存在小波母函數(shù)選擇和階次選擇較難等問題； 優(yōu)選濾波法實(shí)現(xiàn)了指定頻段或尺度的局部異常的分離，但實(shí)際應(yīng)用中延拓高度往往很難正確選擇； 優(yōu)化濾波法異常分離不受延拓高度的困擾，但其前提是目標(biāo)層與其他深度層場源信息互不相關(guān)。實(shí)踐表明，空間域方法有一定的技術(shù)優(yōu)勢。管志寧等[6]基于多年的重磁數(shù)據(jù)處理經(jīng)驗(yàn)認(rèn)為，由于受測區(qū)范圍有限、數(shù)據(jù)離散化，以及不同埋深、不同位置場源體的重磁異常頻譜特征部分重疊等諸多因素的影響，在波數(shù)域變換過程中會(huì)產(chǎn)生明顯的誤差和畸變，只能用作初步的定性估計(jì)； Spector[21]也認(rèn)為，波數(shù)域方法只能對(duì)深源(區(qū)域場)和淺源(局部場)異常作定性的劃分，不能用于定量研究，無法取代空間域中異常曲線擬合類方法。因而，本文以空間域方法為研究對(duì)象。

空間域方法主要包括趨勢面擬合法[22-24]、經(jīng)驗(yàn)?zāi)B(tài)分解法[25]、最小曲率法[26-28]、曲率濾波法[29-30]、有限單元法[31-37]、卡爾曼濾波法[38-39]、矩陣底秩分解法[40-41]和插值切割法[42-53]等。目前，空間域方法中較為前沿且應(yīng)用較廣的重力場分離方法是插值切割法，該方法基于多次切割法發(fā)展而來。程方道等[5]首次提出多次切割法，其基本原理是以四點(diǎn)圓周平均值為切割算子，通過連續(xù)切割得到切割區(qū)域場，并將其從總異常中減去，便得到切割局部場。為了壓制多次切割法中的切割發(fā)散效應(yīng)，文百紅等[42]基于新型算子首次提出了插值切割法； 之后，秦葆瑚[43]、段本春等[44]、徐世浙等[45]、邢怡等[46]、劉東甲[47]、葛粲等[48-49]圍繞優(yōu)化切割算子、探尋切割半徑與異常埋深的關(guān)系展開了系統(tǒng)的研究工作，形成了以插值切割法為核心技術(shù)的方法系列，汪炳柱等[50]、楊金玉等[51]、萬勝榮等[52]和趙文舉等[53]將該方法廣泛應(yīng)用于油氣勘探、礦產(chǎn)勘查、基礎(chǔ)地學(xué)研究等諸多領(lǐng)域。

以上提及的插值切割法及其改進(jìn)算法很少涉及重力異常產(chǎn)生的物理屬性，多針對(duì)異常曲線本身進(jìn)行插值切割計(jì)算。如果忽視異常產(chǎn)生的地質(zhì)屬性，往往會(huì)造成分析問題不夠全面，甚至產(chǎn)生一定的計(jì)算誤差。故本文在分析插值切割法基本原理基礎(chǔ)上，證明了插值切割法的切割算子僅考慮了局部異常與區(qū)域異常彎曲方向相同的疊合情況，不能完全體現(xiàn)實(shí)際復(fù)雜的地下地質(zhì)情況。因此，本文結(jié)合重力局部異常與區(qū)域異常的四種疊合關(guān)系，提出了一種基于異常約束條件的插值切割法，并將其命名為約束插值切割法。通過典型理論模型對(duì)比試算和實(shí)際測量數(shù)據(jù)應(yīng)用，驗(yàn)證了約束插值切割法的準(zhǔn)確性和優(yōu)越性，可為大區(qū)域的基礎(chǔ)地學(xué)研究提供技術(shù)支撐，在具備局部高精度重力資料的情況下，也可為油氣遠(yuǎn)景區(qū)的預(yù)測提供一定技術(shù)參考。

1 插值切割法

插值切割法是一種簡單、有效的空間域非線性濾波方法，以當(dāng)前計(jì)算點(diǎn)場值與周圍多點(diǎn)平均值的插值運(yùn)算為切割算子，通過連續(xù)切割得到切割區(qū)域場，并將其從位場異常中減去，便得到切割局部場[42]。根據(jù)插值切割算子的不同，形成了基于四點(diǎn)圓周插值切割算子[42]、八點(diǎn)窗口插值切割算子[46]和動(dòng)態(tài)改進(jìn)型插值切割算子[49]的插值切割方法。

1.1 計(jì)算思路

插值切割法雖衍生出不同的切割算子，但其計(jì)算原理是一致的。下面以經(jīng)典的四點(diǎn)圓周插值切割算子為例介紹插值切割法的計(jì)算思路。

設(shè)重力異常Z(x,y)由區(qū)域場R(x,y)和局部場L(x,y)組成，即

Z(x,y)=R(x,y)+L(x,y)

(1)

令M(x,y)表示與點(diǎn)(x,y)距離為r的4個(gè)點(diǎn)的重力異常平均值，即

Z(x,y+r)+Z(x,y-r)]

(2)

式中r也稱為切割半徑。區(qū)域場R(x,y)是Z(x,y)與M(x,y)的加權(quán)平均，即

R(x,y)=aM(x,y)+bZ(x,y)

(3)

式中a和b是加權(quán)系數(shù)，且a+b=1。令

(4)

(5)

ΔHx=Z(x+r,y)+Z(x-r,y)

(6)

ΔHy=Z(x,y+r)+Z(x,y-r)

(7)

(8)

(9)

e=c+d(0≤e≤2)

(10)

(11)

用上述方法得到的區(qū)域場稱為第1次切割的區(qū)域場，用R1(x,y)表示； 對(duì)R1(x,y)重復(fù)使用以上方法，得到第2次切割的區(qū)域場R2(x,y)； 依次迭代，當(dāng)

(12)

有

(13)

上式得到的區(qū)域場即為切割半徑的區(qū)域場，將其與重力異常Z(x,y)相減，得到局部場

L(x,y)=Z(x,y)-R(x,y)

(14)

一般地，切割半徑r越大，得到的局部場所反映的密度體的深度和規(guī)模就越大。根據(jù)前人大量模擬實(shí)驗(yàn)和實(shí)測資料[42-53]分析，以切割半徑r得到的局部場可近似代表深度r以上地層的重力效應(yīng)，但尚未得出切割半徑r與密度體埋深的確切計(jì)算公式。由于切割算子中的加權(quán)系數(shù)引入了半二階差分量，它們與負(fù)二階水平導(dǎo)數(shù)即曲率成正比，因此插值切割對(duì)不同非線性異常有不同的作用，可提高對(duì)不同特征異常場的分離效果。

1.2 誤差傳遞分析

插值切割法采用當(dāng)前計(jì)算點(diǎn)值與四點(diǎn)圓周平均值之間的插值關(guān)系作為切割算子G，將切割算子G作用于重力場Z(x,y)，即可得到區(qū)域場

(15)

設(shè)原始觀測數(shù)據(jù)各點(diǎn)的均方誤差δZ(x,y)=ε0，根據(jù)誤差傳遞公式[54]得到區(qū)域場的誤差傳遞關(guān)系為

[δM(x,y)]2

(16)

那么，第1次切割后區(qū)域場的誤差傳遞系數(shù)k1有

(17)

第i次切割后區(qū)域場的誤差傳遞系數(shù)ki有

(18)

由上式可見，多次切割后區(qū)域場的傳遞誤差不會(huì)放大。

2 改進(jìn)的插值切割法

分析插值切割法基本理論公式，該方法主要對(duì)異常曲線本身進(jìn)行插值切割計(jì)算，而對(duì)重力異常產(chǎn)生的地質(zhì)屬性考慮不足。忽略產(chǎn)生異常的地質(zhì)屬性得到的計(jì)算結(jié)果難免產(chǎn)生一定的誤差。下面首先在理論層面分析插值切割法的不足之處，再從實(shí)際地質(zhì)屬性出發(fā)，結(jié)合重力場局部異常與區(qū)域異常的四種疊合關(guān)系，對(duì)插值切割法加以改進(jìn)，提出一種基于異常約束條件的插值切割法，即約束插值切割法。

2.1 插值切割法不足的理論分析

根據(jù)插值切割法的基本理論可知，切割區(qū)域場計(jì)算公式(式(15))中，0≤e≤2，不難證明R(x,y)介于Z(x,y)與M(x,y)之間，這就意味著經(jīng)插值切割法計(jì)算得到的區(qū)域場值不能超出總場異常值與周圍點(diǎn)平均值的取值區(qū)間。這種情形僅代表局部異常與區(qū)域異常彎曲方向相同的疊合情況，并不能完全代表復(fù)雜多變的實(shí)際地下地質(zhì)情況，下面進(jìn)一步分析局部異常與區(qū)域異常不同疊合關(guān)系的情形。

2.2 局部異常與區(qū)域異常疊合關(guān)系分類

在勘探地球物理領(lǐng)域，一般認(rèn)為區(qū)域異常是由埋藏較深、水平延伸范圍較大的地質(zhì)體或構(gòu)造所產(chǎn)生的，在異常圖中表現(xiàn)出形態(tài)寬緩、幅值較大、范圍較廣的特征； 局部異常則是由埋藏較淺、水平延伸范圍較小的地質(zhì)體所產(chǎn)生的，在異常圖中表現(xiàn)出形態(tài)陡窄、幅值較小、異常范圍小的特征。在遵循這一基本原則的前提下，本文針對(duì)四種實(shí)際地質(zhì)情形，總結(jié)了相應(yīng)的局部異常與區(qū)域異常的四種疊合關(guān)系分類情形，詳見表1和圖1。不難發(fā)現(xiàn)，式(11)僅適用于分類情形1和2，即區(qū)域異常介于總場異常值與周圍點(diǎn)平均值的取值區(qū)間。然而，對(duì)于區(qū)域異常與局部異常曲線彎曲方向相反的情形，即分類情形3和4，如仍使用式(11)計(jì)算切割區(qū)域場，計(jì)算結(jié)果雖可無限向周圍點(diǎn)平均值逼近，但依然忽略了圖1c和圖1d中黑色豎線所示部分的異常，這部分異常代表計(jì)算得到切割區(qū)域場與理論區(qū)域場的差異。

表1 局部異常與區(qū)域異常疊合關(guān)系分類表

圖1 對(duì)應(yīng)表1中分類情形1～4(a～d)的局部異常與區(qū)域異常疊合關(guān)系分類

2.3 算法改進(jìn)方案

為彌補(bǔ)傳統(tǒng)插值切割法的不足，本文提出了一種基于異常約束條件的改進(jìn)方案，其基本思路是：首先，給定切割半徑r，分析區(qū)域異常和局部異常曲線彎曲方向； 其次，對(duì)區(qū)域異常與局部異常曲線彎曲方向相同的情形(即分類情形1和2)，采用式(11)計(jì)算切割區(qū)域場，對(duì)彎曲方向相反的情形(即分類情形3和4)，在式(11)基礎(chǔ)上，增加切割區(qū)域異常修正量ΔZ(x,y)(圖2)，即

圖2 切割區(qū)域異常修正量示意圖

(19)

首先考慮x方向異常剖面，假設(shè)局部異常附近的區(qū)域異常是曲率半徑為N、曲率中心為O的曲率圓的部分弧段(圖2中紅色曲線與藍(lán)色曲線重合部分)，則x方向切割區(qū)域的異常修正量為

(20)

(21)

(22)

式中Z′D，x、Z″D，x分別為Z(x,y)在D點(diǎn)對(duì)x的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)。用一階差分替代一階導(dǎo)數(shù)，并忽略一階余量，得到

(23)

(24)

同理，可得到y(tǒng)方向切割區(qū)域的異常修正量ΔZ(y)，并取ΔZ(x,y)=[ΔZ(x)+ΔZ(y)]/2。

約束插值切割法詳細(xì)計(jì)算步驟見附錄A。以上算法基于MATLAB(R2008a)編程平臺(tái)實(shí)現(xiàn)程序和界面編制。

3 理論模型試算

3.1 理論模型設(shè)計(jì)

本文理論模型設(shè)計(jì)參考了程方道等[5]、Kaftan等[31]、Zhu等[40]及劉東甲[47]的設(shè)計(jì)思路，即深部大尺度密度體產(chǎn)生的異常視為區(qū)域異常，淺部小尺度密度體產(chǎn)生的異常視為局部異常，不同之處是在同一個(gè)理論模型中綜合考慮了局部異常與區(qū)域異常的四種疊合關(guān)系。鑒于2．2部分中分類情形1、3分別與2、4呈鏡像對(duì)稱關(guān)系，可合并處理，故本文設(shè)計(jì)了圖3所示的理論模型。三個(gè)不同深度的淺部球體(編號(hào)①、②、③)產(chǎn)生的是局部重力異常，其中淺部球體①、②產(chǎn)生局部重力高，淺部球體③產(chǎn)生局部重力低； 三個(gè)不同深度的深部球體(編號(hào)④、⑤、⑥)產(chǎn)生的是區(qū)域重力異常，其中球體④、⑥引起的是區(qū)域重力高異常，球體⑤引起的是重力低異常。

圖3 理論模型

3.2 理論模型正演

理論模型的重力異常為各球體剩余質(zhì)量的引力位沿重力方向(定義為z方向)的導(dǎo)數(shù)之和[55]，即

(25)

基于式(25)，通過正演計(jì)算可得模型3的正演重力異常(圖4、圖5)。由圖4可見，局部異常基本淹沒在區(qū)域異常之中。由圖5可見，1、2、3號(hào)異常分別對(duì)應(yīng)2．2部分中分類情形1、4、3。

圖4 圖3模型重力異常正演結(jié)果

圖5 圖3模型重力異常正演剖面(y=0)

3.3 理論模型重力場分離

分別采用插值切割法、約束插值切割法、小波多尺度分解法[13]對(duì)理論模型正演結(jié)果進(jìn)行重力場分離，分離得到的區(qū)域場和局部異常見圖6?；谇懈畎霃脚c密度體埋深近似相等的對(duì)應(yīng)關(guān)系，計(jì)算過程中以r=400m為初始值，以100m為增減量，采用多個(gè)切割半徑r，經(jīng)對(duì)比篩選，確定取r=500m。根據(jù)小波多尺度分解的尺度場源深度轉(zhuǎn)換公式[14]，因網(wǎng)格間距為100m，故采用1～4階小波細(xì)節(jié)代表淺部球體產(chǎn)生的局部重力異常，用4階小波逼近代表深部球體產(chǎn)生的區(qū)域重力異常。對(duì)比約束插值切割法和小波多尺度分解法得到的局部異常，可見前者更接近于理論局部異常(圖4c)。由于區(qū)域異常幅值較大，基本淹沒了局部異常，故通過三種方法分離出的區(qū)域異常(圖6上)與理論區(qū)域異常(圖4b)基本一致。為了更直觀地分析本文方法的優(yōu)越性，繪制了沿y=0剖面的三種方法分離得到的局部異常曲線(圖7)。由圖可見，由于1號(hào)異常為區(qū)域重力高疊加局部重力高(分類情形1)，三種方法的結(jié)果基本一致，都能較準(zhǔn)確地體現(xiàn)分類情形1所產(chǎn)生的局部重力異常； 關(guān)于2號(hào)和3號(hào)異常，本文方法得到的分離結(jié)果明顯改善了插值切割法的計(jì)算結(jié)果，且與小波多尺度分解的結(jié)果吻合。三種方法均在局部異常兩翼附近產(chǎn)生了小幅度假異常，這是因?yàn)榉椒ū旧泶嬖谝欢ǖ挠?jì)算誤差。以上分析證明了本文提出的約束插值切割法基本理論的正確性和分離效果的正確性。

圖6 圖3模型重力場不同方法分離得到的區(qū)域異常場(上)和局部異常場(下)

圖7 不同方法分離的局部重力異常場(y=0)

為了定量對(duì)比分析插值切割法、約束插值切割法和小波多尺度分解法的重力場分離效果，用下式計(jì)算局部異常ΔgⅠ與理論局部異常ΔgⅡ的相關(guān)系數(shù)

(26)

式中：Cor(·)表示協(xié)方差； Var(·)表示均方差。經(jīng)計(jì)算，插值切割法、約束插值切割法和小波多尺度分解法得到的局部異常與理論局部異常的相關(guān)系數(shù)分別為0.42、0.75、0.70，可見本文提出的約束插值切割法得到的局部異常與理論局部異常具有最強(qiáng)的相關(guān)性。

另外，為了評(píng)價(jià)約束插值切割法對(duì)插值切割法的改善效果，引入評(píng)價(jià)航磁動(dòng)態(tài)補(bǔ)償效果的“改善率”(Improvement Ratio，IR)[56]評(píng)價(jià)改善效果，其計(jì)算公式為

(27)

式中：Δgcz表示插值切割法得到的局部異常； Δgys表示約束插值切割法得到的局部異常； Std(·)表示求標(biāo)準(zhǔn)差。

因?yàn)楸疚奶岢龅男路椒▋H對(duì)2號(hào)、3號(hào)異常進(jìn)行了優(yōu)化，故僅對(duì)2號(hào)、3號(hào)局部異常的改善率進(jìn)行評(píng)估。經(jīng)計(jì)算，約束插值切割法對(duì)插值切割法的改善率IR為1.48，可見相比于插值切割法，新方法計(jì)算的局部異常效果改善了48%，優(yōu)化效果明顯。

4 應(yīng)用實(shí)例

4.1 中國大陸及周邊區(qū)域重力數(shù)據(jù)分離

本文以中國大陸及周邊區(qū)域的重力數(shù)據(jù)為試驗(yàn)對(duì)象，范圍為東經(jīng)63°～136°，北緯15°～58°。重力數(shù)據(jù)來自美國國家地理空間情報(bào)局地球重力場研發(fā)小組發(fā)布的超高階地球重力場模型EGM2008(Earth Gravitational Model 2008)，該模型的階次完全至2159次，空間分辨率約為5′，可用于小比例尺重力編圖和構(gòu)造研究[57]。本實(shí)例使用的布格重力異常的網(wǎng)格間距為4km。

分別利用插值切割法和本文提出的約束插值切割法對(duì)布格重力異常進(jìn)行分離。為了展示新方法的正確性和改進(jìn)效果，以葛粲等[48]的研究結(jié)果為參照。根據(jù)徐世浙等[8]、葛粲等[48]提出的“切割半徑應(yīng)近似等于異常體埋深”的層切割思想，本例的切割半徑r以網(wǎng)格間距4km為基準(zhǔn)，并以其整數(shù)倍的形式逐步遞增，每次切割得到的區(qū)域場作為下次切割的輸入，并結(jié)合前人經(jīng)驗(yàn)總結(jié)的“切割半徑r得到的局部場可以近似代表深度r以上的地層重力效應(yīng)”，當(dāng)切割半徑為19倍網(wǎng)格間距(r=76km)時(shí)，計(jì)算得到的局部異?？梢暈?2～76km參考深度層的重力異常值，并將其視為上地幔頂部附近介質(zhì)產(chǎn)生的重力異常。插值切割法和約束插值切割法計(jì)算得到的上地幔頂部附近重力異常分別見圖8和圖9。

對(duì)比圖8與葛粲等[48]的計(jì)算結(jié)果，可見重力異常分布特征基本一致，驗(yàn)證了本文編程實(shí)現(xiàn)的插值切割法的正確性。對(duì)比插值切割法(圖8)與約束插值切割法(圖9)得到的上地幔頂部附近重力異常，可見整體上重力異常形態(tài)相似，但約束插值切割法得到的正值重力異常幅度和范圍相對(duì)更大，說明在負(fù)值背景下疊加正異常的情況下，約束插值切割法對(duì)異常分布進(jìn)行了合理調(diào)整，這與約束插值切割法的基本原理和理論模型結(jié)論相一致。

圖8 插值切割法得到的中國大陸及周邊地區(qū)上地幔頂部附近重力異常場

圖9 約束插值切割法得到的中國大陸及周邊地區(qū)上地幔頂部附近重力異常場

為了定量對(duì)比約束插值切割法和插值切割法得到的重力異常，計(jì)算了兩者的差值(圖10)，即前文所述的修正量。由圖可見，在準(zhǔn)噶爾盆地、喜馬拉雅地塊等局部地區(qū)，本文方法的修正達(dá)50%以上。

圖10 圖9與圖8的差值

根據(jù)圖9可得到初步的地質(zhì)認(rèn)識(shí)，即在喜馬拉雅地塊邊界局部地區(qū)(圖中綠色虛線橢圓區(qū)域)，深部高密度支撐物質(zhì)強(qiáng)烈北移侵入，而柴達(dá)木盆地和準(zhǔn)噶爾盆地正高值范圍的擴(kuò)大(圖中白色虛線所示)，則意味著穩(wěn)定盆地下方可能存在高密度物質(zhì)。

楊文采等[58]利用小波多尺度分解法對(duì)中國陸域及周邊布格重力異常進(jìn)行了分離，獲取了包括上地幔頂部異常(D7+D8)在內(nèi)的4個(gè)深度界面所產(chǎn)生的重力異常。對(duì)比約束插值切割法與小波多尺度分解法對(duì)上地幔頂部異常的分離結(jié)果，可見在青藏高原等區(qū)域的異常分布特征與楊文采等[58]的研究成果相近，但約束插值切割法得到的異常細(xì)節(jié)更加豐富，主要是因?yàn)樾〔ǘ喑叨确纸夥▋H將布格重力異常分離為4層，而約束插值切割法可分離為20層，對(duì)異常的分辨率更高[48]。

由于篇幅限制，結(jié)合地震、地質(zhì)等資料的綜合地質(zhì)解釋將另文詳細(xì)闡述。據(jù)本例分析可見，約束插值切割法在實(shí)際應(yīng)用中實(shí)現(xiàn)了對(duì)插值切割法計(jì)算結(jié)果的修正，可得到更準(zhǔn)確的重力異常、更明確的重力異常梯度帶。

4.2 中國南海萬安盆地及周邊區(qū)域重力數(shù)據(jù)分離

為了進(jìn)一步驗(yàn)證約束插值切割法的實(shí)際應(yīng)用效果，選取中國南海萬安盆地及周邊區(qū)域布格重力異常數(shù)據(jù)為研究對(duì)象，范圍為東經(jīng)104°～111°，北緯5°～11°，數(shù)據(jù)來源同4.1部分重力數(shù)據(jù)。陳玲等[59]認(rèn)為南海萬安盆地新生界沉積層厚度最大為12.5km，馮旭亮等[60]的研究結(jié)果顯示該區(qū)域最大厚度大于10km，故本例約束插值切割法中選取切割半徑為3倍網(wǎng)格間距(12km)，經(jīng)計(jì)算得到萬安盆地及周邊區(qū)域新生界沉積層的重力異常(圖11a)。圖11b為選取的一條地震剖面所對(duì)應(yīng)的新生界沉積層重力異常，圖11c為該地震剖面構(gòu)造解釋剖面[60]。

為了說明本文方法計(jì)算結(jié)果的可靠性，與馮旭亮等[60]采用最小曲率位場分離方法得到的計(jì)算結(jié)果(文獻(xiàn)[60]中的圖7、圖8)進(jìn)行了對(duì)比，可以發(fā)現(xiàn)，圖11a中重力異常的整體分布特征與文獻(xiàn)[60]中的圖7較吻合，圖11b中的重力異常曲線走勢與文獻(xiàn)[60]中的圖8a大致相同，但本文的計(jì)算結(jié)果更能反映沉積基底的起伏變化細(xì)節(jié)，特別是更加直觀地顯示了F1斷層右側(cè)沉積基底的局部凸起(圖11c中的虛線圈所示)。因此，基于本文方法分離得到的局部重力異常，可更清晰地解釋新生界沉積層底界面的深度變化，為該區(qū)域新生界深度反演研究提供了較可靠的數(shù)據(jù)基礎(chǔ)，同時(shí)也驗(yàn)證了約束插值切割法的實(shí)用性。

圖11 中國南海萬安盆地及周邊地區(qū)重力數(shù)據(jù)及地震剖面

5 結(jié)論

本文通過基本理論推導(dǎo)，并結(jié)合實(shí)際重力異常的疊合分類情形，實(shí)現(xiàn)了對(duì)插值切割法的優(yōu)化和改進(jìn)，提出了基于異常約束條件的插值切割法(約束插值切割法)，并開展了理論模型試驗(yàn)和衛(wèi)星重力數(shù)據(jù)處理、分析。得出以下幾點(diǎn)認(rèn)識(shí)。

(1)通過對(duì)插值切割法切割算子的公式推導(dǎo)，結(jié)合四種情形的異常疊合關(guān)系，證明了插值切割法的切割算子僅考慮了局部異常與區(qū)域異常彎曲方向相同的疊合情況?；诖?，給出了插值切割算子的補(bǔ)償方法，提出了基于異常約束的插值切割法，制定了新方法的計(jì)算流程，并基于Matlab平臺(tái)編程實(shí)現(xiàn)。

(2)理論模型試算結(jié)果顯示，插值切割法、約束插值切割法、小波多尺度分解法中約束插值切割法得到的局部異常與理論局部異常具有最強(qiáng)的相關(guān)性(相關(guān)系數(shù)為0.75)，與插值切割法相比，對(duì)局部異常效果的改善達(dá)48%，優(yōu)化效果明顯。

(3)對(duì)中國陸域及周邊區(qū)域的布格重力異常進(jìn)行了分離，獲取了72～76km參考深度層產(chǎn)生的重力異常，約束插值切割法獲取了該層位更真實(shí)可靠的重力異常值信息，在準(zhǔn)噶爾盆地、喜馬拉雅地塊等局部地區(qū)的修正達(dá)50%以上，可為上地幔頂部構(gòu)造特征研究提供參考。

(4)利用約束插值切割法分離得到了萬安盆地及周邊區(qū)域新生界沉積層重力異常，與地震剖面顯示的基底起伏特征基本一致，說明該方法可用于基底深度變化的反演研究。

總之，通過理論模型數(shù)據(jù)和實(shí)際重力數(shù)據(jù)分析，約束插值切割法計(jì)算結(jié)果均明顯優(yōu)化了插值切割法的計(jì)算結(jié)果，驗(yàn)證了前者理論的正確性及實(shí)用性，可為實(shí)際重力數(shù)據(jù)的地質(zhì)解釋提供更準(zhǔn)確、可靠的重力異常信息。本文提出的新方法也可用于磁場數(shù)據(jù)的異常分離。

在構(gòu)思本文的初期，與插值切割法的原創(chuàng)者——中國石油勘探開發(fā)研究院文百紅老師進(jìn)行了有益探討，特別是在插值切割法的細(xì)節(jié)上得到了文百紅老師的悉心指導(dǎo)，在此深表謝意； 感謝中國自然資源航空物探遙感中心駱遙教授為文章的修改提出了建設(shè)性的意見！

附錄A 約束插值切割法計(jì)算步驟

(1)根據(jù)目標(biāo)層位深度信息，對(duì)重力異常Z(x,y)設(shè)置切割半徑r；

(2)利用式(2)計(jì)算與點(diǎn)(x,y)相距r的某4個(gè)點(diǎn)重力異常的平均值M(x,y)；

(3)利用式(10)計(jì)算加權(quán)參數(shù)e;

(4)比較區(qū)域異常與局部異常曲線彎曲方向，具體可細(xì)分為2步：

①計(jì)算下列參數(shù)

②當(dāng)比較參數(shù)不能同時(shí)滿足下列條件時(shí)，表明區(qū)域異常和局部異常彎曲方向相同(即分類情形1和2)

Ix1Iy1>0

Iy1Iy2<0

Ix1Ix2<0

Iy2Ix2>0：

當(dāng)比較參數(shù)同時(shí)滿足下列條件時(shí)，表明區(qū)域異常和局部異常彎曲方向相反(即分類情形3和4)

Ix1Iy1>0

Iy1Iy2<0

Ix1Ix2<0

Iy2Ix2>0：

(5)對(duì)于彎曲方向相同的情形，區(qū)域場等于Z(x,y)與M(x,y)的加權(quán)平均，采用式(11)計(jì)算切割區(qū)域場；

(6)對(duì)于彎曲方向相反的情形，區(qū)域場為Z(x,y)與M(x,y)的加權(quán)平均的基礎(chǔ)上增加切割區(qū)域異常修正量ΔZ(x,y)，此種情況可細(xì)分為以下4步：

①首先考慮x方向異常剖面，假設(shè)局部異常附近的區(qū)域異常為曲率半徑為N、曲率中心為O的曲率圓的部分弧段，利用式(20)計(jì)算x方向切割區(qū)域異常修正量ΔZ(x)。

②同理可得y方向切割區(qū)域異常修正量ΔZ(y)；

③取ΔZ(x,y)=[ΔZ(x)+ΔZ(y)]/2；

④計(jì)算修正量ΔZ(x,y)后，采用式(19)計(jì)算切割區(qū)域場；

(7)用上述方法得到第1次切割的區(qū)域場R1(x,y)，對(duì)R1(x,y)重復(fù)使用步驟(1)～(6)，得到第2次切割的區(qū)域場R2(x,y)，依次迭代，直到滿足式(12)，迭代終止，得到切割區(qū)域場；

(8)利用式(14)計(jì)算最終的局部場R(x，y)。