张明华
（北京科技大学资源工程学院，北京 100083）
管志宁
（中国地质大学物探系，北京 100083）
摘要 在分析以往重磁界面反演方法的基础上，本文提出了拟合地壳内部磁化强度分布的磁化强度函数和进一步的反演居里深度和磁性地层界面深度的正反演迭代算法。该方法克服了以往方法的不足，节省时间，而且具有很高的反演精度和准确度。在以寻找大型和超大型矿床为目标的国家“攀登计划-B”项目的课题研究工作中，于华北地台北缘地区获得了对地质研究和贵金属找矿都有价值的结果。
关键词 磁化强度函数 高精度 地层界面
1 磁化强度函数和磁性层界面模型
关于地壳区域磁场场源的形成机制和分布范围，不同的研究工作者提出了不同的假设。目前的知识来自对岩石圈和大陆科学钻探的研究。根据目前所得到的地壳综合深部岩石物理剖面参数和磁性矿物的居里点数据[1,2,3,6,8]，我们认为，岩石磁性是随深度变化的，在温度高于磁性矿物居里点的深度（简单称为居里深度）以下变为顺磁性。同时，岩石磁化强度在横向上随岩性和地质构造单元的不同而变化。这就是地壳磁化强度分布的基本特征。我们提出一个能够较好地拟合这一地壳岩石磁性特征的磁化强度函数，见图1a。表达式为：
第30届国际地质大会论文集 第20卷 地球物理
其中，J（ξ，η，ζ）是磁化强度函数，n=1或2,a（ξ，η）、b（ξ，η）＞0和c（ξ，η）≥0是在横向上随不同磁性地质构造单元而变化的变量，ζ是垂向上的深度变量。
我们采用如图1b所示的磁性地层模型用于区域磁异常解释。这是一个上、下界面均为起伏的磁性层模型。模型的上界面以上是弱磁性或无磁性的沉积岩等介质，下界面（居里深度）以下的岩石变为顺磁性。磁性层内存在着由磁化强度函数所描述的纵、横向上的磁性差异。
2 反演磁性界面深度的迭代法
图1 磁化强度函数与磁性地层模型示意图
a—地壳结构、磁化强度变化及磁化强度函数；b—磁性层模型
为了利用航磁资料对深部地质构造进行解释、了解居里深度和太古宇顶面深度的变化，结合国家“攀登计划-B”项目课题的工作，我们研究了以往所有这方面的反演方法。传统方法或者由于其对地壳磁性假设的近似程度较低，或者由于在反演过程中对磁异常或深度值施加滤波[2,9,4]，反演结果准确性和精度受到限制。为此，我们基于前述的磁化强度函数和磁性层模型，研究并提出一种正反演结合迭代反演磁性界面深度的方法，英文简写为MIDI。为节约篇幅，这里给出取n=1和c（ξ，η）＝0时的正演计算公式。
利用平面磁场垂直分量Z（x,y）的化极场频谱[2]，
第30届国际地质大会论文集 第20卷 地球物理
其中，Z（）uv和F{}表示傅里叶变换或频谱。J（ξ，η，ζ）为磁性层磁化强度。h（ζ，η）=h＋Δh（ζ，η），H（ζ，η）＝H＋△H（ζ，η），分别为上下界面深度。h,H为上下界面的平均深度。Δh（ζ，η）和△H（ζ，η）分别为上下界面的起伏幅度。R（x—ξ，y—η，ζ）＝[（x—ζ）2＋（y—η）2＋ζ2]-1/2。F{R}=exp（-2πf）exp[-2πj（uξ＋vη）]/f。
利用指数函数的泰勒级数展开，有：
第30届国际地质大会论文集 第20卷 地球物理
我们将磁化强度表达式（1）代入（2）式，最终导出磁性按磁化强度函数变化、上下界面都起伏的磁性层磁场频谱表达式为：
第30届国际地质大会论文集 第20卷 地球物理
磁场垂直分量Z（x,y）可以由实测得到，或由ΔT（x,y）转换得到。
由正演公式，若其中一个界面的深度已知，则可反演另一界面深度变化。不计级数项的部分可用于直接反演，称为直接反演公式。由直接反演公式结合级数项，即可构筑迭代反演。迭代计算在第n次反演之界面深度值与第n—1次反演之深度值之差，或相邻两次中间磁场频谱值之差满足精度要求时，迭代终止。
3 计算方法的精度
常磁化强度情况是本文算法（MIDI）的一个特例。此情况下，MIDI方法与已有方法结果一致。在变磁化强度情况下，理论模型上的反演计算表明，无论对于连续起伏界面，还是对于跳跃起伏界面（尤其是大断裂所致的界面起伏），也不论起伏的形式和幅度如何，MIDI反演都很好地，可以说准确而且相对省时地收敛到已知的理论真值，只是要求反演结果的精确愈高，计算迭代的时间愈长而已。图2a是一个迭代精度要求达到5%即终止的例子。在4M内存的486微机上迭代355次，花费22min。由图明显可见，反演深度等值线（虚线）与理论模型深度（实线）几乎重合。
MIDI方法的一个特点是，不需要对面积性场值进行滤波处理，也不对反演过程中的深度或场值作滤波等可能改变场源特征的处理。所以，反演计算速度快，而且对深度的分辨能力强。使用图2b所示的起伏界面，MIDI方法与Parker法[4,10]进行了对比。用同样的磁化强度计算出磁场值，再由磁场反演界面深度起伏，反演结果分别见图2c、2d。二者的计算时间都限制为20min。显然，MIDI方法的深度分辨能力较强。
4 关于平均深度的讨论
平均深度的选择一直是重磁界面反演的一个难题。在MIDI算法，当一个界面已知时，由于磁化强度的非线性变化是事先约束给定的，所以另一个要反演界面的平均深度，可以由反演计算来调整确定。给出一个模型计算的例子如图3。给定上界面的平均深度值h=2.0km，其深度起伏为图2b所示的模型。下界面的平均深度H=5.0km，下界面水平。计算这一磁性层模型的磁场Z（x,y）。然后，下界面不变，反演上界面。当h给以不同的值时，反演出的深度起伏幅度值Δh（x,y）的平均值（Ria）不同。给定的值大于真值2.0时，Ria为负；给定的值小于真值2.0时，Ria为正。给定的h值越接近真值，Ria越小。当h取真值时，Ria趋于零（由于计算误差而不为零）。Ria对h的真值具有明显的指向性。因此，我们可用Ria做为指示因子来求解平均深度。
图2 理论模型上的方法精度检验
a—实线是给定模型的深度等值线，虚线是MIDI反演结果等值线，MIDI方法与Parker方法对比结果。实例2:b—界面深度模型；c—Parker方法反演结果；d—在与Parker方法计算时间相同时的MIDI反演结果。深度单位：km
图3 界面平均深度的调整计算
当平均深度接近真值时，深度起伏平均值（Ria）趋于0
5 呼和浩特—张家口地区居里深度反演
地壳磁性层下界深度即居里深度，是由地壳温度场决定的，是地下热状态的一个重要指标。它对深部地质构造研究、地震学研究和矿产预测都有重要意义。我们选择深部构造及上部矿产都具有重要价值的呼和浩特—张家口地区，53200km2范围，进行了居里深度分布的反演研究。图4a是该区太古宙地层、岩浆岩、侵入玄武岩分布图。太古宇和玄武岩具有较强磁性。图4b是该区航磁异常图。
图4 呼和浩特—张家口地区太古宙地层分布图（a）与航磁异常图（b）
异常值单位：nT
（1）综合该地区岩石及地层磁性统计结果及有关的已知地层深度结果，大致得磁化强度于不同深度变化的数据，由遗传算法[5]选择出的磁化强度函数为：
第30届国际地质大会论文集 第20卷 地球物理
（2）为提取居里深度的磁场信息，采用正则化滤波（因子选L=40Km）方法消除浅层磁性变化影响。滤波结果作为居里深度界面起伏和磁性基底顶界起伏二者的综合异常。
（3）根据已知的地层深度资料和人工地震剖面解释[7,11,13]，由已知的太古宙地层（Ar）顶界面的位置，通过插值得到全区磁性基底顶界起伏的大致深度。其平均深度在114°以西为3.2km,114°以东为3.5km。这一深度分布作为磁性层上界面。
（4）通过计算将居里深度调整为33.0km。按反演迭代中间结果中，前后两次的深度值相差5%为迭代终止条件，所得反演结果见图5。可以看到，居里深度的主要凸起与本区主要的深断裂构造相一致。
图5 呼和浩特—张家口地区居里面深度变化图
深度等值线单位：km
6 蔡家营地区太古宇顶面反演
华北地台北缘地区金、铅、锌等多金属矿产一般与太古宙地层（有时与元古宙地层）顶面起伏和断裂构造及岩浆活动（尤其燕山期侵入岩）关系密切。太古宙地层巨厚，且与上覆地层相比具有强磁性。因而利用航磁异常研究隐伏构造和太古宙地层的顶界面具有较好物理前提。我们选择正在进一步深入和扩大找矿的张家口市蔡家营地区，26112km2范围，进行了研究。对1/20万航磁异常进行消去浅地表干扰的滤波处理[12]，利用前述的居里深度结果参与计算，我们得到的蔡家营地区太古宇顶面深度起伏结果见图6。其深度变化与一些剖面上布格重力异常遗传算法反演的结果[5]几乎完全一样。结合该区地质构造、岩浆岩分布和重力异常，我们对反演结果解释于下。
（1）在蔡家营、土城子的南面都存在老地层隆起。蔡家营矿位于老地层四周隆起的凹陷之中，目前的蔡家营矿位于凹陷的次级平缓地段。这一凹陷对应着重力相对低异常，由于该区燕山期岩浆岩密度较低，此凹陷上应有一燕山期隐伏岩体存在。地表已见岩株。与此类似，土城子南面的隆起南侧，张北的西面有一类似的平缓地段。这两个地段应是寻找金、铅、锌等多金属矿床的有望地段。
（2）尚义—张北—土城子之间及其西北的新生代玄武岩层应是比较薄的。从航磁图上以及滤波处理结果上可以清楚说明这一点。土城子西南的太古宇隆起应是在该层玄武岩之下。
（3）尚义到康保之间应该存在一条深大断裂（图6）。这一点从航磁异常图及反演的太古宇顶界起伏可以清楚地看到。断裂两侧的磁场特征和太古宇深度截然不同。
致谢 研究工作得到中国科学院院士刘光鼎先生支持与指导，刘士毅教授级高工对蔡家营地区地球物理解释给予极大帮助，谭承泽教授和中国地质大学古地磁实验室其他老师、姚长利讲师、中科院地球物理研究所郝天珧副研究员为研究工作提供很大帮助，在此一并衷心感谢。
图6 菜家营地区太古宇顶面起伏形态图
虚线是推测的深断裂，深度数据单位：km
参考文献
[1] Kelochihovskarya，3．A．.Magnetic nonhomogeneity of the crust.Journal of Geophysics（in Russian）,1986,8（5）:3～23.
[2] 管志宁，安玉林．区域磁异常定量解释．北京：地质出版社，1991.
[3] 郭武林．磁异常区域分量的岩石学解释．物化探译丛，1987,51（4）:26～32.
[4] Parker,R.L..The rapid calculation of potential anomalies.Geophysics,J.R.Soc.,1973,31:445～447.
[5] 张明华．重磁数据处理新方法研究及其在华北地台北缘地区应用[博士学位论文]．中国地质大学研究生院（北京），1995.
[6] Hanson,R.O..Nationwide Curie point depth analysis of Japan.53th annual meeting,SEG,1983.
[7] 地质矿产部华北地台北缘物化探遥感编图委员会航磁、重力解释组．华北地台北缘航磁、重力解释成果报告（1:1000000）.1993,12.
[8] Hanson,et al..Magnetic anomalies in Japan and its adjacent regions.Special issue,Journal of eomagnetism and Geo-electricity,1994,46（6）.
[9] 穆石敏，申宁华，孙运生．区域地球物理数据处理方法及其应用．吉林：吉林科学技术出版社，1990.
[10] Parker,R.L.and Huestis,S.P..The inversion of magnetic anomalies in the presence of topography.Journal of Geophysical Research,1974,79（11）:1587～1592.
[11] 王宫福，张文治．华北陆台北缘重磁异常的解释．物探与化探，1992,16（1）:31～37.
[12] 王继伦．二维最佳线性数字滤波器的设计原理．地球物理学报，1977,20（2）:159～172.
[13] 余钦范，马杏垣．华北地区航磁图像处理结果和地震构造解释．地震地质，1989,11（4）:5～14．