重力場歐拉反褶積最優解提取

周 勇， 曹書錦， 侯萍萍， 楊子龍

(1. 湖南文理學院 資源環境與旅游學院，常德 415000；2. 湖南科技大學 頁巖氣資源利用湖南重點實驗室，資源環境與安全工程學院，湘潭 411201；3. 中南大學 地球科學與信息物理學院，長沙 410083;4.湖南文理學院 洞庭湖生態經濟區建設與發展湖南省協同創新中心，常德 415000)

重力場歐拉反褶積最優解提取

周 勇1，3，4， 曹書錦2,3， 侯萍萍2， 楊子龍2

(1. 湖南文理學院 資源環境與旅游學院，常德 415000；2. 湖南科技大學 頁巖氣資源利用湖南重點實驗室，資源環境與安全工程學院，湘潭 411201；3. 中南大學 地球科學與信息物理學院，長沙 410083;4.湖南文理學院 洞庭湖生態經濟區建設與發展湖南省協同創新中心，常德 415000)

在傳統重力場歐拉反褶積中，單一構造指數難于表征多個異弱常源；而枚舉構造指數易于導致歐拉解過度發散、歐拉解解集龐大且其中謬解極多。為此，在采用非預測歐拉反褶積方法，避免枚舉構造指數的基礎上，采用截斷奇異值分解法解歐拉齊次方程，以評價歐拉解的不穩定性，在此基礎上，構建估計函數以保留位場總場效應相對較弱的異常源。最后通過數值例子進一步表明了方法的有效性和可靠性。

歐拉反褶積； 最優解； 估計函數； 截斷誤差； 奇異值分解

0 引言

歐拉反褶積是一種快速高效且適應性強的位場解釋方法。其以使用先驗信息少，不需要準確的解釋模型，在某一尺度大小的滑動窗口下，通過人為確定/枚舉構造指數從而快速而有效地圈定出異常源的基本輪廓而著稱。當多個不同形態的異常場并存時，難于使用單一常數或枚舉構造指數表征多個異常源或復雜地質構造異常源，導致通過人工經驗選取構造指數的難度加大。由于歐拉超定方程組條件數極大，易于導致歐拉解發散。因而通過枚舉構造指數，常使得歐拉解解集龐大且其中謬解極多，導致后續分析處理工作量急劇增大[1-2]。因而如何從眾多的歐拉解中分辨和確定最優解一直是困擾人們的一個難題[3-9]。

為有效地剔除歐拉解解集中的謬解和提取最優解[2-3, 10-17]，諸多學者開展了富有成效地研究。姚長利等[12]提出水平梯度濾波準則、距離約束評價準則和聚集度約束評價準則等方法，促進歐拉反演方法進入實用化階段;魯寶亮等[2, 13]利用歐拉反褶積建立的超定方程組，求解出一組構造指數，計算出構造指數的可信范圍和最佳構造指數，然后利用最佳構造指數歐拉反演計算;Gerovska等[1]用微分相似變換對歐拉解奇異點處的空間坐標、構造指數和線性背景場進行分析以提取繆解，并在此基礎上，以歐拉解標準差構造評價標準過濾歐拉解解集中的謬解;Beiki等[18]利用截斷奇異值分解方法對誤差相對較大的歐拉解進行剔除，以提升歐拉反褶積對磁源異常的確定精度;曹書錦等[5]將截斷誤差與核密度估計進行相關分析，確定構造指數的一維核密度估計曲線的第二個峰值處的構造指數為最優構造指數，但該方法并沒有將歐拉解解集作為多維數據來處理。這導致存在相鄰或相似異常源混疊時，難于確定哪些歐拉解標示著同一個異常源。

在歐拉反褶積中，使用單一構造指數，易于導致歐拉解解集趨于發散；而枚舉構造指數易使歐拉解解集過于龐大，導致后續分析、處理工作難度增大。針對如何選取最優解的問題，筆者從歐拉超定方程組的條件數過大的問題出發，計算歐拉解的截斷誤差，分析了不同噪聲水平與截斷誤差的關系，并進一步確定了在異常源的邊界上及其異常源中心歐拉解的截斷誤差最小，即為最優解，通過數值例子進一步表明了方法的有效性。

1 歐拉反褶積反演

假設笛卡爾坐標系于點(x,y,z)存在一個孤立異常源，在觀測點(x0,y0,z0)處的重力異常響應f可以寫為：

(1)

式(1)滿足歐拉齊次方程，則

(2)

式中：?f/?x、?f/?y和?f/?z為f沿笛卡爾坐標系X、Y和Z三軸向的梯度，一般由離散傅里葉變換或離散余弦變換獲得；B為消除區域背景場的影響，而引入的一個代表區域背景異常的參數。

將式(1)帶入式(2)得到式(3)。

Nf(x-x0,y-y0,z-z0)+NB

(3)

特定地質構造有特定的構造指數，在求解式(3)時，可將構造指數作為一個預設參數，采用步長△N在0～3范圍內枚舉構造指數，獲得歐拉解解集。隨著步長△N不斷地變小，歐拉解解集急劇增大，導致后續數據處理工作難度增大。為此，將構造指數N作為待求參數，將式(3)重寫為式(4)。

(4)

式(1)基于剩余密度計算，背景場B數量級很小，故此這里將B略去。利用某一尺度wx的方形滑動窗口在勘察區域滑動，將滑動窗口內n個觀測點的重力異常響應及其梯度，代入到式(4)中構建超定方程組，利用數值方法獲取歐拉解，即異常源的空間位置信息(x,y,z)及構造指數N兩部分。但由于式(4)本質上仍然是歐拉齊次方程，仍然存在對歐拉反褶積優解提取和繆解剔除的工作。

FitzGerald等[3]系統地給出了傳統提取最優解的策略，但在實施上存在諸多限制或需要考慮的因素，如沒有給出可靠的可信度評價原則或奇異值的上下限，也不能估計整個歐拉解解集的優良性或確定哪些相類似的解標示同一個異常源。由于系數矩陣條件數很大，歐拉解解集受噪聲干擾及滑動窗口大小影響很大。為此，筆者引入誤差估計方法對歐拉解進行評價，以對比在不同噪聲和不同位置處的歐拉解的分布情況。

2 基于截斷奇異值分解法的歐拉解

將滑動窗口內n個觀測值f帶入到式(4)中，把歐拉反褶積方程構成的方程組寫為矩陣的形式：

Gm=d

(5)

由于式(5)為超定方程組，使得歐拉解易受如觀測數據噪聲、滑動窗口大小、異常源深度等干擾因素影響而擾動。為對歐拉解的誤差進行估計，利用截斷奇異值分解法對m的誤差進行估計。首先對系數矩陣G進行奇異值分解

G=USVT

(6)

式中：U為n×n的酉矩陣；V為4×4的對角矩陣，其對角線上的四個元素為G的特征值；S=diag(σ1,…,σp,…,0,…)，為n×4的矩陣，且對角線上的奇異值有σ1≥σ2≥…≥σp≥0，p=4；向量u1、…、up和v1、…、vp分別對應于矩陣U與V中各列向量。

將式(6)帶入到式(5)中，歐拉解可寫為

m=V/SUd

(7)

此時，m的誤差可由下式估計得到：

(8)

其中

(9)

其中：b=Ge=m-m*；e為殘差；m*=G-1d。

為評價歐拉解穩定性，特設定一個估計函數以過濾歐拉解解集中的繆解[19-20]

(10)

當歐拉反演對深度估計不準時，如兩個異常源相互干擾，易于使得σr所表示的位置偏離異常源中心。基于截斷誤差的改進的估計函數寫為式(11)。

(11)

當σr大于某一特定的閾值σ0，便認為其過濾得到的歐拉解為最優解。

3 模型試驗

3.1 算法驗證

以重力正演解析解明確且構造指數易于確定的簡單幾何體(如立方體，其構造指數為2)構建地球物理模型，驗證本文歐拉反演算法的正確性。在均勻半空間構造模型I：一個大小為200 m × 200 m × 300 m、頂板埋深為150 m、底板為350 m且剩余密度為0.3 g/cm3的異常體。該模型沿笛卡爾坐標系X、Y和Z三軸向的剖分數分別為32、32和16；測網大小為50 m × 50 m，觀測點個數為32 × 32 = 1 024，觀測點高度為地面上25 m。利用解析解計算簡單規則模型重力異常響應f，進一步利用快速傅里葉變換計算重力響應f沿著笛卡爾坐標系沿三軸向的導數fx、fy、fz。在此基礎上，應用歐拉反褶積方法計算異常源的空間位置信息(x,y,z)及構造指數N兩部分。

如圖 1所示，在歐拉反演解釋中，一般將歐拉解逐點畫于觀測數據上，以標示異常源的分布情況。該立方體異常源的歐拉解均位置異常源中心，且構造指數的值分布于1.91～2.13之間，這與立方體異常源的構造指數的理論值一致，表明了本算法的正確性。

圖1 孤立異常源歐拉解分布示意圖(未添加噪聲)Fig.1 Euler deconvolution solutions of isolated anomaly (without noise)

3.2 孤立異常源

在此基礎上，為提取歐拉反褶積最優解，特從孤立異常源、同一深度和不同深度異常源等三個方面研究歐拉解的分布情況。由于歐拉齊次方程為超定方程，易于受觀測數據所含噪聲的影響,特在模型I正演結果內添加3%的高斯白噪聲。

圖 2為孤立異常源正演數據添加3%高斯白噪聲后的歐拉解分布示意圖。通過與圖 1進行對比，雖然分布于異常體內部的構造指數值趨近于立方體構造指數的理論值，但分布范圍從原有的1.91～2.13變為1.15～2.04。相比于圖 1中在異常源外部沒有任何繆解，而在圖 2中出現大量繆解，表明觀測數據的噪聲將極大地導致歐拉反褶積的發散。

圖2 孤立異常源歐拉解分布示意圖 (添加3%噪聲，未過濾)Fig.2 Euler deconvolution solutions of isolated anomaly (with 3% white Gaussian noise, before filtered)

圖 3為基于截斷誤差估計的孤立異常源歐拉解的閾值示意圖。從圖3中可以看出，在異常源中心和異常源外部σr的差異性非常小，不利于過濾辨別及過濾歐拉反褶積繆解。主要由于當隨著歐拉解遠離異常源中心時(圖 2)，其深度z逐漸變小，這導致由式(10)計算得到估計函數σr在異常源內部和外部的差異非常小。

圖3 估計函數閾值σr示意圖Fig.3 Threshold value of estimation function σr

圖4 改進的估計函數閾值示意圖Fig.4 Threshold value of improved estimation function

圖 4為基于截斷誤差估計的孤立異常源歐拉解的閾值示意圖。相比于圖 3，改進的σr在異常源外部和異常源內部的差異很大，這便于對歐拉解解集中謬解的剔除。

圖5 構造指數及截斷奇異值分解閾值的概率密度估計Fig.5 Probability density estimation of struct index and threshold of truncated singular value decomposition

為確定σ0，利用核密度估計分析截斷奇異值分解閾值曲線和構造指數,并進行相關分析(圖 5)，確定采用閾值σ0>0.2作為過濾歐拉解解中繆解的標準。如圖 6所示，表明歐拉積解恰處于正演模型的內部中。

圖6 歐拉解分布示意圖 (添加3%高斯白噪聲，過濾后)Fig.6 Euler deconvolution solutions (with 3% white Gaussian noise, before filtered)

從圖 1和圖 2中可以看出，在未對歐拉解進行過濾時，淺部的歐拉解受噪聲干擾比較明顯。因此采用估計函數σr以σ0=0.2作為閾值過濾歐拉反褶積繆解。過濾后的歐拉解恰位于異常源中心，這表明了估計函數σr的正確性。相比于傳統基于構造指數的過濾策略，如以N<0、0.2、0.5等為標準，過濾歐拉解解中的繆解，采用估計函數σr以σ0=0.2作為閾值過濾歐拉解繆解，則能保留能標示異常源中心但N不接近2的,因而，本文算法更具有優勢。

3.3 同一深度異常源

由于改進的估計函數σr與深度無關，故而需要對同一深度和不同深度模型進行對比分析。特構如下地球物理模型 II：由兩個均為200 m × 200 m × 200 m，且剩余密度分別為0.3 g/cm3，質心分別為(-400 m, 0 m, 250 m)與(400 m, 0 m, 250 m)的立方體組成。該模型沿著笛卡爾坐標系三軸向的剖分數分別為32、32和16；測網大小為50 m × 50 m，觀測點個數為32 × 32 = 1 024，觀測點高度為地面上25 m。應用歐拉反演計算異常源異常空間(x,y,z)及構造指數N兩部分。

圖7 歐拉解分布示意圖(未過濾)Fig.7 Euler deconvolution solutions (before filtered)

圖8 歐拉解分布示意圖(過濾后)Fig.8 Euler deconvolution solutions (after filtered)

圖 7和圖 8分別為歐拉解過濾前及過濾后分布示意圖。由于兩異常源的位場總場效應相當，歐拉解的分布呈現出一定的對稱趨勢。通過圖7和圖8的對比分析可以發現，采用估計函數σr以σ0=0.2作為閾值過濾歐拉反褶積繆解，能有效地標示兩相鄰異常源中心。

3.4 不同深度異常源

在上一算例研究的基礎上，將其中一個異常源的質心移到(400 m, 0 m, 350 m)，而其他參數不變，以研究不同深度異常源對估計函數σr影響。

從圖9可以看出，深部異常源受到淺部異常源的干擾，歐拉解有上延的趨勢；未對歐拉解進行過濾時，表層的歐拉解受噪聲干擾比較明顯，淺部的歐拉解的分布很雜亂。

圖10為不同深度異常源歐拉解過濾后分布示意圖。相比于圖1、圖2和圖8，由于深部異常源受到淺部異常源的干擾，構造指數N有逐漸變小的趨勢，因而僅通過構造指數過濾歐拉反褶積繆解，可能會遺漏深部位場總場效應相對較弱的異常源。采用估計函數σr以σ0=0.2作為閾值過濾歐拉反褶積繆解，能有效地標示兩相鄰異常源中心。

圖9 歐拉解分布示意圖(未過濾)Fig.9 Euler deconvolution solutions (before filtered)

圖10 歐拉解分布示意圖(過濾后)Fig.10 Euler deconvolution solutions (after filtered)

4 結論與建議

通過對孤立異常源、同一深度和不同深度異常源基于估計函數σr的歐拉反演解釋，得到如下結論：

1)由于歐拉齊次方程為超定方程，易于受觀測數據所含噪聲水平的影響，且深部異常源受到淺部異常源的干擾，構造指數N有逐漸變小的趨勢，僅通過構造指數過濾歐拉反褶積繆解，可能會遺漏深部位場總場效應相對較弱的異常源。

2)在異常源中心和異常源外部原有估計函數σr的差異性非常小，不利于過濾辨別及過濾歐拉反褶積繆解。

3)改進估計函數σr能有效地標示同一深度和不同深度兩相鄰異常源中心，適應性強。

[1] GEROVSKA D, ARAúZO-BRAVO M. Automatic Interpretation of magnetic data based on euler deconvolution with unprescribed structural index[J]. Computers & Geosciences, 2003,29(8):949-960.

[2] 魯寶亮, 范美寧, 張原慶. 歐拉反褶積中構造指數的計算與優化選取[J]. 地球物理學進展, 2009,24(3): 1027-1031.

LU B L,FAN M N,ZHANG Y Q. The calculation and optimization of structure index in euler deconvolution[J]. Progress in Geophysics, 2009, 24(3): 1027-1031.(In Chinese)

[3] FITZGERALD D,REID A,MCINERNEY P.New discrimination techniques for euler deconvolution[J].Computers and Geosciences,2004,30(5):461-469.

[4] DEWANGAN P,RAMPRASAD T,RAMANA M,et al.Automatic interpretation of magnetic data using euler deconvolution with nonlinear background[J].Pure and Applied Geophysics,2007,164(11):2359-2372.

[5] 曹書錦, 朱自強, 魯光銀. 基于自適應模糊聚類分析的重力張量歐拉解[J]. 中南大學學報(自然科學版), 2012,43(3): 1033-1039.

CAO SH J, ZHU Z Q, LU G Y.Gravity tensor euler deconvolution solutions based on adaptive fuzzy cluster analysis[J]. Journal of Central South University, 2012, 43(3): 1033-1039. (In Chinese)

[6] UGALDE H, MORRIS W. Cluster analysis of euler deconvolution solutions: New filtering techniques and geologic strike determination[J]. Geophysics, 2010,75(3): L61-L70.

[7] QUARTA T. Euler homogeneity equation along ridges for a rapid estimation of potential field source properties[J]. Geophysical Prospecting, 2009,57(4): 527-542.

[8] AGARWAL B, SRIVASTAVA S. Analyses of self-potential anomalies by conventional and extended euler deconvolution techniques[J]. Computers & Geosciences, 2009,35(11): 2231-2238.

[9] STAVREV P,REID A.Degrees of homogeneity of potential fields and structural indices of euler deconvolution[J].Geophysics,2007,72(1):1-12.

[10] HANSEN R, SUCIU L. Multiple-source euler deconvolution[J]. Geophysics, 2002,67: 525-535.

[11] MUSHAYANDEBVU M F, LESUR V, REID A B, et al. Grid euler deconvolution with constraints for 2D structures[J]. Geophysics, 2004,69(2): 489-496.

[12] 姚長利, 管志寧, 吳其斌, 等. 歐拉反演方法分析及實用技術改進[J]. 物探與化探, 2004,28(2): 150-155.

YAO C L, GUAN Z N, WU Q B, et al.An analysis of euler deconvolution and its improvement[J]. Geophysical and Geochemical Exploration, 2004, 28(2): 150-155. (In Chinese)

[13] 范美寧, 孫運生, 田慶君. 關于歐拉反褶積方法計算中的一點改進[J]. 物探化探計算技術, 2005, 27(2): 171-174.

FAN M N, SUN Y S, TIAN Q J. An improvement on calculation of euler decovolution[J]. Computing techniques for geophysical and geochemical exploration, 2005, 27(2): 171-174. (In Chinese)

[14] 史輝, 劉天佑. 利用歐拉反褶積法估計二度磁性體深度與位置[J]. 物探與化探, 2005, 29(3): 230-233.

SHI H,LIU T Y. The application of euler deconvolution to estimating depth and location of the 2-D magnetic body[J].Geophysical and Geochemical Exploration, 2005, 29(3): 230-233. (In Chinese)

[15] 范美寧, 江裕標, 張景仙. 不同數據用于歐拉方程的模型計算[J]. 地球物理學進展, 2008,23(04): 1250-1253.

FAN M N, JIANG Y B, ZHANG J X. Model calculation of euler’s equation for different data types[J]. Progress in Geophysics, 2008,23(4):1250-1253.(In Chinese)

[16] 習宇飛, 劉天佑, 楊坤彪, 等. 歐拉反褶積法用于井中磁測數據反演與解釋[J]. 工程地球物理學報, 2008,5(2): 181-186.

XI Y F,LIU T Y,YANG K B,et al. The application of euler deconvolution to borehole magnetic surveys inversion and interpretation[J]. Chinese Journal of Engineering Geophysics, 2008, 5(2): 181-186.(In Chinese)

[17] 卞光浪, 翟國君, 劉雁春,等. 應用改進歐拉算法解算磁性目標空間位置參數[J]. 測繪學報, 2011,40(03): 386-392.

BIAN G L, ZHAI G J,LIU Y C, et al. Computation of spatial position parameters of magnetic object with improved euler approach[J]. Acta Geodaetica et Cartographica Sinica,2011, 40(3): 386-392.(In Chinese)

[18] BEIKI M.TSVD analysis of euler deconvolution to improve estimating magnetic source parameters: An example from the sele area, sweden[J].Journal of Applied Geophysics,2013,90:82-91.

[19] STAVREV P, GEROVSKA D, ARAUZO-BRAVO M. 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.

[20] THOMPSON D T. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data[J]. Geophysics, 1982, 47(1): 31-37.

Extractionoptimalsolutionofeulerdeconvolutionforgravitydata

ZHOU Yong1,3,4, CAO Shujin2,3, HOU Pingping2, YANG Zilong2

(1.College of Resource Environment and Tourism, Hunan University of Arts and Science, Changde 415000, China;2.Hunan Provincial Key Laboratory of Shale Gas Resource Utilization, School of Resource Environment and Safety Engineering, Hunan University of Science and Technology, Xiangtan 411201, China;3.School of Geosciences and Info-physics, Central South University, Changsha 410083, China;4. Construction and Development of Dongting Lake Ecological Economic Zone of Hunan Synergy Innovation Center, Hunan University of Arts and Science, Changde 415000, China)

It is difficult to characterize multiple anomalies by a single structural index. Euler deconvolution's solutions is hugeous with divergent trend caused by enumerated structural index in traditional Euler decomposition. To overcome those problems, unprescribed Euler deconvolution method is introduced to avoid enumerate structural index, and a novel technology roadmap is proposed for estimating instability of Euler solution based on Euler equation of gravity data solving by truncated singular value decomposition method. A new estimation function is presented based on truncation error to filter the spray solutions. Numerical results show that the technology roadmap of Euler deconvolution is feasible, effective, and has better adaptability in this study.

optimal euler solution; Euler deconvolution; estimate function; truncation error; singular value decomposition

P 631.4

A

10.3969/j.issn.1001-1749.2017.05.03

2016-05-01 改回日期： 2017-06-09

國家自然科學基金(4147114，41374120);湖南省自然科學基金項目(14JJ7038, 2017JJ3069)；湖南科技大學博士科研啟動基金(E51651)；湖南省教育廳科研項目(2010JD42)；湖南科技大學企業技術服務項目(D11646)；湖南科技大學SRIP項目(SYZ2017010)

周勇(1971-)，男，博士，副教授,主要從事地球物理勘探數據處理與反演研究等，E-mail：zhouyong_csu@163.com。

曹書錦(1983-)，男，博士,講師,主要從事淺層地球物理正反演方法研究，Email:shujin.cao@163.com。

1001-1749(2017)05-0598-07