石甲強(qiáng) 肖 鋒 鐘 煬

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,吉林長(zhǎng)春 130026)

0 引言

重力歸一化總梯度法由前蘇聯(lián)學(xué)者Березкин[1]于1967年提出，因其先對(duì)重力異常進(jìn)行向下解析延拓，再計(jì)算歸一化總梯度而得名。它是一種利用高精度重力異常確定場(chǎng)源、斷裂位置及密度界面的方法。肖一鳴[2]首次將該方法引入中國(guó)，給出了利用泰勒級(jí)數(shù)展開(kāi)計(jì)算重力歸一化總梯度的方法，并探討了影響該方法應(yīng)用效果的因素，如諧波數(shù)、隨機(jī)噪聲、測(cè)線長(zhǎng)度等。隨后該方法在中國(guó)得到了廣泛應(yīng)用。肖一鳴等[3]將該方法成功應(yīng)用于油氣勘探；王家林等[4]利用該方法分析斷層和密度分界面；吳燕岡等[5]提出將重力歸一化總梯度與相位疊置處理，用于確定深大斷裂的空間位置；張鳳旭等[6-7]利用Hilbert變換改進(jìn)重力歸一化總梯度法，并提出在位場(chǎng)轉(zhuǎn)換的向下延拓和求導(dǎo)過(guò)程中分別引入圓滑濾波因子，提高了該方法的分辨率、可靠性和穩(wěn)定性，增大了延拓深度；肖鵬飛等[8]利用向下延拓的迭代算法替代波數(shù)域的直接向下延拓算子，提高了重力歸一化總梯度法的穩(wěn)定性；張鳳琴等[9]采用DCT變換計(jì)算重力歸一化總梯度，增加了下延深度；王彥國(guó)等[10]和郭燦燦等[11]提出了基于泰勒級(jí)數(shù)迭代法的快速穩(wěn)定向下延拓的重力歸一化總梯度法；蘇超等[12]對(duì)歸一化函數(shù)的分母計(jì)算幾何平均，提高了重力歸一化總梯度法的抗噪能力；王選平等[13]將正則化因子引入重力場(chǎng)的下延計(jì)算，提出了利用正則化方法計(jì)算重力歸一化總梯度的方法；王彥國(guó)等[14]提出基于冪次平均的離散歸一化總梯度法，可以有效識(shí)別疊加場(chǎng)源的位置信息。由此可見(jiàn)，對(duì)重力歸一化總梯度法的改進(jìn)主要集中在向下延拓和求導(dǎo)這兩步計(jì)算過(guò)程，其穩(wěn)定性和精度直接影響重力歸一化總梯度法的穩(wěn)定性和精度。

張沖等[15-17]針對(duì)向下延拓問(wèn)題，提出了基于求解微分方程的Adams-Bashforth公式法和基于Simpson求積公式的Milne法。前者下延深度能達(dá)到5倍點(diǎn)距，后者甚至可超過(guò)10倍點(diǎn)距。Milne法利用觀測(cè)面上的重力垂向一階導(dǎo)數(shù)、向上延拓不同高度的重力場(chǎng)和重力垂向一階導(dǎo)數(shù)求解下延深度的重力場(chǎng)。向下延拓Milne法具有下延深度大、相對(duì)誤差小的特點(diǎn)，能獲得穩(wěn)定準(zhǔn)確的結(jié)果。但該方法需要計(jì)算垂向?qū)?shù)，且求導(dǎo)的精度直接影響向下延拓的精度。

早在2001年，F(xiàn)edi等[18]基于重力位的二階導(dǎo)數(shù)在場(chǎng)源外滿足拉普拉斯方程的原理，提出了對(duì)隨機(jī)噪聲有壓制作用的積分垂向二階導(dǎo)數(shù)法(Integrated Second Vertical Derivative,ISVD)。之后，崔瑞華等[19]將該方法引入中國(guó)，首先對(duì)重力異常積分，計(jì)算得到重力位，然后計(jì)算重力位的水平二階導(dǎo)數(shù)，最后根據(jù)拉普拉斯方程計(jì)算重力異常的垂向一階導(dǎo)數(shù)。相對(duì)于常規(guī)的波數(shù)域求垂向?qū)?shù)的方法，ISVD法由于進(jìn)行了積分運(yùn)算，對(duì)隨機(jī)噪聲有抑制作用，因而更穩(wěn)定。

本文將這兩種穩(wěn)定算法，即向下延拓Milne法和ISVD法引入重力歸一化總梯度的計(jì)算，通過(guò)理論模型試算和實(shí)際資料處理，驗(yàn)證了方法的有效性。

1 方法原理

1.1 重力歸一化總梯度法的基本原理

以重力剖面觀測(cè)數(shù)據(jù)為例，歸一化總梯度的表達(dá)式[2]為

(1)

式中：GH(x,z)為計(jì)算點(diǎn)(x,z)處的重力歸一化總梯度值；M為重力剖面上的測(cè)點(diǎn)總數(shù)；Vzx(x,z)和Vzz(x,z)為計(jì)算點(diǎn)的重力位二階偏導(dǎo)數(shù)，同時(shí)也是計(jì)算點(diǎn)(x,z)處重力異常沿x方向和z方向的一階偏導(dǎo)數(shù)。

1.2 向下延拓Milne法的基本原理

重力場(chǎng)向下延拓四階Milne公式[15]為

Vzz(x,z-h)+2Vzz(x,z0)]

(2)

式中：h為向下延拓深度；Vz(x,zh)表示由觀測(cè)面z0向下延拓h深度的重力場(chǎng)；Vz(x,z-3h)表示由觀測(cè)面向上延拓3h高度的重力場(chǎng)；Vzz(x,z-h)和Vzz(x,z-2h)分別表示由觀測(cè)面向上延拓h和2h高度的重力垂向一階導(dǎo)數(shù)；Vzz(x,z0)表示觀測(cè)面上的重力垂向一階導(dǎo)數(shù)。由式(2)可知，計(jì)算下延重力場(chǎng)時(shí)需要利用向上延拓所得到的重力場(chǎng)及其垂向一階導(dǎo)數(shù)。

這里選擇在波數(shù)域中進(jìn)行向上延拓。向上延拓是一種穩(wěn)定算法，對(duì)高波數(shù)成分(包括噪聲)有壓制作用。將觀測(cè)面上的重力異常向上延拓nh高度的表達(dá)式為

(3)

式(2)右端涉及向上延拓和垂向?qū)?shù)兩個(gè)計(jì)算過(guò)程，其中向上延拓具有壓制隨機(jī)噪聲的作用，而垂向求導(dǎo)有提高分辨率的作用，因此該算法比較穩(wěn)定，也具有較高的精度。根據(jù)張沖等[15]的理論模型試驗(yàn)結(jié)果，下延超過(guò)5倍以上點(diǎn)距深度時(shí)，Milne法的計(jì)算誤差顯著低于積分迭代法和傳統(tǒng)波數(shù)域向下延拓方法，即使下延10倍點(diǎn)距時(shí)，下延結(jié)果也不會(huì)發(fā)散。

式(2)中的Vzz(x,z-2h)和Vzz(x,z-h)是重力場(chǎng)經(jīng)過(guò)向上延拓后再計(jì)算其垂向一階導(dǎo)數(shù)，因而對(duì)隨機(jī)噪聲不敏感；Vzz(x,z0)是觀測(cè)面的重力垂向一階導(dǎo)數(shù)，一般情況下，它是由觀測(cè)面上的重力異常直接求垂向?qū)?shù)所得，因此受隨機(jī)噪聲影響最大。故本文選擇對(duì)噪聲有一定壓制作用的垂向?qū)?shù)算法，即ISVD法。

1.3 ISVD法的基本原理

二度體引起的重力位二階導(dǎo)數(shù)，在場(chǎng)源外滿足拉普拉斯方程

Vxx+Vzz=0

(4)

式中：Vxx表示重力位的x方向二階導(dǎo)數(shù)；Vzz表示重力位的垂向二階導(dǎo)數(shù)。通常實(shí)測(cè)數(shù)據(jù)是重力異常，為了獲得重力位，可先在波數(shù)域?qū)χ亓Ξ惓z積分得到重力位V

(5)

(6)

式中：Vx(0)表示三點(diǎn)滑動(dòng)窗口內(nèi)中心點(diǎn)重力位的水平一階導(dǎo)數(shù)；V(1)和V(-1)分別表示滑動(dòng)窗口右端點(diǎn)和左端點(diǎn)的重力位； Δx表示點(diǎn)距。將第一次使用式(6)計(jì)算得到的Vx再次代入式(6)，可求得重力位V在x方向的二階導(dǎo)數(shù)Vxx。最后將Vxx代入式(4)即得到重力異常垂向一階導(dǎo)數(shù)Vzz。由于該方法先進(jìn)行積分運(yùn)算，在一定程度上平滑了重力異常中的隨機(jī)噪聲。值得注意的是，該方法并不能完全抑制數(shù)據(jù)中的非隨機(jī)噪聲，例如淺層地質(zhì)體引起的干擾。當(dāng)重力異常中含有“毛刺”等明顯干擾時(shí)，建議先去除高波數(shù)的干擾成分，然后再做后續(xù)處理。

1.4 基于向下延拓Milne法的重力歸一化總梯度法的實(shí)現(xiàn)步驟

③將第①步和第②步計(jì)算得到的結(jié)果代入式(2)，求出下延深度為h的重力異常Vz(x,zh)。

④再利用第②步方法，計(jì)算Vz(x,zh)的垂向一階導(dǎo)數(shù)Vzz(x,zh)，并應(yīng)用式(6)直接在空間域計(jì)算Vz(x,zh)沿x方向的一階導(dǎo)數(shù)Vzx(x,zh)。

⑤將Vzz(x,zh)和Vzx(x,zh)代入式(1)計(jì)算GH(x,zh)。

⑥修改下延深度，重復(fù)以上步驟，獲得不同深度的重力歸一化總梯度值。

⑦最后繪制GH(x,z)等值線剖面，極大值點(diǎn)即是異常體的中心位置。

2 理論模型實(shí)驗(yàn)

設(shè)計(jì)一個(gè)無(wú)限長(zhǎng)水平圓柱體模型，剩余密度為-0.5×103kg/m3，中心位置為x=100m、z=25m，觀測(cè)剖面長(zhǎng)度為200m，測(cè)量點(diǎn)距為1m，共201個(gè)觀測(cè)點(diǎn)。計(jì)算該模型的重力異常，分別利用本文方法和基于泰勒級(jí)數(shù)展開(kāi)的重力歸一化總梯度法計(jì)算其中心位置，結(jié)果如圖1所示。

圖1 無(wú)限長(zhǎng)水平圓柱體理論模型不同方法重力歸一化總梯度計(jì)算結(jié)果(a)理論重力異常曲線；(b)本文方法；(c)泰勒級(jí)數(shù)展開(kāi)

由圖1可知，采用兩種不同的重力歸一化總梯度法，其極大值點(diǎn)(圖1b和圖1c中的紅色圓點(diǎn))均能準(zhǔn)確地指示出水平圓柱體模型的中心位置，說(shuō)明這兩種方法處理理論模型數(shù)據(jù)的效果較好。

實(shí)測(cè)重力數(shù)據(jù)中通常包含區(qū)域場(chǎng)和隨機(jī)誤差。為檢驗(yàn)該方法的抗干擾能力，對(duì)模型數(shù)據(jù)加入5%的隨機(jī)噪聲和一階趨勢(shì)背景場(chǎng)，處理結(jié)果如圖2所示。

由圖2可知，基于向下延拓Milne法的重力歸一化總梯度法所得到的模型中心位置為x=100m、z=28m(圖2b)，與理論坐標(biāo)位置相比，計(jì)算結(jié)果偏深3m，水平位置無(wú)偏移，說(shuō)明隨機(jī)噪聲和背景場(chǎng)仍然對(duì)該方法有一定的影響。同樣，基于泰勒級(jí)數(shù)展開(kāi)的重力歸一化總梯度法也受影響，計(jì)算埋深偏大1m，水平位置向左偏移2m(圖2c)。對(duì)比圖2b與圖2c 可見(jiàn)，后者出現(xiàn)很多等值線圈閉。在實(shí)際資料處理中，這些虛假異常圈閉容易引起對(duì)真實(shí)場(chǎng)源位置的誤判。引起這些圈閉的原因是基于泰勒級(jí)數(shù)展開(kāi)的重力歸一化總梯度法在計(jì)算過(guò)程中放大了高波數(shù)信號(hào)。相比而言，基于向下延拓Milne法的重力歸一化總梯度法在一定程度上壓制了高波數(shù)噪聲，計(jì)算結(jié)果穩(wěn)定，虛假異常圈閉相對(duì)較少。

選取圖2a剖面兩端的數(shù)據(jù)(0～20m及180～200m)進(jìn)行一階多項(xiàng)式擬合，用擬合值作為區(qū)域異常；再用重力異常減去區(qū)域異常得到局部異常，并對(duì)其進(jìn)行圓滑濾波處理，得到圖3a中虛線所示的異常曲線；最后分別用基于向下延拓Milne法的重力歸一化總梯度法和基于泰勒級(jí)數(shù)展開(kāi)的重力歸一化總梯度法得到圖3b和圖3c所示結(jié)果。由圖3b可見(jiàn)，計(jì)算得到的模型中心位置為x=98m、z=25m，中心埋深與理論模型一致，水平方向有2m的左向偏移；由圖3c可見(jiàn)，計(jì)算得到的模型中心位置為x=98m、z=26m，中心埋深比實(shí)位置偏大1m，水平方向向左偏移2m。這兩種方法均出現(xiàn)水平位置的偏差，分析應(yīng)該與區(qū)域背景場(chǎng)分離不徹底有關(guān)，因?yàn)榉蛛x后的局部重力異常的極小值也有微小的左向偏移。

圖2 對(duì)含5%噪聲和一階趨勢(shì)背景場(chǎng)的理論模型數(shù)據(jù)利用不同方法得到的重力歸一化總梯度剖面(a)重力異常曲線； (b)本文方法； (c)泰勒級(jí)數(shù)展開(kāi)

圖3 對(duì)局部重力異常采用不同方法得到的重力歸一化總梯度剖面(a)局部重力異常曲線； (b)本文方法； (c)泰勒級(jí)數(shù)展開(kāi)

由模型試驗(yàn)可知，利用本文方法計(jì)算重力歸一化總梯度具有中心埋深計(jì)算準(zhǔn)確、虛假異常圈閉相對(duì)較少等優(yōu)點(diǎn)；但若應(yīng)用于含噪、含背景場(chǎng)干擾的重力數(shù)據(jù)處理時(shí)，需要進(jìn)行背景場(chǎng)的剝離和去噪處理，這樣才能得出比較準(zhǔn)確的場(chǎng)源中心位置。

3 實(shí)測(cè)重力剖面數(shù)據(jù)應(yīng)用

將本文方法應(yīng)用于遼寧省葫蘆島市楊家杖子鎮(zhèn)經(jīng)濟(jì)技術(shù)開(kāi)發(fā)區(qū)黑魚溝A礦洞的實(shí)測(cè)重力剖面數(shù)據(jù)。A礦洞位于黑魚溝西山寒武系下方，是鉛鋅礦開(kāi)采礦洞。礦洞斷面為馬蹄形，水平延伸近50m，可近似為一個(gè)無(wú)限長(zhǎng)的水平圓柱體。礦洞中心埋深為5m，高和寬均約2m。在地面采集重力數(shù)據(jù)，重力觀測(cè)剖面垂直于礦洞走向。

觀測(cè)布格重力異常如圖4a中的實(shí)線所示。首先對(duì)該數(shù)據(jù)進(jìn)行異常分離，這里采用非線性濾波分離區(qū)域(背景)場(chǎng)(圖4a中虛線)； 去掉背景場(chǎng)后的局部重力異常如圖4b中實(shí)線所示，可見(jiàn)局部重力異常含有明顯的高波數(shù)干擾成分，故對(duì)其進(jìn)行去噪處理。選用多次圓滑方法壓制隨機(jī)噪聲，去噪后的結(jié)果如圖4b中的虛線所示，對(duì)其分別利用本文方法和基于泰勒級(jí)數(shù)展開(kāi)得到如圖4c和圖4d所示的重力歸一化總梯度剖面。對(duì)比圖4c與圖4d可見(jiàn)，根據(jù)本文方法計(jì)算結(jié)果解釋的礦洞深度與實(shí)際情況基本一致，而根據(jù)圖4d解釋的礦洞埋深較真實(shí)位置偏大1.5m。對(duì)比其他方法計(jì)算的該礦洞中心埋深(5.0～5.6m)[22-25]可以證實(shí)，本文方法處理結(jié)果的準(zhǔn)確度較高。

圖4 實(shí)測(cè)重力數(shù)據(jù)計(jì)算結(jié)果對(duì)比(a)實(shí)測(cè)布格重力異常和分離的區(qū)域(背景)場(chǎng)； (b)局部重力異常(實(shí)線)及去噪后的局部重力異常(虛線)； (c)本文方法得到的重力歸一化總梯度剖面； (d)基于泰勒級(jí)數(shù)展開(kāi)的重力歸一化總梯度剖面 圖中紅色馬蹄形實(shí)線為礦洞位置，紅色圓點(diǎn)為重力歸一化總梯度剖面中極大值點(diǎn)(即礦洞的中心點(diǎn))

4 結(jié)論與建議

理論模型和實(shí)際數(shù)據(jù)計(jì)算結(jié)果表明，基于向下延拓Milne法的重力歸一化總梯度法具有算法穩(wěn)定、準(zhǔn)確性較高的優(yōu)點(diǎn)。相對(duì)于其他重力歸一化總梯度法，不需要設(shè)置最佳諧波數(shù)、正則化參數(shù)及迭代次數(shù)等參數(shù)，減少了數(shù)據(jù)處理中參數(shù)選取引入的誤差。但原始重力數(shù)據(jù)中的隨機(jī)噪聲和背景場(chǎng)對(duì)本方法有一定的影響，故在實(shí)際應(yīng)用中需首先進(jìn)行區(qū)域場(chǎng)的分離和去噪處理。