改進(jìn)的局部波數(shù)法及其在磁場數(shù)據(jù)解釋中的應(yīng)用

李麗麗，杜曉娟，馬國慶

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130021

改進(jìn)的局部波數(shù)法及其在磁場數(shù)據(jù)解釋中的應(yīng)用

李麗麗，杜曉娟，馬國慶

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130021

局部波數(shù)法是一種進(jìn)行磁場數(shù)據(jù)解釋的常用方法。現(xiàn)有的局部波數(shù)法在進(jìn)行反演時(shí)往往需要計(jì)算局部波數(shù)的導(dǎo)數(shù)，會(huì)明顯地增大噪聲的干擾，為解釋結(jié)果帶來誤差。提出3種利用磁異常在不同位置或不同高度上局部波數(shù)的簡單組合來進(jìn)行場源體深度及構(gòu)造指數(shù)的計(jì)算方法，不需要計(jì)算局部波數(shù)的導(dǎo)數(shù)，降低了噪聲的干擾，增強(qiáng)了反演結(jié)果的穩(wěn)定性。通過理論模型試驗(yàn)，證明改過的局部波數(shù)法在有無噪聲的情況下均能很好地完成異常的反演，其反演結(jié)果與理論值之間的差距小于理論值的5%。將其應(yīng)用于四川某地區(qū)磁異常的反演中，其結(jié)果與解析信號(hào)的歐拉反褶積法的反演結(jié)果吻合。

磁場；深度；構(gòu)造指數(shù)；導(dǎo)數(shù)；局部波數(shù)

0 前言

磁數(shù)據(jù)解釋的主要任務(wù)是確定場源體的深度及類型（構(gòu)造指數(shù)）。為此人們提出了很多方法來完成這一任務(wù)，有維納反褶積法［1］、解析信號(hào)法［2］、歐拉反褶積法［3－4］、最小二乘法［5－6］、場源參數(shù)成像［7］及人工神經(jīng)網(wǎng)絡(luò)法［8］等，其中場源參數(shù)法是近期應(yīng)用最為廣泛的一種方法。該方法利用局部波數(shù)快速地進(jìn)行場源體反演；但需人為給定研究區(qū)場源體的構(gòu)造指數(shù)，而一個(gè)地區(qū)場源體的構(gòu)造指數(shù)往往是未知的。為了解決這一問題，人們對(duì)該方法進(jìn)行了改進(jìn)［9－11］，使其能夠同時(shí)完成深度及構(gòu)造指數(shù)的反演。然而，改進(jìn)后的方法需要計(jì)算局部波數(shù)的一階甚至更高階導(dǎo)數(shù)，對(duì)數(shù)據(jù)精度要求較高，且會(huì)增大噪聲干擾；當(dāng)測量點(diǎn)距較大時(shí)，高階導(dǎo)數(shù)的計(jì)算非常不穩(wěn)定。Pilkington等［12］在2006年證明局部波數(shù)為解析信號(hào)［13］的導(dǎo)數(shù)，因此局部波數(shù)在二維情況下具有不受磁化方向干擾的優(yōu)點(diǎn)。Salem等［14－15］利用原始局部波數(shù)和相位轉(zhuǎn)換后的局部波數(shù)完成了異常體位置的反演，并于2008年試驗(yàn)了該方法在三維情況下的應(yīng)用效果。Keating［16］利用局部波數(shù)及其垂直導(dǎo)數(shù)對(duì)異常體的深度及構(gòu)造指數(shù)進(jìn)行反演，但是該方法會(huì)增大噪聲的干擾。

筆者對(duì)局部波數(shù)法進(jìn)行改進(jìn)，采用不同位置和不同高度磁異常的局部波數(shù)直接計(jì)算場源體的深度和構(gòu)造指數(shù)，不需要計(jì)算局部波數(shù)的導(dǎo)數(shù)，且降低了噪聲干擾。

1 改進(jìn)的局部波數(shù)法

對(duì)于磁場T而言，其局部波數(shù)kx被定義為局部相位θ的水平方向的導(dǎo)數(shù)［17］：

局部相位θ被定義為

其中，?T／?x和?T／?z分別表示異常T的水平與垂直導(dǎo)數(shù)。對(duì)式（1）進(jìn)行微分后得到

為異常T的解析信號(hào)。Smith等［11］推導(dǎo)出水平位置為x0、埋深為z0的磁性體在點(diǎn)（x，z）處引起的局部波數(shù)kx的表達(dá)式為

其中，n為與場源體類型對(duì)應(yīng)的構(gòu)造指數(shù)。Stavrev等［18］推導(dǎo)了磁源體類型與構(gòu)造指數(shù)之間的對(duì)應(yīng)關(guān)系，其中接觸帶的構(gòu)造指數(shù)為0，巖墻為1，水平圓柱體為2，點(diǎn)源為3。

異常體的水平位置x0可以通過解析信號(hào)和局部波數(shù)的極大值來確定。局部波數(shù)kx在水平位置x＝x0、垂直位置z＝0處的值為

局部波數(shù)kx在水平位置x＝x0＋c，垂直位置z＝0處的值為

將式（6）與式（7）相除后可以得到

利用方程（8）可獲得場源體的深度，然后利用方程（6）計(jì)算場源體的構(gòu)造指數(shù)：

經(jīng)過以上的推導(dǎo)可以看出，利用局部波數(shù)上2個(gè)不同水平位置的點(diǎn)依據(jù)方程（8）和（9）就可以獲得異常體的深度和構(gòu)造指數(shù)。在反演過程中，計(jì)算點(diǎn)應(yīng)選取在場源體中心位置附近；因?yàn)樵谠摲秶鷥?nèi)異常的信噪比較大，得到的反演結(jié)果更準(zhǔn)確。筆者給出了一種利用不同高度的局部波數(shù)kx進(jìn)行磁異常反演的方法，為了使輸出結(jié)果更加穩(wěn)定，采用上延的方法來獲得不同高度上的局部波數(shù)。局部波數(shù)kx在水平位置x＝x0、垂直位置z＝－d處的值為

將（6）式與（10）式相除后得到

分別計(jì)算方程（6）與（10）的倒數(shù)然后相減得到

從上述公式中可以看出，利用2個(gè)不同高度的局部波數(shù)kx，由（11）和（12）式可以得到異常體的深度和構(gòu)造指數(shù)。除上述2種方法外，筆者還提出了另外一種利用相位轉(zhuǎn)換后的局部波數(shù)進(jìn)行磁異常解釋的方法。

Salem等［13］定義相位轉(zhuǎn)換后的局部波數(shù)kz的表達(dá)式：

并推導(dǎo)出水平位置在x0、埋深為z0磁性體的局部波數(shù)kz在點(diǎn)（x，z）處的表達(dá)式為

局部波數(shù)kz在水平位置x＝x0＋b、垂直位置z＝0處的表達(dá)式為

局部波數(shù)kz在水平位置x＝x0＋c、垂直位置z＝0處的表達(dá)式為

將方程（15）與（16）相除后得到

利用（17）式可以得到場源體的深度，然后根據(jù)方程（15）計(jì)算異常體的構(gòu)造指數(shù)：

利用方程（17）和（18）就可以計(jì)算出場源體的深度和構(gòu)造指數(shù)，考慮反演結(jié)果的穩(wěn)定性，b一般選取為1。

經(jīng)過以上的推導(dǎo)，筆者給出進(jìn)行磁性體深度和構(gòu)造指數(shù)計(jì)算的3種方法，這3種方法不需要局部波數(shù)的導(dǎo)數(shù)參與計(jì)算，有效地提高了結(jié)果的穩(wěn)定性。在利用上述方法進(jìn)行反演時(shí)，為了提高結(jié)果的準(zhǔn)確性，計(jì)算點(diǎn)應(yīng)選取在局部波數(shù)kz兩個(gè)極值點(diǎn)之間；因?yàn)閮蓚€(gè)極值點(diǎn)對(duì)應(yīng)于異常體的邊界，當(dāng)測點(diǎn)在異常體上方時(shí)其反演結(jié)果是最準(zhǔn)確的［4］。

2 模型試驗(yàn)

為了試驗(yàn)方法的適用性，設(shè)計(jì)了如下模型：在長為100km的剖面中點(diǎn)位置處存在一個(gè)頂面埋深為4km的巖脈，頂面寬度為10km，磁化強(qiáng)度為20A／km，地磁傾角為60°，磁偏角為0°，測量點(diǎn)距為1 km。利用筆者提出的3種方法對(duì)磁異常進(jìn)行反演，其反演結(jié)果為2個(gè)計(jì)算點(diǎn)之間的水平距離（c）和上延高度（d）的函數(shù)（圖1），并給出了反演結(jié)果與理論值之間的標(biāo)準(zhǔn)偏差（sd）。

從反演結(jié)果中可以看出：當(dāng)水平距離小于14 km、上延高度小于7km時(shí)，本文3種方法反演得到的構(gòu)造指數(shù)均與理論值之間的偏差小于0.01，反演得到的深度與真實(shí)深度之差小于0.1km，其誤差小于理論值的2%，所以本文方法能很好地完成異常的反演，但隨著水平距離與上延高度的增大，標(biāo)準(zhǔn)偏差增大，反演結(jié)果的精度降低。

在進(jìn)行實(shí)際數(shù)據(jù)處理時(shí)，噪聲是必須考慮的一個(gè)因素。為了試驗(yàn)方法存在噪聲時(shí)的應(yīng)用效果，在圖1a的原始異常中加入均值為0、均方差為1nT的隨機(jī)噪聲，采用筆者提出的3種方法對(duì)異常進(jìn)行反演（圖2）。

從圖2中可以看出，在存在噪聲時(shí)本文方法的反演結(jié)果有一定的波動(dòng)。筆者采用如下準(zhǔn)則來獲得最終結(jié)果：1）選取水平距離為0～10km的所有深度解，計(jì)算其平均值；2）將與平均值差異大于平均值25%（根據(jù)實(shí)際情況設(shè)定）的解去除，重新計(jì)算平均值；3）以新平均值為最后深度，計(jì)算構(gòu)造指數(shù)。依據(jù)該準(zhǔn)則獲得的反演結(jié)果如表1所示。

表1 不同方法的反演結(jié)果Table 1 Inversion results of different methods

從表1中可以看出，在存在噪聲的情況下本文方法的反演結(jié)果與理論值之間差距也較小。

為了獲得更加準(zhǔn)確的結(jié)果，可采用向上延拓技術(shù)先對(duì)異常進(jìn)行處理，然后再進(jìn)行反演。筆者對(duì)異常進(jìn)行向上延拓處理，當(dāng)延拓高度為2km時(shí)局部波數(shù)已足夠圓滑，即局部波數(shù)kx為單峰值（圖3）。

從反演結(jié)果（圖3）可以看出，與圖2相比，向上延拓后異常反演結(jié)果的不規(guī)則波動(dòng)消失，計(jì)算反演結(jié)果的平均值，式（8）和（9）反演得到的深度為4.01 km，構(gòu)造指數(shù)為1.01；式（11）和式（12）反演得到的深度為4.02km，構(gòu)造指數(shù)為1.01；式（17）和式（18）反演得到的深度為4.01km，構(gòu)造指數(shù)為1.01，因此對(duì)異常進(jìn)行向上延拓處理后反演結(jié)果的精度得到提高，其值更加準(zhǔn)確。因此，實(shí)際數(shù)據(jù)處理前對(duì)數(shù)據(jù)進(jìn)行向上延拓處理后會(huì)提高反演的精度。

3 實(shí)際數(shù)據(jù)處理

利用筆者提出的方法解釋四川某段實(shí)測航磁異常（圖4）。數(shù)據(jù)采樣間距為100m。將異常向上延拓200m，根據(jù)解析信號(hào)和局部波數(shù)的最大值得出異常體的水平位置在距起始點(diǎn)2 950m處。

從反演結(jié)果可以看出：利用不同方法反演得到的構(gòu)造指數(shù)均接近1，表明地質(zhì)體為一巖墻，反演得到深度的平均值為543m。

圖1 模型理論異常及不同方法的反演結(jié)果Fig.1 Theoretical anomaly and inversion results of different methods

該地區(qū)沒有其他的數(shù)據(jù)來驗(yàn)證反演結(jié)果的正確性，將筆者的反演結(jié)果與解析信號(hào)歐拉反褶積法［19］的反演結(jié)果進(jìn)行對(duì)比。解析信號(hào)歐拉反褶積法的反演結(jié)果顯示異常體的位置在（2 945±7）m，深度為（554.6±18.2）m，構(gòu)造指數(shù)大小為0.98±0.03。這與本文方法的反演結(jié)果相一致，表明本文方法反演結(jié)果是準(zhǔn)確的。

圖2 加噪后異常及不同方法的反演結(jié)果Fig.2 Noise－corrupted anomaly and the inversion results of different methods

圖3 濾波后異常及不同方法的反演結(jié)果Fig.3 Smoothed anomaly and the inversion results of the methods

圖4 實(shí)際航磁異常的反演結(jié)果Fig.4 Inversion results of real aeromagnetic anomaly

4 結(jié)論

1）筆者提出3種改進(jìn)的局部波數(shù)法進(jìn)行磁性體反演方法，這些方法僅利用不同水平位置和不同高度的局部波數(shù)來進(jìn)行深度和構(gòu)造指數(shù)的計(jì)算，不需要計(jì)算局部波數(shù)的導(dǎo)數(shù)，降低了對(duì)數(shù)據(jù)精度的要求。

2）通過理論模型試驗(yàn)證明，新方法在有無噪聲情況下，均能準(zhǔn)確地計(jì)算出地質(zhì)體的深度和構(gòu)造指數(shù)。

3）在有噪聲情況下，采用向上延拓技術(shù)處理數(shù)據(jù)后可獲得更加精確的結(jié)果，為實(shí)際數(shù)據(jù)解釋提供了理論依據(jù)。

4）將新方法應(yīng)用于四川實(shí)測剖面磁數(shù)據(jù)的解釋中，其反演結(jié)果與解析信號(hào)的歐拉反褶積法的結(jié)果相吻合，表明本文方法在實(shí)際中也有較好的應(yīng)用效果。

（References）：

［1］ Hartman R R，Teskey D J，F(xiàn)riedberg J L.A System for Rapid Digital Aeromagnetic Interpretation［J］.Geophysics，1971，36：891－918.

［2］ Nabighian M N.The Analytic Signal of Two Dimensional Magnetic Bodies with Polygonal Cross Section：Its Properties and Use for Automated Anomaly Interpretation［J］.Geophysics，1972，37：507－517.

［3］ Thompson D T.“EULDPH”：A New Technique for Making Computer－Assisted Depth Estimates fromMagnetic Data［J］.Geophysics，1982，47：31－37.

［4］ Reid A B，Allsop J M，Granser H，et al.Magnetic Interpretation in Three Dimensions Using Euler Deconvolution［J］.Geophysics，1990，55：80－91.

［5］ Abdelrahman E M.Discussion on“A Least－Squares Approach to Depth Determination from Gravity Data”by O P Gupta［J］.Geophysics，1990，55：376－378.

［6］ Salem A，Ravat D，Mushayandebvu M F，et al.Linearized Least－Squares Method for Interpretation of Potential Field Data from Sources of Simple Geometry［J］.Geophysics，2004，69：783－788.

［7］ Thurston J B，Smith R S.Automatic Conversion of Magnetic Data to Depth，Dip，and Susceptibility Contrast Using the SPI Method［J］.Geophysics，1997，62：807－813.

［8］ Elawadi E，Salem A，Ushijima K.Detection of Cavities from Gravity Data Using A Neural Network［J］.Exploration Geophysics，2001，32：75－79.

［9］ Thurston J B，Smith R S，Guillon J.A Multi－Model Method for Depth Estimation from Magnetic Data［J］.Geophysics，2002，67：555－561.

［10］ Salem A，Smith R S.Depth and Structural Index from the Normalized Local Wavenumber of 2DMagnetic Anomalies［J］.Geophysical Prospecting，2005，51：83－89.

［11］ Smith R S，Thurston J B，Dai T，et al.ISPI：The Improved Source Parameter Imaging Method［J］.Geophysical Prospecting，1998，46：141－151.

［12］ Pilkington M，Keating P.The Relationship Between Local Wavenumber and Analytic Signal in Magnetic Interpretation［J］.Geophysics，2006，71（1）：L1－L3.

［13］ 馬國慶，杜曉娟，李麗麗.利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行位場數(shù)據(jù)的邊界識(shí)別［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41（增刊）：345－348.

Ma Guo－qing，Du Xiao－juan，Li Li－li.Edge Detection of Potential Field Using Correlation Coefficient of Horizontal and Vertical Derivatives［J］.Journal of Jilin University：Earth Science Edition，2011，41（Sup.）：345－348.

［14］ Salem A，Ravat D，Smith R，et al.Interpretation of Magnetic Data Using an Enhanced Local Wavenumber（ELW）Method［J］.Geophysics，2005，70（2）：L7－L12.

［15］ Salem A，Williams S，F(xiàn)airhead D.Interpretation of Magnetic Data Using Tilt－Angle Derivatives［J］.Geophysics，2008，73：L1－L10.

［16］ Keating P.Improved Use of the Local Wavenumber in Potential－Field Interpretation［J］.Geophysics，2009，74（6）：L75－L85.

［17］ Bracewell R N.The Fourier Transform and Its Application［M］.New York：McGraw Hill Book Co，1965.

［18］ Stavrev P，Reid A.Degrees of Homogeneity of Potential Fields and Structural Indicates of Euler Deconvolution［J］.Geophysics，2007，72（1）：L1－L12.

［19］ Keating P，Pilkington M.Euler Deconvolution of the Analytic Signal and Its Application to Magnetic Interpretation［J］.Geophysical Prospecting，2004，52：165－182.

Improved Local Wavenumber Methods in the Interpretation of Magnetic Fields

Li Li－li，Du Xiao－juan，Ma Guo－qing

College of GeoExploration Science and Technology，Jilin University，Changchun 130026，China

Local wavenumber method is a widely used method in the interpretation of magnetic fields，the previous local wavenumber methods need to compute the derivatives of local wavenumber in the inversion，which increase the effect of noise obviously，and bring the errors to inversion results.The authors proposed using the simple combinations of different locations and heights of local wavenumber of magnetic anomaly to compute the depth and structural index of the causative source，and these methods do not need to compute the derivatives of the local wavenumber，which reduce the effect of the noise，and increase the stability of the inversion results.The new methods are demonstrated on theoretical magnetic field data without and with random noise，and the difference between the inversion results and theoretical values is less than 5%of the theoretical values.The practical application of the new methods is tested on real magnetic anomaly from Sichuan basin，Southwest China，the inversion results are consistent with the results computed by the Euler deconvolution of analytic signal.The proposed local wavenumber methods can finish the inversion of magnetic data effectively.

magnetic fields；depth；structural index；derivatives；local wavenumber

book=2012,ebook=558

P631.2

A

1671－5888（2012） 04－1179－07

2011－10－23

國家深部探測技術(shù)與實(shí)驗(yàn)研究專項(xiàng)項(xiàng)目（SinoProbe－09－01）

李麗麗（1983－），女，博士研究生，主要從事地球物理數(shù)據(jù)處理與解釋方面的研究，E－mail：lilili10＠m(xù)ails.jlu.edu.cn。