居里深度和地层界面深度的高精度反演
(北京科技大学资源工程学院，北京100083)
关志宁
(中国地质大学地球物理勘探系，北京100083)
本文在分析以往重磁界面反演方法的基础上，提出了一种拟合地壳磁化强度分布的磁化强度函数，并进一步提出了反演居里深度和磁性地层界面深度的正反向迭代算法。该方法克服了以往方法的缺点，节省时间，具有较高的反演精度和准确度。在寻找大型和超大型矿床的国家“攀登计划-B”项目的研究工作中，在华北地台北缘地质研究和贵金属找矿中取得了有价值的成果。
磁化函数；高精度地层界面
1磁化函数和磁性层界面模型
关于地壳区域磁场源的形成机制和分布范围，不同的研究者提出了不同的假说。目前的知识来源于岩石圈研究和大陆科学钻探。根据综合深部地壳的物理剖面参数和磁性矿物的居里点数据，我们认为岩石磁性随深度而变化，当温度高于磁性矿物的居里点(简称居里深度)时变为顺磁性。同时，岩石的磁化强度在横向上随岩性和地质构造单元而变化。这是地壳磁化分布的基本特征。我们提出一个磁化函数，可以很好地拟合这种地壳岩石的磁性特征，如图1a所示。表达式是:
第30届国际地质大会会议录第20卷地球物理学
其中J(ξ，η，ζ)是磁化强度的函数，n=1或2，a(ξ，η)，b(ξ，η) > 0，c(ξ，η)≥0是横向上随不同磁性地质构造单元而变化的变量。
我们使用如图1b所示的磁性地层模型来解释区域磁性异常。这是一个上下界面起伏的磁性层模型。模型上界面以上的介质为弱磁性或无磁性沉积岩，下界面(居里深度)以下的岩石变为顺磁性。磁性层中存在由磁化函数描述的垂直和水平磁差。
2反演磁性界面深度的迭代法
图1磁化函数和磁性地层模型示意图
a——地壳结构、磁化变化和磁化函数；b磁层模型
为了利用航磁资料解释深部地质结构，了解居里深度和太古宙顶深的变化，我们结合国家“攀登计划-B”项目的工作，研究了以往所有的反演方法。由于传统方法对地壳磁假设的近似程度低，或者在反演过程中对磁异常或深度值进行了滤波，从而限制了反演结果的准确性和精度。因此，基于上述磁化函数和磁层模型，我们研究并提出了一种正反迭代相结合反演磁界面深度的方法，英文缩写为MIDI。为了节省篇幅，这里给出了n=1，c(ξ，η) = 0时的正演计算公式。
使用平面磁场的垂直分量Z(x，y)的极化场谱，
第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次，耗时22分钟。从图中可以明显看出，反演深度等值线(虚线)与理论模型深度(实线)几乎一致。
MIDI方法的一个特点是不需要对面积场值进行滤波，在反演过程中也不需要对深度或场值进行滤波，这可能会改变场源的特性。因此反演计算速度快，深度分辨能力强。使用图2b所示的波状界面，将MIDI方法与Parker方法[4，10]进行比较。用相同的磁化强度计算磁场值，然后从磁场反演界面深度起伏。反演结果分别如图2c和2d所示。两种方法的计算时间都限制在20分钟以内。显然，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— MIDI反演结果与Parker方法相同。深度单位:千米
图3平均界面深度的调整和计算
当平均深度接近真实值时，平均深度波动(Ria)趋于0。
5呼和浩特-张家口地区居里深度的反演
地壳磁层下界的深度，即居里深度，由地壳温度场决定，是地下热状态的重要指标。这对研究深部地质构造、地震学和矿产预测具有重要意义。我们选择了53200km2范围内具有重要深部构造和上部矿产价值的呼和浩特-张家口地区进行居里深度分布的反演研究。图4a显示了该地区太古宙地层、岩浆岩和侵入玄武岩的分布。太古代和玄武岩有很强的磁性。图4b是该区域的航磁异常图。
图4呼和浩特-张家口地区太古宙地层分布图(a)和航磁异常图(b)。
异常值单位:新台币
(1)根据本区岩石和地层磁性统计结果及相关已知地层深度结果，大致得出不同深度的磁化变化数据，遗传算法[5]选取的磁化函数为:
第30届国际地质大会会议录第20卷地球物理学
(2)为了提取居里深度的磁场信息，采用正则化滤波(因子L=40Km)消除浅层磁变化的影响。滤波结果被认为是居里深度界面起伏和磁基顶界起伏的综合异常。
(3)根据已知的地层深度数据和人工地震剖面解释[7，11，13]，从已知的太古代地层顶界面(Ar)位置插值得到全区磁性基底顶界面起伏的大致深度。其平均深度为114以西3.2公里，114以东3.5公里。这种深度分布充当磁性层上界面。
(4)居里深度经计算调整为33.0km。根据反演迭代的中间结果，迭代前后深度值之差为5%作为迭代终止条件，得到的反演结果如图5所示。可以看出，居里深度的主要凸起与该区主要的深大断裂构造相一致。
图5呼和浩特-张家口地区居住深度变化图。
深度等值线单位:千米
6蔡家营地区太古界顶面反转
华北地台北缘的金、铅、锌等多金属矿产一般与太古宙地层(有时与元古代地层)的起伏、断裂构造和岩浆活动(特别是燕山期侵入岩)密切相关。太古宙地层极厚，与上覆地层相比具有很强的磁性。因此，利用航磁异常研究隐伏构造与太古宙地层的顶界面具有良好的物理前提。我们选择了正在进一步深化和扩大找矿的张家口市蔡家营地区，范围26112km2，进行了研究。对1/20万的航磁异常进行滤波，消除浅层地表干扰[12]，用上述居里深度结果参与计算。蔡家营地区太古宙顶部深度波动结果见图6。深度变化在某些剖面上与布格重力异常遗传算法反演结果几乎相同[5]。结合本区地质构造、岩浆岩分布和重力异常，我们对反演结果解释如下。
(1)蔡家营、土城子南部有老地层隆起。蔡家营矿位于老地层周围的隆起凹陷内，现蔡家营矿位于凹陷的次级平缓段。这个凹陷对应于相对低重力的异常。由于本区燕山期岩浆岩密度较低，该凹陷应存在燕山期隐伏岩体。在地表可以看到岩石变形。同样，在土城子以南和张北以西的隆起南部也有类似的平缓地段。这两个地区应该是寻找金、铅和锌等多金属矿床的有希望的地区。
(2)尚义-张北-土城子之间及其西北部的新生代玄武岩层应较薄。这可以从航磁图和滤波结果中得到清晰的解释。土城子西南的太古代隆起应该在这层玄武岩之下。
(3)尚义和康保之间应该有很深的断层(图6)。这从航磁异常图和太古代天顶起伏反演中可以看得很清楚。断层两侧的磁场特征与太古代深度完全不同。
在中国科学院刘光鼎院士的支持和指导下，教授对蔡家营地区的地球物理解释工作给予了很大的帮助，中国地质大学古地磁实验室谭成泽教授等老师，中国科学院地球物理研究所姚长礼讲师和郝副研究员对研究工作给予了很大的帮助。在此表示衷心的感谢。
图6蔡家营地区太古代地形图
虚线为假定的深断裂，深度数据单位为km。
参考
基洛奇霍夫斯卡娅，3。A..地壳的磁性不均匀性。地球物理学报(俄文)，1986，8(5):3~23。
关之宁，安玉林。区域磁异常的定量解释。北京:地质出版社，1995438+0。
郭五林。磁异常区域分量的岩石学解释。物化探翻译系列，1987，51 (4): 26 ~ 32。
[4]帕克..潜在异常的快速计算。地球物理学，J.R.Soc，1973，31:445~447。
张明华。华北地台北缘重磁数据处理新方法及其应用研究[博士学位论文]。中国地质大学研究生院(北京)，1995。
[6]汉森先生..日本全国居里点深度分析。第53届年会，SEG，1983。
地矿部华北地台物化探遥感填图委员会航磁重力解释组。华北地台(1: 1000000)航磁和重力解释成果报告。36860.68666686666
[8]汉森等人..日本及其邻近地区的磁异常。《磁学与地电杂志》特刊，1994，46(6)。
穆时敏，沈宁华，孙。区域地球物理数据处理方法及其应用。吉林:吉林科学技术出版社，1990。
[10]帕克和休斯蒂斯..地形存在时磁异常的反演。地球物理研究学报，1974，79(11):1587 ~ 1592。
张文质王国夫。华北大陆架北缘重磁异常解释。物化探，1992，16 (1): 31 ~ 37。
王继伦。二维最优线性数字滤波器的设计原理。地球物理学报，1977，20 (2): 159 ~ 172。
余，马兴元。华北地区航磁图像处理结果及地震构造解释。地震地质，1989，11 (4): 5 ~ 14。