iTilt-Euler法在重力數據處理及斷裂解釋中的應用

陳青，孫帥，丁成藝，黃小宇，陳浩，申鵬，羅港，魏耀聰

(1.重慶科技學院 石油與天然氣工程學院，重慶 401331；2.復雜油氣田勘探開發重慶市重點實驗室，重慶 401331；3.重慶市二零八地質環境工程勘查設計院有限公司，重慶 400700)

0 引言

重力異常場包含了地表及地下諸多密度不均勻體場源產生的綜合信息，直觀地反映了地下地質體的分布位置、深部構造以及斷裂展布等信息。重力資料解釋中最重要的目的是定性和定量地推斷地下客觀存在的異常體的位置、深度、幾何形態及物性參數的過程。然而，受各種密度不均勻體疊加效應的影響，重力異常平面等值線特征往往不能較好地標識深、淺部地質體的信息，無疑增加了重力資料解釋的難度。因此，迫切需要一種能夠將重力場異常進行自動化或半自動化處理和解釋的方法和技術，以提取更多的有效信息。歐拉反褶積作為一種能自動、快速估算場源位置和深度的方法便應運而生。該方法以歐拉齊次方程為基礎，運用位場異常、空間導數，根據場源形狀選擇構造指數，通過歐拉齊次方程的求解，自動或半自動的圈定地質體的邊界及深度[1-3]。該方法的優點是無需已知場源物性的先驗信息，只需要事先確定與場源性質有關的構造指數，便可快速、有效地推算出地質體的具體位置，尤其適用于大面積位場數據的分析和解釋[3]。理論上，歐拉反演結果對應于重磁異常平面上的特征區域，如線性梯級帶、規則性扭曲的等值線或水平錯動的異常軸等，而斷裂往往出現在這些區域，因此，歐拉反褶積法在圈定場源體邊界和識別斷裂中具有重要意義。

盡管歐拉反褶積方法在位場數據自動反演中得到了廣泛應用，但在實際應用中，受地質體的復雜性和場源體之間的疊加影響，該方法存在反演結果發散、虛假解，以及復雜情況下構造指數的確定等問題，制約了該方法實際應用。為了改善歐拉解的發散性，提高反演精度，國內外眾多學者從不同數據源、構造指數選取、消除虛假解等方面做了大量的工作。

在利用不同數據源改善反演結果的發散性方面，Hsu[4]提出了利用重力高階導數進行歐拉反演，大大減小了反演位置和深度解的擴散；Salem等[5]將歐拉反褶積方法與解析信號法相結合，利用解析信號幅度的極大值直接估算出場源的深度及構造指數；Salem等[6]將歐拉反褶積方法與Tilt 梯度相結合，該方法無需構造指數就能推斷出場源邊界和深度分布，從而避免了因構造指數選取導致歐拉反演解的發散和不穩定；范美寧等[7]通過模型計算，證明了用高階導數或解析信號進行歐拉反演計算，可提高反演結果的收斂性；Ma等[8]推導出了歸一化傾斜角(TDX)和Theta圖的歐拉反褶積形式，該形式與Tilt梯度法的歐拉方程形式類似，不受構造指數選取的影響，提高了反演結果的收斂性。

在提高算法精度方面，Neil 等[9]討論了如何獲得相對穩定的最小二乘解；Fairhead等[10]提出利用拉普拉斯方程濾波法來消除歐拉方程的發散解；Keating[11]通過對誤差函數的歐拉方程進行加權計算來消除假頻信號的干擾；范美寧等[12]把最小二乘問題轉化成了解線性方程的問題，從而減少了計算量，減小了誤差傳遞；Davis和Li[13]討論了利用異常數據的振幅信息來降低噪聲和發散解對場源深度估算的干擾。周勇等[14]采用截斷奇異值分解法解歐拉齊次方程，將異常源邊界及中心歐拉解的截斷誤差最小的解作為最優解。

在有效剔除虛假解方面，Gerovska和Araúzo-Bravo[15]通過微分相似變換提取歐拉解奇異點處的虛假解，并利用歐拉解標準差構造評價標準剔除解集中的虛假解；姚長利等[3]提出水平梯度濾波準則、距離約束評價準則和聚集度約束評價準則等方法，推動了歐拉反褶積實用化研究；Ugalde和Morris[16]提出采用聚類分析和內核密度估算技術來解決歐拉解中虛假解的問題；Salem等[17]采用傾斜角水平總梯度峰值來約束反演數據范圍，有效提高了反演結果的收斂性；Beiki等[18]利用截斷奇異值分解方法對誤差相對較大的歐拉解進行剔除，以提升歐拉反褶積對磁源異常的確定精度。

構造指數的確定在歐拉反褶積法求解過程中至關重要，需要根據場源形狀或有關異常性質的先驗知識來選擇。在實際應用中，利用經驗獲取的構造指數進行歐拉反演，往往導致反演結果不準確或解的發散。在構造指數選取方面，Neil 等[9]提出利用統計方法來推斷出構造指數；Salem等[5]利用解析信號幅度的極大值來自動推算構造指數；Barbosa等[19]提出利用歐拉方程中位場總場值f與背景場B之間最小線性相關來獲取構造指數最佳解；姚長利等[3]采用變構造指數對數據滑動窗口重復掃描，從而使得場源解更集中；郭志宏[20]提出將自動估算場源的構造指數的辦法與異常位置及范圍的自動識別方法相結合，形成了一套智能型的改進歐拉反褶積方法；Salem等[6]提出利用不依賴于構造指數的Tilt-Euler法快速推斷出場源邊界和深部，在此基礎上自動估算出構造指數的分布；魯寶亮等[21]提出通過建立歐拉反褶積的超定方程組，求解出最佳構造指數；曹書錦等[22]將截斷誤差與核密度估計進行相關分析，確定最優構造指數。

為減小由經驗獲取的構造指數計算結果的發散性，本文利用不依賴于構造指數的Tilt-Euler法和改進Tilt-Euler法，進行了正演模型計算和對比分析；同時，采用水平總梯度傾斜角峰值約束反演解，優化計算結果。最后，將其應用于肯尼亞ANZA盆地中部地區重力數據處理中，劃分出研究區的斷裂構造體系，為ANZA盆地的進一步地球物理工作提供依據。

1 改進Tilt-Euler法基本原理

1.1 傾斜角法

傾斜角(tilt Derivative)是為了識別不同埋深的場源體邊界提出的方法，該方法利用總場強f的垂直梯度(VDR)比水平總梯度(THDR)的絕對值的反正切角度[23]，定義為：

(1)

式(1)中：

(2)

其中：?f/?x、?f/?y和?f/?z分別為總場強f沿x、y和z方向的一階導數。根據反正切函數的性質，傾斜角的變化范圍為(-π/2，π/2)，在場源體內為正值，外圍為負值，邊界處為零。對于深部場源，在其垂向導數和水平導數都很小的情況下，兩者的比值仍然可以很大，因此，傾斜角反演結果受地質體埋深的影響很小[23-24]。然而，在區域背景場下，位于傾斜角分母上的水平總梯度(THDR)可能趨近于0，導致傾斜角存在“解析奇點”[25]。因此，劉鵬飛等[26]對傾斜角進行了改進，其數學定義式為：

(3)

其中：分母為場強的解析信號振幅A。式(3)中，利用解析信號振幅替換了分母水平總梯度模，提高了計算結果的穩定性[26]。同時，改進的傾斜角繼承了傾斜角對弱異常提取的優勢，從而能很好識別出埋深不同的隱伏場源信息。

1.2 Tilt-Euler法及iTilt-Euler法

歐拉反褶積的基本公式表示為[1]：

=-N(f-B)，

(4)

式中：f為場源位場異常；x、y、z為觀測點的坐標；x0、y0、z0為場源坐標；N為構造指數(N=1，2，…)，與場源的幾何構造有關，是場源異常強度隨深度變化的衰減率；B稱為區域場或背景場。由此可見，該方法是以歐拉齊次方程為基礎，運用位場異常及其空間導數，結合地質體特定的“構造指數”來確定異常場源的位置和深度。然而，由于實際地質構造的場源類型復雜多樣，構造指數值選擇的正確與否直接影響到了場源深度反演解的準確性和穩定性。

Salem等[17]根據傾斜角公式推導出了基于傾斜角的歐拉反褶積，即Tilt-Euler。該方法的優勢是不依賴于已知場源的構造指數，可方便估算出場源的位置分布，提高了歐拉反褶積方法的實用性。Tilt-Euler方程為[17]：

kxx0+kyy0+kzz0=kxx+kyy+kzz，

(5)

其中：kx，ky，kz分別是沿x、y、z方向的傾斜角的導數，分別表示為[17]：

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

1.3 水平總梯度傾斜角峰值約束法

為了提高反演結果的收斂性，Salem等[17]提出利用一定范圍內的傾斜角水平總梯度峰值來約束網格點，從而有效減少了發散解。本文提出利用水平總梯度傾斜角TAHG峰值對反演數據點進行約束，其表達式為[28]：

(14)

由于反正切的特點，TAHG變換范圍也為(-π/2，π/2)。該方法不僅有效平衡了來自淺源和深源的信息，同時最大值位于場源體邊界。因此，本文利用TAHG的最大值范圍對歐拉反褶積的數據點進行約束，從而提高反演結果的收斂性。

本文提出利用傾斜角水平總梯度峰值來約束數據點，通過改進傾斜角沿x、y、z方向導數的峰值建立歐拉方程組，求解出場源體位置參數。具體步驟如下：

1)計算重力異常在x，y和z方向進行一階、二階導數，即Vzx、Vzy、Vzz、Vzxx、Vzxy、Vzyy、Vzzx、Vzzy、Vzzz；

3)利用式(14)，求得水平總梯度傾斜角TAHG，并提取其峰值；

5)設置一定的滑動窗口，在窗口內建立歐拉方程組，求解出場源體的位置參數x0、y0和z0結果，即水平位置和深度；

6)按一定步長滑動子窗口，重復第5步，直至覆蓋全區；

7)將歐拉反演結果成圖并進行解釋。

2 模型試算

為了驗證iTilt-Euler在位場數據處理中的實際應用效果，本文選取3個剩余密度均為0.5×103kg/m3，但頂面埋深和厚度不同的立方體進行理論模型正演計算，模型三維立體圖如圖1a所示，模型參數見表1。圖1b為組合模型的理論重力異常，網格間距為0.8 km×0.8 km。

表1 組合模型參數

a—組合模型三維立體圖；b—組合模型重力異常；c—添加2%高斯噪聲的組合模型重力異常

由圖1b可以看出，模型1埋深較淺，異常幅值基本能勾勒出場源邊緣位置，而對于埋深較深的模型3，因受場源疊加的影響，異常幅值偏向場源邊界外側，且發生扭動。圖2a和2c分別是傾斜角和改進傾斜角的計算結果，由零值反映的邊界位置來看，兩種方法都增強了埋深較深的弱異常信息，但改進傾斜角有效降低了傾斜角的高幅值畸變成分，使得計算結果更為穩定。水平總梯度傾斜角(圖2e)不僅有效增強了弱異常信息，并且其最大值反映的地質體邊界與實際位置也吻合較好。圖2b和2d分別為傾斜角和改進傾斜角在不加峰值約束下的歐拉反演結果，滑動窗口大小選擇11×11?？梢钥闯?，兩種方法的計算結果較相似，但由于改進傾斜角的邊界位置相對收斂，其歐拉解的連續性更好，虛假解也較少。然而，從iTilt-Euler反演結果來看，淺源地質體的歐拉解雖然連續性較好，但由于受淺層高頻異常的影響，估算深度明顯比實際深度小，同時使得深源地質體的歐拉結果出現了大量小于真實值的虛假解。相對而言，水平總梯度傾斜角及其約束之下的iTilt-Euler法反演結果(圖2f)較為收斂，場源解集中分布在場源體邊界處，較好地反映出了模型體的水平位置和深度。約束下的iTilt-Euler法反演的3個模型體深度結果分別為0.94±0.01 km，2.14±0.02 km和3.27±0.02 km，誤差分別為6%、7%和9%，與理論埋深值吻合度較高。因此，利用水平總梯度傾斜角峰值約束反演數據，可以有效減小場源異常相互疊加的影響，使得深源位置也得到較好的反映。

a—傾斜角；b—Tilt-Euler法反演結果；c—改進傾斜角；d—iTilt-Euler法反演結果；e—水平總梯度傾斜角；f—水平總梯度傾斜角峰值約束下的iTilt-Euler法反演結果

由于實際位場數據包含一定的噪聲，為了進一步檢驗反演結果的穩定性，對組合模型加入了2%高斯噪聲(圖1c)進行計算，并與Tilt-Euler法反演結果進行了比較。加噪模型傾斜角(圖3a)的結果反映的模型邊界較為模糊，且明顯偏向外側，而改進傾斜角(圖3b)的結果顯示出較好的收斂性，其不加峰值約束的歐拉反演結果分別位于圖3b和3d，滑動窗口大小選擇11×11。計算過程中，因導數計算對噪聲較為敏感，因此將各階導數向上延拓1.6 km以平滑噪聲的影響。從計算結果看出，iTilt-Euler法反演的邊界位置比Tilt-Euler法的結果更加清晰，但因受場源體異常相互疊加的影響，淺源和深源地質體的結果均出現較多雜解。尤其對深源地質體來講，因受到淺源高頻異常的干擾，鄰近淺源體一側深度解明顯偏大，而另一側則明顯偏小。圖3e和3f分別為水平總梯度傾斜角及其峰值約束之下的歐拉反演結果?？梢钥闯?，約束下的iTilt-Euler反演結果收斂性較好，虛假解明顯減少，深度結果分別為1.12±0.01 km、2.28±0.02 km和3.33±0.03 km，誤差分別為12%、14%和11%。

a—傾斜角；b—Tilt-Euler法反演結果；c—改進傾斜角；d—iTilt-Euler法反演結果；e—水平總梯度傾斜角；f—水平總梯度傾斜角峰值約束下的iTilt-Euler法反演結果

3 實際數據應用

為了驗證iTilt-Euler法的實際應用效果，本文采用肯尼亞ANZA盆地某區塊的重力數據進行了處理分析。從研究區布格重力異常圖(圖4)可以看出，重力高、重力低成排相間分布，整體呈NW走向。研究區東北角發育一NW向展布、未閉合的重力高，幅值約為-50×10-5m/s2。MATASADE與BARCHUMA之間，向北西延伸為一條帶狀展布的重力高，有兩個重力高中心點，異常最大幅值位于MATASADE以西，約為-30×10-5m/s2。DUMA井-MATASADE一線以東至NDOVU井，發育一NW向的重力低，分別在NDOVU井以西和DUMA井東南有兩個重力低圈閉中心。此外，研究區西部為一NW向展布的不規則形重力低，該重力低圈閉中心異常值大約為-90×10-5m/s2。

圖4 研究區布格重力異常

重力異常平面圖上的線性梯級帶、等值線的規則性扭曲或異常軸的水平錯動、串珠狀異常等的分布規律為斷層結構的解釋提供了依據。研究區垂向二階導數(圖5a)和水平總梯度異常(圖5b)清晰反映出研究區異常主體呈NW走向，且異常高、低成排相間分布。圖5c～e分別為利用傳統歐拉反褶積、Tilt-Euler和iTilt-Euler法計算所得的結果，滑動窗口均選擇11×11。反演結果表明，Tilt-Euler和iTilt-Euler法所得的結果更為收斂，且解的分布與垂向二階導數和水平總梯度異常圖中斷裂識別標志較吻合。解的分布主要呈NW向，其次是NE向，反映了研究區斷層的主要展布方向。圖5f為采用水平總梯度傾斜角峰值約束下的iTilt-Euler法反演結果，與無約束的iTilt-Euler結果(圖5e)相比，由于數據點的減少，歐拉解的連續性有所降低，但是研究區中部的局部小斷裂得到了較好的體現，且較深的部分斷裂連續性得到增強。

a—垂向二階導數；b—水平總梯度；c—常規歐拉反褶積反演結果(N=0.5)；d—Tilt-Euler反演結果；e—iTilt-Euler反演結果；f—約束下的iTilt-Euler法反演結果

根據ANZA盆地某區塊重力數據的歐拉反演結果，結合解在重力異常平面圖上的線性梯級帶、等值線的規則性扭曲或異常軸的水平錯動、串珠狀異常等標志帶上的分布規律，推斷出研究區斷裂15條(圖6)。這些斷裂的走向主要可分為近NW向和近NE向兩組，其中，近NW向斷裂6條，分別為F1、F2、F3、F4、F5、F6，近NE向斷裂6條，為F7、F8、F9、F10、F11、F12。此外，還劃分出區域內相互切割的一組次級小斷裂F13、F14和F15。斷裂走向大體分為NW向和NE向兩組，其中以NW向為主，與區域構造走向一致。

圖6 研究區斷裂展布

F1斷裂：該斷裂位于研究區東北角，NNW走向。布格重力異常圖上反映為斷續分布的重力梯級帶，垂向二階導數圖上表現為重力高與重力低的過渡帶，水平總梯度圖上表現為異常的極值連線。約束下的iTilt-Euler法反演結果顯示，斷裂北西段埋深較淺，約3 km，東南段埋深較深，達到8 km左右。斷裂東側埋深較淺，西側埋深較深，故傾向近WS向。

F2-F5斷裂：斷裂位于研究區中東部，NNW—NW走向，是研究區東北部凹陷西邊界的主控斷裂。斷裂的重力場異常標志非常明顯，布格重力異常圖上反映為密集的重力梯級帶，垂向二階導數圖上均表現為明顯的重力高與重力低的過渡帶，在水平總梯度圖上表現為異常的極值連線?；诓煌瑪祿臍W拉反演解連續性均較好，反映的埋深相對較淺。其中，F3、F4、F5斷裂被近NE走向的F10斷裂所切斷，推測其斷裂形成時間可能均早于NE向構造。

F6斷裂：該斷裂是研究區西部重力低值區的東部邊界斷裂，呈NNE走向。布格重力異常圖、垂向二階導數圖、水平總梯度圖等圖件上均表現為明顯的斷裂構造特征顯示。歐拉反演結果顯示，斷裂西南側埋深較深，東北側埋深較淺，故傾向WS。F6斷裂被NE走向的F7、F9斷裂所切斷，因此，推測其斷裂形成時間可能早于NE向構造。該斷裂深度較大，約6～7 km，為控凹斷層，它們控制了西部凹陷內沉積層的厚度和范圍。

F9斷裂：該斷裂位于研究區中部，呈近EW走向。布格重力異常圖上表現為異常等值線扭曲，垂向二階導數圖上表現為兩重力高的鞍部，水平總梯度圖上表現為異常的極值連線。歐拉反演結果顯示，斷裂北側埋深較深，南側埋深較淺，故傾向N。該斷裂連續性較差，且被NW走向的F6斷裂所切割，因此，推測該斷裂形成時間可能晚于NW向構造。

F10斷裂：該斷裂位于研究區東南部，由南端的近SN向NE向延伸。布格重力異常圖上反映為斷續分布的重力梯級帶，垂向二階導數圖上表現為重力高與重力低的過渡帶，水平總梯度圖上表現為極大值連線。歐拉反演結果顯示，該斷裂埋深較大，約在6～9 km的范圍內。該斷裂西北側埋深較淺，南東側較深，故傾向SE。該斷裂連續性較好，切割了NW走向的F3、F4和F5斷裂和NE向的F9斷裂，因此，推測其可能是一條最晚形成、控制區域構造單元東南邊界的深大斷裂。

據20世紀90年代的少量地震反射剖面[29-30]資料來看，NW—SE向展布的地塹系是Anza盆地的主導構造，且均截切基底。此外，據前人[31]對穿過研究區的一條二維地震測線(KAISUT井—NDOVU井一線)的解釋說明，研究區前寒武基底深度差異較大，整體呈現東、西兩側凹陷，而中部隆起的特征。東部凹陷基底深度約為10 km，西部凹陷約為7 km，中間隆起僅約3 km。同時，該地震測線的資料也說明研究區NW向正斷層非常發育，且截切基底，而大量的次生斷裂發育于前寒武結晶基底之上。這一特征與本文歐拉反褶積反演的深度結果一致，研究區NW向基底斷裂在東部和西部區域的截切深度較大，約為6～9 km，而中部區域則較淺，約3～4 km左右。

4 結論

iTilt-Euler法不依賴于場源體的構造指數，能快速估算出場源體的邊界和深度位置，為位場資料的處理和解釋提供了重要的方法手段。相較于Tilt-Euler法，iTilt-Euler法避免了各方向導數求解中存在的“解析奇點”，保證了計算結果的穩定性。同時，本文采用水平總梯度傾斜角峰值約束法，有效約束了反演數據點，使得反演結果更加收斂、準確。模型試算和實際數據應用均反映出，在約束法控制下的iTilt-Euler法的反演結果，由于減少了場源異常相互疊加的影響，不僅有效壓制解的發散性，也使得深部場源的位置和深度信息得到了更好的反映。

研究區斷裂發育，斷裂走向主要可以分為NW向和NE向兩組。NW向斷裂規模大，延伸距離長，切穿深度大，大多為控制區內構造單元邊界的基底斷裂。大多NE向斷裂規模相對小，它們一般切斷NW向斷裂，為蓋層(沉積層內部)斷裂。而在研究區東南部發育的一條近NNE向展布的斷裂，同時切割了NW向和NE向斷裂，推測其可能是控制區域構造單元東南邊界的一條深斷裂。

致謝：感謝西安石油大學袁炳強教授提供的重力數據。