重力梯度歐拉反褶積及其在文頓鹽丘的應用

侯振隆 王恩德 周文納 吳國超

(①東北大學深部金屬礦山安全開采教育部重點實驗室，遼寧沈陽 110819； ②東北大學資源與土木工程學院，遼寧沈陽 110819； ③蘭州大學地質科學與礦產資源學院，甘肅蘭州 730000； ④浙江大學地球科學學院，浙江杭州 310058)

0 引言

重力勘探是重要的物探方法之一，觀測數據可用于地下空間物性分布[1]和場源位置[2]的計算。歐拉反褶積(Euler Deconvolution,ED)[3]是一種重磁勘探數據處理方法，可以計算地質體的埋深、位置等參數，被廣泛應用于各類資源勘查工作中。相比于傳統重力異常數據，重力梯度數據精度較高，信噪比也更高，包含更多的地質信息，在對地質體的邊界[4]、形狀識別[5]上有較好的應用效果。Zhang等[6]提出利用重力全張量梯度數據進行歐拉反褶積，能夠更充分利用梯度數據。針對傳統歐拉反褶積中反演結果多解性強和空間分辨率較低等問題，一些學者進行方法改進，包括引入基于水平梯度率[7]、均衡邊界檢測濾波器[8]、梯度反褶積[9]、垂向一階導數[10]和解析信號比值[11]等；聯合不同高度觀測面上數據的歐拉反褶積，也被應用于實際勘探中[12-13]。以上研究均一定程度上改善了歐拉反褶積的計算效果。

本文在傳統方法基礎上提出重力梯度數據聯合的歐拉反褶積(Joint Euler Deconvolution of Gravity Gradiometry Data,JEDGG)法，通過計算垂向導數重新推導了反演方程，聯合多個梯度分量建立了新的方程組，避免了構造指數的選擇以及不同分量間的換算。模型試驗驗證了該方法的有效性，并將其應用于文頓鹽丘地區的實測重力數據，在蓋巖東部發現了一處小型構造，進一步探明了該區地下地質結構。

1 重力梯度聯合歐拉反褶積原理

根據歐拉齊次方程的定義[3]，傳統歐拉反褶積公式為

=N(B-f)

(1)

式中： (x,y,z)和(x0,y0,z0)分別代表觀測點和地質體中心坐標， 其中z0表示異常體中心深度；f和B分別表示觀測異常場和背景場；N為構造指數，如薄板的N取0.5，而球體取2.0。將上式轉化為向量積的形式

(2)

用引力位V在x、y和z三個方向的一階導數分別替代式(1)中的f，有

(3)

由于實際數據中缺少Vx和Vy，只能通過分量轉換得到，故式(3)不適于高精度重力梯度數據的計算。對式(3)兩端分別計算z方向上的導數，可得

(4)

整理式(4)，合并同類項，將三階導數轉化為計算二階導數的一階偏導數

(5)

式中Vxz、Vyz和Vzz為重力梯度數據，計算其導數后利用式(5)可求解地質體坐標。需要注意的是，理論上計算z方向導數不會對結果造成影響，但計算實際數據時，求導數相當于濾波，將會放大噪聲。因此，本文將二階導數的一階垂向偏導數轉化為一階導數的二階偏導數。由于引力位滿足拉普拉斯方程，將式(5)中引力位在z方向的三階偏導數轉化為

(6)

將式(6)代入式(5)，轉化為矩陣形式

(7)

由式(7)可知，本文方法避免了構造指數的選擇以及計算垂向導數引起的誤差，計算形式更簡單。聯合Vxz、Vyz和Vzz，減少了分量轉換步驟。與傳統方法一樣，梯度聯合歐拉反褶積也需要使用滑動窗口，窗口中所有測點的方程構成方程組。

2 模型測試

為測試方法的有效性，在笛卡爾坐標系下建立理論模型并進行試驗。試驗包含不加噪聲、加入5%高斯噪聲兩種情況，并對比傳統ED方法的結果。試驗中長方體密度均設為1.0g/cm3。根據式(2)和式(7)，可分別計算ED和JEDGG的結果。

2.1 單長方體模型

首先設地下一長方體的頂面埋深為100m，三維尺寸為800m×800m×200m，觀測點距為20m，異常體水平位置及理論異常見圖1。滑動窗口的大小為19×19，即窗口內包含19×19個測點。試驗結果見圖2～圖5。

從圖2、圖3可以看出，在無噪數據試驗中，JEDGG結果與理論模型基本一致，長方體的水平邊界被清晰地顯示，三維結果可較準確地揭示長方體的埋深；相對地，ED結果的聚焦程度低于JEDGG，其反映的目標體的水平和深度方向上的范圍均比真實值大，并且水平方向上邊界不連續，僅能反映異常體頂點的位置。在含高斯噪聲數據試驗中，其結果(圖4、圖5)與無噪數據計算結果(圖2、圖3)相似，說明JEDGG適用于含噪聲數據的反演。

2.2 雙長方體模型

為了驗證方法對不同埋深地質體模型的計算效果,建立一個雙異常體模型進行試驗。設地下有兩個完全相同的長方體，大小均為500m×500m×100m，頂面埋深分別為100m和200m，觀測點距為20m，這兩個異常體的平面位置及理論重力梯度異常場見圖6。如果JEDGG能夠同時準確地反演出這兩個長方體的空間位置，則驗證了方法的有效性。計算中滑動窗口大小設為11×11，對該模型分別應用JEDGG和ED方法計算，結果見圖7～圖10。

圖1 單長方體模型產生的理論重力梯度異常場(上)及含5%高斯噪聲的重力梯度異常場(下)

圖2 單長方體模型正演數據(無噪聲)的JEDGG結果

圖3 單長方體模型正演數據(無噪聲)的ED結果

圖4 單長方體模型正演(含噪聲)的JEDGG結果

圖5 單長方體模型(含噪聲)的ED結果

圖6 雙長方體模型產生的理論重力梯度異常場(上)及含5%高斯噪聲的重力梯度異常場(下)

圖7 雙長方體模型正演數據(無噪聲)的JEDGG結果

從圖7和8可見，JEDGG法能夠同時反演得到兩個長方體的水平位置和埋深，其中深部長方體的范圍比真實值略小。與ED法計算結果相比，JEDGG計算結果的空間分辨率更高，ED反演得到的異常體范圍偏大、邊界不連續。對于含噪聲(圖9、圖10)情況，JEDGG法計算結果與理論值基本一致，僅在頂面頂點處出現了些許虛假解，即頂點的上方出現“發散狀”的解，但不影響對長方體位置的判斷； ED計算雖然沒有這樣的解，但是長方體范圍仍偏大。這個試驗也驗證了JEDGG處理含噪數據的有效性。

理論數據和含噪數據試驗結果都證明了JEDGG反演效果優于ED，結果更準確，且更適用于較復雜的地質情況。

圖8 雙長方體模型正演數據(無噪聲)的ED結果

圖9 雙長方體模型正演數據(含噪聲)的JEDGG結果

圖10 雙長方體模型正演數據(含噪聲)的ED結果

3 應用實例

基于上述對JEDGG方法的理論研究與模型試驗，將其應用于文頓鹽丘的實測重力全張量梯度數據。文頓鹽丘位于美國路易斯安那州西南，由于富含油氣資源而成為研究熱點。文頓鹽丘包含巖鹽和上方的蓋巖，實測異常主要來源于蓋巖[14-16]。根據Ennen等[15]和Oliveria等[16]的研究，蓋巖埋深約為160～650m。實測數據的x和y軸的正方向分別代表東、北。數據已經密度為2.2g/cm3的地形改正和低通濾波，選取Vxz、Vyz和Vzz三個分量進行JEDGG計算，結果見圖11。

從圖11可以判斷蓋巖頂的水平范圍，并與其他學者對該蓋巖的研究結果[17-19]對比，發現基本吻合。還可以看出，蓋巖的埋深在西部較小，約為125m，較深的位置出現在南部，約550m；頂部深度由東南分別向南部、東部增加，而蓋巖兩端的深度在減少，使蓋巖形態在空間上呈彎曲的“W”形狀。根據已有的密度分布資料[16-20]可知，蓋巖的密度高值位于東南方向，即圖11c中東西方向442.5～443.0km、深度方向275～425m的區域，以及圖11d中南北向3334.0～3334.5km、深度250～550m的區域。結果表明JEDGG法對高密度區域反演效果較好。一般將剩余密度大于0.15g/cm3的區域解釋為蓋巖。從反演結果可見，JEDGG可以反演出剩余密度為0～0.15g/cm3的異常，即蓋巖外圍部分，蓋巖呈“南高北低”、“東高西低”的形態。因此，JEDGG反演結果進一步完善了蓋巖位置的解釋結果。

另外，在圖11中可見，蓋巖的東部可見一處小規模的、呈塊狀的解，常規密度反演[16-20]中沒有該顯示，初步認為這是一處低剩余密度、塊狀或脈狀的異常體的頂端。另外，在埋深0～150m范圍內還存在少量零星分布的解，分析認為是結果中的干擾成分，不具實際意義。

從圖11d可見，在南北方向3334km處、深度范圍0～100m的區域為蓋巖西部解的投影，不能解釋為蓋巖整體深度范圍。

綜上所述，利用JEDGG反演結果可以判斷文頓鹽丘地下蓋巖的水平范圍：東西方向441.7～442.9km， 南北方向3333.8～3334.8km， 深度125～550m。

蓋巖東部還存在一小型構造，有待在以后的勘探工作中做進一步詳查。

圖11 文頓鹽丘實測重力梯度數據JEDGG反演結果

4 結束語

本文提出了多分量重力梯度數據的歐拉反褶積方法，該方法避免了構造指數的選擇和分量轉換，充分利用了梯度數據的高精度特點，提高了地質體空間位置的反演精度。模型數據試驗證明了該方法能夠較準確地識別模型位置。將該方法應用于文頓鹽丘的實際數據，圈定了蓋巖的范圍并發現蓋巖東部還存在一小規模的地質體，豐富了該區域的地下構造信息。

感謝Bell Geospace Inc提供的實測重力全張量梯度數據。