重力場歐拉反褶積最優(yōu)解提取

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

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

重力場歐拉反褶積最優(yōu)解提取

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

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

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

歐拉反褶積； 最優(yōu)解； 估計函數(shù)； 截斷誤差； 奇異值分解

0 引言

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

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

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

1 歐拉反褶積反演

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

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

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

Gm=d

(5)

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

為評價歐拉解穩(wěn)定性，特設(shè)定一個估計函數(shù)以過濾歐拉解解集中的繆解[19-20]

(10)

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

(11)

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

3 模型試驗

3.1 算法驗證

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

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

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

3.2 孤立異常源

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

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

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

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

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

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

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

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

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

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

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

3.3 同一深度異常源

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

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

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

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

3.4 不同深度異常源

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

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

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

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

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

4 結(jié)論與建議

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

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

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

3)改進估計函數(shù)σ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] 魯寶亮, 范美寧, 張原慶. 歐拉反褶積中構(gòu)造指數(shù)的計算與優(yōu)化選取[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] 姚長利, 管志寧, 吳其斌, 等. 歐拉反演方法分析及實用技術(shù)改進[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] 范美寧, 孫運生, 田慶君. 關(guān)于歐拉反褶積方法計算中的一點改進[J]. 物探化探計算技術(shù), 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] 范美寧, 江裕標, 張景仙. 不同數(shù)據(jù)用于歐拉方程的模型計算[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] 習宇飛, 劉天佑, 楊坤彪, 等. 歐拉反褶積法用于井中磁測數(shù)據(jù)反演與解釋[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] 卞光浪, 翟國君, 劉雁春,等. 應用改進歐拉算法解算磁性目標空間位置參數(shù)[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)；湖南科技大學企業(yè)技術(shù)服務(wù)項目(D11646)；湖南科技大學SRIP項目(SYZ2017010)

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

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

1001-1749(2017)05-0598-07