馬國慶，杜曉娟，李麗麗，李文成，王鳳剛

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

2.首鋼地質(zhì)勘查院地質(zhì)研究所，北京 100144

0 引言

重磁異常的梯度信息能更好地描述地下地質(zhì)體的細(xì)節(jié)特征，在近幾年應(yīng)用十分廣泛，國內(nèi)外已經(jīng)陸續(xù)開展重磁梯度測量。歐拉反褶積法是一種常用的位場自動(dòng)解釋方法［1－2］，人們通過對常規(guī)歐拉反褶積公式求導(dǎo)實(shí)現(xiàn)了梯度數(shù)據(jù)的解釋［3－5］，獲得了較好的實(shí)際應(yīng)用效果，但是該方法仍需要構(gòu)造指數(shù)參與計(jì)算。由于在實(shí)際數(shù)據(jù)解釋中測區(qū)內(nèi)地質(zhì)體的類型是未知的，構(gòu)造指數(shù)無法準(zhǔn)確地獲得，且構(gòu)造指數(shù)的微小誤差會(huì)使反演結(jié)果產(chǎn)生較大的影響［6－8］，筆者提出梯度反褶積法完成梯度數(shù)據(jù)的反演。

梯度反褶積法是在常規(guī)歐拉反褶積法的基礎(chǔ)上推導(dǎo)出來的，且在公式中消去了構(gòu)造指數(shù)，避免了因構(gòu)造指數(shù)選取不當(dāng)造成的誤差，使反演結(jié)果更加準(zhǔn)確。通過理論模型證明梯度反褶積法在位場梯度反演中的有效性，且相對常規(guī)歐拉反褶積法其精度較高，適用性更強(qiáng)。將本方法應(yīng)用于某鐵礦區(qū)的實(shí)測磁梯度數(shù)據(jù)，確定了鐵礦的賦存狀態(tài)。

1 梯度反褶積法

梯度反褶積法以梯度數(shù)據(jù)為基礎(chǔ)進(jìn)行場源體的反演，是在歐拉反褶積法基礎(chǔ)上推導(dǎo)出來的。常規(guī)歐拉反褶積法［9－10］的公式為

其中：T表示重力異常或磁異常；（x0，z0）分別表示異常體的埋藏位置；N為地質(zhì)體的構(gòu)造指數(shù)；B表示背景異常場。為了使式（1）可適用于梯度數(shù)據(jù)，分別計(jì)算其在x，z方向的導(dǎo)數(shù)：

對比發(fā)現(xiàn)，式（4）與式（5）右側(cè)相等，因此左側(cè)也相等，因此可以得到

利用式（7）通過重磁梯度異常可以獲得異常體的位置信息（x0，z0）。

從式（7）中可以看出，梯度反褶積公式完成了梯度數(shù)據(jù)的反演，并且去除了構(gòu)造指數(shù)的影響。

同理，可以推導(dǎo)出三維情況下梯度反褶積法的表達(dá)式：

2 理論模型試驗(yàn)

直角坐標(biāo)系下，存在水平位置為（55m，55m）、埋深為50m、半徑為7m的球體，當(dāng)計(jì)算點(diǎn)距為1 m時(shí)引起的梯度異常如圖1所示。為了試驗(yàn)梯度反褶積法的應(yīng)用效果，分別利用常規(guī)歐拉反褶積法和梯度反褶積法對梯度異常進(jìn)行反演，反演結(jié)果如圖2所示。

從圖2a、b中可以看出，常規(guī)歐拉反褶積法和梯度歐拉反褶積法均能很好地給出異常體的真實(shí)水平位置。圖2c和圖2d分別為利用這2種方法計(jì)算得到的深度統(tǒng)計(jì)結(jié)果，可以看出，2種方法均能很準(zhǔn)確地完成異常體深度的反演，但梯度歐拉反褶積法得到的結(jié)果發(fā)散程度較小，且誤差小于5%。

為了試驗(yàn)方法的穩(wěn)定性，將計(jì)算點(diǎn)距加大到5 m。因?yàn)榍蝮w的半徑為7m，因此在其正上方僅存在4個(gè)點(diǎn)。在點(diǎn)數(shù)比較少的情況下反演結(jié)果的準(zhǔn)確度會(huì)降低。分別利用常規(guī)歐拉反褶積法和梯度反褶積法對球體引起的異常進(jìn)行反演，反演結(jié)果如圖3所示。從圖3中可以看出：常規(guī)歐拉反褶積法在計(jì)算點(diǎn)距較大的情況下，反演得到的異常體的水平位置分散，無法準(zhǔn)確給出場源體的位置信息，其結(jié)果與理論值之間偏差較大；梯度反褶積法仍能較好完成異常的反演，能很好地描述地質(zhì)體的真實(shí)位置。

圖1 球體引起的重力異常及梯度異常Fig.1 Gravity and gradient anomaly of sphere

圖2 球體重力異常反演結(jié)果Fig.2 Inversion results of gravity anomaly caused by sphere

圖3 球體重力異常反演結(jié)果Fig.3 Inversion results of gravity anomaly caused by sphere

下面試驗(yàn)一下2種方法存在多個(gè)異常體時(shí)的情況。在水平位置（35m，35m）和（65m，65m）處存在埋深均為15m的球體，并在引起的重力異常中加入均值為0、方差為1mGal的高斯噪聲。其原始重力異常及其梯度異常如圖4所示。

分別利用常規(guī)歐拉反褶積法和梯度歐拉反褶積法對圖4中所示異常進(jìn)行反演，結(jié)果如圖5所示。

圖4 球體引起的重力異常及梯度異常Fig.4 Gravity and gradient anomaly of sphere

圖5 球體重力異常的反演結(jié)果Fig.5 Inversion results of gravity anomaly caused by sphere

從圖5的反演結(jié)果中可以看出：由于噪聲的影響，常規(guī)歐拉反褶積法的反演結(jié)果在異常周圍產(chǎn)生了一定的干擾，即產(chǎn)生多余的干擾異常，為解釋帶來困難；梯度反褶積法的反演結(jié)果受噪聲干擾較小，所得到的結(jié)果能較好地反映異常體的位置和深度。為了試驗(yàn)梯度反褶積法對不同形狀地質(zhì)體的應(yīng)用效果，采用了上頂埋深為10m、邊長為20m的正方體，其引起的原始重力及其梯度異常如圖6所示。

分別利用常規(guī)歐拉反褶積法和梯度反褶積法對正方體進(jìn)行反演，在利用常規(guī)歐拉反褶積進(jìn)行反演時(shí)構(gòu)造指數(shù)選取0，2種方法的反演結(jié)果如圖7所示。

圖6 正方體引起的重力異常及梯度異常Fig.6 Gravity and gradient anomaly caused by square

圖7 正方體引起異常的反演結(jié)果Fig.7 Inversion results of gravity anomaly caused by square

從圖7可以看出：常規(guī)歐拉反褶積法難以確定異常體的水平位置，且深度與理論值之間也存在一定的誤差；梯度反褶積法能很好地圈定異常體的范圍，且反演得到的深度與理論值一致，其均值與理論值之間的誤差小于5%。

通過上述試驗(yàn)可以看出，梯度反褶積法反演結(jié)果的精度明顯優(yōu)于常規(guī)歐拉反褶積法的反演結(jié)果，能更加有效地完成位場異常的反演，且適用性較強(qiáng)。

3 實(shí)際數(shù)據(jù)應(yīng)用

應(yīng)用梯度反褶積法解釋實(shí)際航磁梯度數(shù)據(jù)。圖8為河北遷安地區(qū)實(shí)測航磁及梯度異常。遷安地區(qū)主要發(fā)育沉積變質(zhì)型鐵礦，經(jīng)過50多年的開采，地形地貌遭到嚴(yán)重破壞，磁異常形態(tài)已經(jīng)發(fā)生很大改變。磁梯度測量具有不受外界環(huán)境及日變的干擾，數(shù)據(jù)更加穩(wěn)定，因此采用磁梯度異常來確定地下鐵礦的分布特征，為下一步開采提供基礎(chǔ)性依據(jù)。本次航磁梯度測量的總誤差為0.002 1nT／m，分別利用常規(guī)歐拉反褶積法和梯度反褶積法對該地區(qū)數(shù)據(jù)進(jìn)行反演。根據(jù)異常的形態(tài)特征，在利用常規(guī)歐拉反褶積進(jìn)行反演時(shí)構(gòu)造指數(shù)選取1，反演結(jié)果如圖9所示。

從圖9可以看出：常規(guī)歐拉反褶積法的結(jié)果較發(fā)散，對于詳查測線的布設(shè)及井口的布設(shè)不能提供強(qiáng)有力的支持；梯度反褶積法計(jì)算結(jié)果明顯優(yōu)于常規(guī)歐拉反褶積結(jié)果，該方法的結(jié)果比較集中，離散點(diǎn)比較少，能較準(zhǔn)確地給出鐵礦的分布范圍，為開采計(jì)劃提供明確的方向。

4 結(jié)論

圖8 鐵礦區(qū)航磁異常及梯度異常Fig.8 Aeromagnetic and gradient anomaly of an iron ore area

1）筆者提出的位場反演方法——梯度反褶積法，以歐拉反褶積法為基礎(chǔ)，通過求導(dǎo)被推廣到梯度數(shù)據(jù)處理中，并在反演公式中消去了構(gòu)造指數(shù)，避免了構(gòu)造指數(shù)選取不當(dāng)給反演結(jié)果帶來的誤差 。

2）通過模型試驗(yàn)，驗(yàn)證了梯度反褶積法在位場反演中的有效性，其反演精度高于常規(guī)歐拉反褶積法。

3）梯度反褶積法對鐵礦區(qū)實(shí)測異常的反演結(jié)果集中、離散點(diǎn)少，能較準(zhǔn)確地給出鐵礦的分布范圍，為下一步開采計(jì)劃提供基礎(chǔ)性依據(jù)。

圖9 鐵礦區(qū)磁異常反演結(jié)果Fig.9 Inversion results of aeromagnetic anomaly in the iron ore area

黃大年教授在論文修改過程中提出了寶貴意見，工作組的同仁們?yōu)楸疚奶峁┝藢氋F的資料，在此一并致謝。

（References）：

［1］馬國慶，孟令順，李麗麗.龍門山及鄰區(qū)斷裂分布及地震前后斷裂形態(tài)差異［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2012，42（2）：519－525.Ma Guoqing，Meng Lingshun，Li Lili.Fault Distribution of Longmenshan and Adjacent Regions and Fault Morphological Differences Before and After Earthquake［J］.Journal of Jilin University：Earth Science E－dition，2012，42（2）：519－525.

［2］李麗麗，杜曉娟，馬國慶.改進(jìn)的局部波數(shù)法及其在磁場數(shù)據(jù)解釋中的應(yīng)用［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2012，42（4）：1179－1185.Li Lili，Du Xiaojuan，Ma Guoqing.Improved Local Wavenumber Methods in the Interpretation of Magnetic Fields［J］.Journal of Jilin University：Earth Science Edition，2012，42（4）：1179－1185.

［3］Huang D，Gubbins D，Clark R A，et al.Combined Study of Euler’s Homogeneity Equation for Gravity and Magnetic Field［C］／／57th EAGE Conference.Glasgow：European Association of Geoscientists and Engineers，1995：144.

［4］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（3）：783－788.

［5］Schimidt P W，Clark A.The Magnetic Gradient Tensor：Its Properties and Uses in Source Characterization［J］.The Leading Edge，2006，25（1）：75－78.

［6］Barbosa V C F，Silva J B C，Medeiros W E.Stability Analysis and Improvement of Structural Index Estimation in Euler Deconvolution［J］.Geophysics，1999，64（1）：48－60.

［7］Daniela G，Marcos J，Arauzo B.Automatic Interpretation of Magnetic Data Based on Euler Deconvolution with Unprescribed Structural Index［J］.Computers ＆Geosciences，2003，29（8）：949－960.

［8］Hansen R O，Laura S.Multiple－Source Euler Deconvolution［J］.Geophysics，2002，67（2）：525－535.

［9］Thompson D T.‘EULDPH’：A New Technique for Making Computer－Assisted Depth Estimates from Magnetic Data［J］.Geophysics，1982，47（1）：31－37.

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