基于向下延拓Milne法的重力歸一化總梯度法

石甲強 肖 鋒 鐘 煬

(吉林大學地球探測科學與技術學院,吉林長春 130026)

0 引言

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

張沖等[15-17]針對向下延拓問題，提出了基于求解微分方程的Adams-Bashforth公式法和基于Simpson求積公式的Milne法。前者下延深度能達到5倍點距，后者甚至可超過10倍點距。Milne法利用觀測面上的重力垂向一階導數、向上延拓不同高度的重力場和重力垂向一階導數求解下延深度的重力場。向下延拓Milne法具有下延深度大、相對誤差小的特點，能獲得穩定準確的結果。但該方法需要計算垂向導數，且求導的精度直接影響向下延拓的精度。

早在2001年，Fedi等[18]基于重力位的二階導數在場源外滿足拉普拉斯方程的原理，提出了對隨機噪聲有壓制作用的積分垂向二階導數法(Integrated Second Vertical Derivative,ISVD)。之后，崔瑞華等[19]將該方法引入中國，首先對重力異常積分，計算得到重力位，然后計算重力位的水平二階導數，最后根據拉普拉斯方程計算重力異常的垂向一階導數。相對于常規的波數域求垂向導數的方法，ISVD法由于進行了積分運算，對隨機噪聲有抑制作用，因而更穩定。

本文將這兩種穩定算法，即向下延拓Milne法和ISVD法引入重力歸一化總梯度的計算，通過理論模型試算和實際資料處理，驗證了方法的有效性。

1 方法原理

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

以重力剖面觀測數據為例，歸一化總梯度的表達式[2]為

(1)

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

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

重力場向下延拓四階Milne公式[15]為

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

(2)

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

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

(3)

式(2)右端涉及向上延拓和垂向導數兩個計算過程，其中向上延拓具有壓制隨機噪聲的作用，而垂向求導有提高分辨率的作用，因此該算法比較穩定，也具有較高的精度。根據張沖等[15]的理論模型試驗結果，下延超過5倍以上點距深度時，Milne法的計算誤差顯著低于積分迭代法和傳統波數域向下延拓方法，即使下延10倍點距時，下延結果也不會發散。

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

1.3 ISVD法的基本原理

二度體引起的重力位二階導數，在場源外滿足拉普拉斯方程

Vxx+Vzz=0

(4)

式中：Vxx表示重力位的x方向二階導數；Vzz表示重力位的垂向二階導數。通常實測數據是重力異常，為了獲得重力位，可先在波數域對重力異常Vz積分得到重力位V

(5)

(6)

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

1.4 基于向下延拓Milne法的重力歸一化總梯度法的實現步驟

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

④再利用第②步方法，計算Vz(x,zh)的垂向一階導數Vzz(x,zh)，并應用式(6)直接在空間域計算Vz(x,zh)沿x方向的一階導數Vzx(x,zh)。

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

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

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

2 理論模型實驗

設計一個無限長水平圓柱體模型，剩余密度為-0.5×103kg/m3，中心位置為x=100m、z=25m，觀測剖面長度為200m，測量點距為1m，共201個觀測點。計算該模型的重力異常，分別利用本文方法和基于泰勒級數展開的重力歸一化總梯度法計算其中心位置，結果如圖1所示。

圖1 無限長水平圓柱體理論模型不同方法重力歸一化總梯度計算結果(a)理論重力異常曲線；(b)本文方法；(c)泰勒級數展開

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

實測重力數據中通常包含區域場和隨機誤差。為檢驗該方法的抗干擾能力，對模型數據加入5%的隨機噪聲和一階趨勢背景場，處理結果如圖2所示。

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

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

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

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

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

3 實測重力剖面數據應用

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

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

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

4 結論與建議

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