朱自強，王燦，魯光銀，曹書錦

?

重力梯度張量解析信號的歐拉反褶積

朱自強，王燦，魯光銀，曹書錦

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

利用重力梯度張量數(shù)據(jù)高精度的特點以及解析信號在確定異常體位置上的優(yōu)勢，將重力梯度張量解析信號代替位場的導數(shù)完成重力梯度張量解析信號的歐拉反褶積；通過在1個窗口內(nèi)對1組數(shù)據(jù)點解3個歐拉方程來自動識別構(gòu)造指數(shù)，從而規(guī)避了傳統(tǒng)歐拉反褶積方法中需要事先確定構(gòu)造指數(shù)的問題，同時減少了背景場的影響。研究結(jié)果表明：使用重力梯度張量的解析信號，其歐拉反褶積的解收斂性很好，能準確地判斷地下異常體源的位置，有效規(guī)避背景場的影響，反演效果較好。

重力梯度張量；解析信號；歐拉反褶積

重力梯度張量是重力位的二階導數(shù)。傳統(tǒng)的重力測量主要依靠密度的變化，重力梯度張量對微小的密度梯度的變化很敏感，因此，重力梯度張量能夠更加靈敏地反映地下密度異常體，能夠更加直接地反映場源體的邊界范圍[1]。解析信號是梯度值[2]，其振幅能夠反映地下異常體的位置。將解析信號應用于重力場中，可以減少背景場的影響[3]，從而能夠清晰地確定地下異常體的位置[4]。歐拉反褶積是一種能夠自動或半自動地估算場源體的位置的反演方法，由于其在反演時所需使用的參數(shù)少，不需要知道場源密度或磁性等物性的先驗信息，不需要精確的解釋模型[5]，能快速而有效地反映地下異常源的邊界。傳統(tǒng)的歐拉反褶積需要事先確定或者列舉出與場源性質(zhì)有關的構(gòu)造指 數(shù)，但當異常場中有多個源時，構(gòu)造指數(shù)可能是變化的，而且在實際情況中，地下異常場源是未知的，因此，其構(gòu)造指數(shù)難以確定，這樣往往使得歐拉反褶積的解具有多解性和發(fā)散性。為了消除歐拉方程解的發(fā)散性問題，提取正確的歐拉解，F(xiàn)itzGerald[6]等采用了一種新的識別歐拉解方法，但是該方法不能評價歐拉解是否為優(yōu)解。Mikhailov等[7]采用人工智能的方式提取歐拉方程的解。Gerovska等[8]在Stavrev等[9?10]提出的算法的基礎上使用非聚類和聚類的方式處理歐拉解。Ugalde等[11]采用模糊聚類的方法來提取歐拉解，但該法過度地依賴人工來確定地下異常體源的個數(shù)。范美寧[12]提出水平梯度濾波法、基于主體異常距離準則的有效性統(tǒng)計篩選和構(gòu)造指數(shù)的有效性篩選等改進措施[13?14]，將三維歐拉反褶積方法進行推廣應用。曹書錦等[5]通過利用邊緣檢測張量不變量數(shù)據(jù)對異常體邊界進行提取，來輔助剔除歐拉解中的發(fā)散解，取得了很好的效果。本文1個窗口內(nèi)對1組數(shù)據(jù)點解3個歐拉方程[15]來最后通過模型算例進一步對比說明本文方法的有效性。

1 重力梯度張量及其解析信號

重力梯度張量即在笛卡爾坐標下重力位在，和方向的二階導數(shù)，1個體積為、剩余密度為的質(zhì)量體，其重力位表示為

其中：為觀測點到場源之間的距離；為萬有引力常量。重力梯度張量表示為

其二維振幅表示為

Nabighian[16]將二維解析信號推廣到三維，Roset等[17]將位場的解析信號的定義進行擴展如下：

則其振幅可以表示為

對于重力梯度全張量，其解析信號可用矩陣表示為

因此，振幅可以表示為

Debeglia等[18]提出解析信號振幅的一階導數(shù)能更有效地分離目標異常體，對解析信號振幅求一階導數(shù)表示為：

2 重力梯度張量歐拉反褶積

傳統(tǒng)的歐拉反褶積是使用的正交梯度，表示為

由于重力梯度張量是重力異常沿坐標系軸向的導數(shù)，因此，將傳統(tǒng)歐拉方程推廣到重力梯度張量的歐拉反褶積方程用矩陣表示為

其中：，和為背景場，由于背景場很小，故可忽略；，和為重力梯度；，，，，，，，和為重力梯度張量。

3 解析信號的歐拉反褶積

Keating等[19]使用歐拉反褶積有效地去除了解析信號的背景場。重力梯度張量解析信號的歐拉反褶積可以表述為

4 模型算例

以規(guī)則的球體和立方體作為研究模型，分別用重力梯度張量及其解析信號進行歐拉反褶積計算。

4.1 算例1：球體模型

球體模型設置：球體半徑為1 km，球心坐標為 (0, 0, 1.5) km。剩余密度設為1 kg/m3，測網(wǎng)高度為0 m；在測區(qū)，設為?10~10 km，為?10~10 km；測網(wǎng) 長×寬為0.1 km×0.1 km。圖1 所示為球體重力梯度張量解析信號的歐拉反褶積解，圖1中黑點表示歐拉反褶積的解，黑色實線表示球體模型的輪廓線。圖2所示為球體重力梯度張量的歐拉反褶積解。

圖1 球體重力梯度張量解析信號的歐拉反褶積解

圖2 球體重力梯度張量的歐拉反褶積解

對比圖1和圖2可見：對于同樣的球體模型，使用重力梯度張量能較好地反映球體的球心位置，重力梯度張量的解析信號具有唯一解，能準確地描繪出球體場源的中心位置；歐拉解完全處于球體的中心部位。

從球體重力梯度張量解析信號歐拉反褶積結(jié)果可知：每個解在深度上都與模型設定值極吻合，構(gòu)造指數(shù)與理論值完全相符，由此證明了本文方法的正確性。

4.2 算例2：立方體模型

立方體模型設置：立方體的長×寬×高為 0.6 km×0.6 km×0.6 km，質(zhì)心坐標為(0, 0, 1.0) km。剩余密度設為1 kg/m3，測網(wǎng)高度為0 km。在測區(qū)，設為?10~10 km，為?10~100 km；測網(wǎng)長×寬為0.1 km×0.1 km。圖3所示為立方體重力梯度張量解析信號的歐拉反褶積解，圖4所示為立方體重力梯度張量解析信號的歐拉反褶積解，其中黑點表示歐拉反褶積的解，黑色實線表示立方體模型的輪廓線。

對比圖3和圖4可見：對于立方體模型，重力梯度張量的歐拉反褶積解大體匯聚在立方體的場源中心，歐拉反褶積解的底部較準確地反映了立方體的中心位置，但張量的歐拉反褶積解有一些向上延拓的發(fā)散解，對于圈定立方體場源中心可能帶來不利的影響。從圖3可以看出：使用重力梯度張量的解析信號之后，立方體的歐拉反褶積解很好地匯聚在立方體的場源中心，歐拉解完全處于立方體的中心部位。

圖3 立方體重力梯度張量的歐拉反褶積解

圖4 立方體重力梯度張量解析信號的歐拉反褶積解

立方體張量解析信號進行歐拉反褶積的埋深和構(gòu)造指數(shù)分別見圖5和圖6。從圖5和圖6可見：在前100個測點中，埋深大體在1.00 km的位置，上、下浮動在4.12 m之內(nèi)，之后，變化很小；同樣，前100個測點中，構(gòu)造指數(shù)有極小的偏差，隨后與理論值極吻合。這說明了該方法在確定單個模型位置和構(gòu)造指數(shù)中的準確性。

圖5 立方體重力張量解析信號的歐拉反褶積埋深

圖6 立方體重力張量解析信號的歐拉反褶積構(gòu)造指數(shù)

4.3 算例3：立方體疊加模型

將5個立方體疊加，立方體模型的質(zhì)心坐標用 (0，0，0)表示，剩余密度設為1 kg/m3，測網(wǎng)高度為0 m。在測區(qū)，為?10~10 km，為?10~10 km；測網(wǎng)長×寬為0.1 km×0.1 km。模型幾何參數(shù)設置如表1所示。

表1 模型幾何參數(shù)

Table 1 Model geometrical parameters km

立方體序號X0Y0Z0長寬高 11.000 1.0000.1500.0500.0500.025 20.650 0.6500.2500.1000.1000.100 30.300 0.3000.4000.2000.2000.200 4?0.500?0.5000.9000.4000.4000.600 5?3.100?3.1003.1001.2001.2002.400

立方體疊加模型重力梯度張量分量的等值線圖見圖7。從圖7可看出該疊加模型正演的正確性。

圖7 立方體疊加模型重力梯度張量gzz分量的等值線圖

對于規(guī)模較大的異常體(模型中的立方體5)，其解在總體歐拉解集中占有比較大的部分，相對較小規(guī)模尤其是靠近大規(guī)模異常體的小規(guī)模異常體(例如組合模型中的立方體3)，其異常源被淹沒。圖8和圖9所示分別為立方體疊加模型的重力梯度張量解析信號歐拉反褶積解、重力梯度張量歐拉反褶積解。從圖8可以看出，根據(jù)重力梯度張量的解析信號進行的歐拉反褶積解能準確、充分地判斷出地下各個異常體源的位置情況。

圖8 立方體疊加模型重力張量解析信號的歐拉反褶積解

圖9 立方體疊加模型重力張量歐拉反褶積解

據(jù)圖9中僅可大致估計異常體的深度。

5 結(jié)論

1) 簡單規(guī)則模型的重力梯度張量解析信號的歐拉解及其構(gòu)造指數(shù)表明本文方法的正確性和可行性。

2) 重力梯度張量的解析信號能大大提高歐拉反褶積解的收斂性，改善了反演效果。這是由于本文算法中使用了重力位的三階導數(shù)，而傳統(tǒng)的重力張量歐拉反褶積只用到重力位的二階導數(shù)。

3) 使用重力梯度張量解析信號能減少異常源間歐拉解的相互影響和干擾，因此，使用重力梯度張量的解析信號能有效規(guī)避背景場的影響。

4) 使用重力梯度張量的解析信號在1個窗口解3個歐拉方程，能在解的過程中自動識別構(gòu)造指數(shù)，規(guī)避了傳統(tǒng)歐拉反褶積方法中需要事先確定構(gòu)造指數(shù)的問題，同時大大減少了背景場的影響。

5) 由于我國目前很難得到實測的梯度，解析信號數(shù)據(jù)需要經(jīng)過計算得到，若能夠獲得直接的實測梯度和解析信號，或者使用某些技術(shù)手段來求取梯度和解析信號，壓制其誤差，則反演效果會得到更大改善。為了獲得直接實測梯度數(shù)據(jù)和解析信號數(shù)據(jù)，需要加大對重力梯度儀的改進和推廣使用。

[1] 曾思紅. 重力梯度張量正演研究及邊界提取[D]. 長沙: 中南大學地球科學與信息物理學院, 2010: 1?3.
ZENG Sihong. Gravity gradient tensor forward research and edge extraction[D]. Changsha: Central South University. School of Geosciences and Info-Physics, 2010: 1?3.

[2] 王萬銀. 位場解析信號振幅極值位置空間變化規(guī)律研究[J]. 地球物理學報, 2012, 55(4): 1288?1299.
WANG Wanyin. Spatial variation law of the extreme value position of analytical signal amplitude for potential field data[J]. Chinese J Geophys, 2012, 55(4): 1288?1299.

[3] 張季生, 高銳, 李秋生, 等. 歐拉反褶積與解析信號相結(jié)合的位場反演方法[J]. 地球物理學報, 2011, 54(6): 1634?1641.
ZHANG Jisong, GAO Rui, LI Qiusen, et al. A combined Euler and analytic signal method for an inversion calculation of potential data[J]. Chinese J Geophys, 2011, 54(6): 1634?1641.

[4] Beiki M. Analytic signals of gravity gradient tensor and their application to estimate source location[J]. Geophysics, 2010, 75(6): 59?74.

[5] 曹書錦, 朱自強, 魯光銀. 基于自適應模糊聚類分析的重力張量歐拉反褶積解[J]. 中南大學學報(自然科學版), 2012, 43(3): 1033?1039.
CAO Shujin, ZHU Ziqiang, LU Guangyin. Gravity tensor Euler deconvolution solutions based adaptive fuzzy cluster analysis[J]. Journal of Central South University (Science and Technology), 2012, 43(3): 1033?1039.

[6] FitzGerald D, Reid A, McInerney P. New discrimination techniques for Euler deconvolution[J]. Computers and Geosciences, 2004, 30(5): 461?469.

[7] Mikhailov V, Galdeano A, Diament M, et al. Application of artificial intelligence for Euler solutions clustering[J]. Geophysics, 2003, 64(2): 168?180.

[8] Gerovska D, Araúzo-Bravo M J. Automatic interpretation of magnetic data based on Euler deconvolution with unprescribed structural index[J]. Computers & Geosciences, 2003, 29(8): 949?960.

[9] Stavrev P Y. Euler deconvolution using differential similarity transformations of gravity or magnetic anomalies[J]. Geophysical Prospecting, 1997, 45(2): 207?246.

[10] Stavrev P Y, Gerovska D, Arauzo-Bravo M J. Depth and shape estimates from simultaneous inversion of magnetic fields and their gradient components using differential similarity transforms[J]. Geophysical Prospecting, 2009, 57(4): 707?717.

[11] Ugalde H, Morris W A. Cluster analysis of euler deconvolution solutions: new filtering techniques and geologic strike determination[J]. Geophysics, 2010, 75(3): L61?L70.

[12] 范美寧. 歐拉反褶積方法的研究及應用[D]. 長春: 吉林大學地球探測科學與技術(shù)學院, 2006: 62?63.
FAN Meining. The study and application of Euler deconvolution method[D]. Changchun: Jilin University. Geo-Exploration Science and Technology Institute, 2006: 62?63.

[13] 姚長利, 管志寧, 吳其斌, 等. 歐拉反演方法分析及實用技術(shù)改進[J]. 物探與化探, 2004, 28(2): 150?155.
YAO Changli, GUAN Zhining, WU Qibin, et al. An analysis of Euler deconvolution and its improvement[J]. Geophysical & Geochemical Exploration, 2004, 28(2): 150?155.

[14] 王家林. 對我國石油重磁勘探發(fā)展的幾點思考[J]. 勘探地球物理進展, 2006, 29(2): 82?86.
WANG Jialin. Views on the domestic situation and progress of gravity and magnetic petroleum exploration[J]. Progress in Exploration Geophysics, 2006, 29(2): 82?86.

[15] 范美寧, 江裕標, 張景仙. 不同數(shù)據(jù)用于歐拉方程的模型計算[J]. 地球物理學進展, 2008, 23(4): 1250?1253.
FAN Meining, JIANG Yubiao, ZHANG Jingxian. Model calculation of Euler’s equation for different data types[J]. Progress in Geophysics, 2008, 23(4): 1250?1253.

[16] Nabighian M N. Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: Fundamental relations[J]. Geophysics, 1984, 49(5): 780?786.

[17] Roest W R, Verhoef J, Pilkington M. Magnetic interpretation using the 3-D analytic signal[J]. Geophysics, 1992, 57(1): 116?125.

[18] Debeglia N, Corpel J. Automatic 3-D interpretation of potential field data using analytic signal derivatives[J]. Geophysics, 1997, 62(1): 87?96.

[19] Keating P, Pilkington M. Euler deconvolution of the analytic signal and its application to magnetic interpretation[J]. Geophysical Prospecting, 2004, 52(3): 165?182.

Euler deconvolution of analytic signals of gravity gradient tensor

ZHU Ziqiang, WANG Can, LU Guangyin, CAO Shujin

(School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Depending on high resolution of gravity gradient tensor data as well as the advantage of determining the location of gravity anomalies by the analytic signal, Euler deconvolution of the analytic signal of gravity gradient tensor was solved by analytic signal instead of gravity field derivatives. Structural index was automatically identified by a set of data points in a window solving Euler equations, so there was no need to determine structural index in advance and denoise background field which traditional Euler deconvolution method had. The results show that Euler deconvolution solution converges quickly, the disturbance of background field can be eliminated effectively, and the results of this method are useful and robust.

gravity gradient tensor; analytic signal; Euler deconvolution

P631.1

A

1672?7207(2015)01?0217?06

2014?01?10；

2014?03?12

國家自然科學基金項目資助(41174061) (Project(41174061) supported by the National Natural Science Foundation of China)

魯光銀，博士，教授，從事地震及重磁相關領域的研究；E-mail: csulgy@163.com

10.11817/j.issn.1672?7207.2015.01.029

(編輯 陳燦華)