0 引言

磁源体深度反演是地球物理位场数据解释中十分重要的工作之一(Barnett, 1976; 谢汝宽等, 2016).准确地反演场源埋深有助于了解磁源体结构以及地下地质信息，提高解释质量(Portniaguine and Zhdanov, 2002; Thurston et al., 2002; Lee et al., 2010).此外，磁源深度的反演在石油勘查中确定沉积基底的深度、矿产勘查中确定磁性矿体或岩体的埋深等应用中发挥重要作用(史辉等，2005).

场源几何参数反演和物性反演是常见的两类磁场反演方法(杨文采，1997).几何参数反演是在给定物性参数大小的基础上，利用地面观测异常来拟合几何形态大小，通过几何体的形状大小来模拟目标体的分布规律(Talwani and Ewing, 1960; Barnett, 1976).物性反演则是将观测区域地下半空间离散化成规则的网格单元，通过反演方法确定各离散单元的物性值，由物性的分布确定场源的实际分布情况(Li et al., 1996).但是此类方法常常存在多解性，需要添加合适的约束条件才能获取较好的结果.几何反演和物性反演方法需要同时考虑场源的深度和形态分布特征，利用场源的分布获取深度结果.此外，还有一些无需考虑场源具体形态，直接求取场源深度的较为简单的方法.例如，特征点法、切线法、沃纳(Werner)反褶积方法(Hartman et al., 1971)、功率谱法(Spector and Grant, 1970)、欧拉(Euler)反褶积法(Thompson, 1982)、局部波数法及改进方法(Salem and Smith, 2005; Ghasemiannia and Oskooi, 2017)等.其中，欧拉反褶积是一种快速反演方法，无需对地质模型作任何假设，也无需进行化极处理，具有较广的适用范围，解决了众多实际问题，为我国地质问题研究提供了帮助，已经成功地应用在塔里木盆地区域磁异常反演以及磁源体分布研究中(杨文采等, 2012).但是，该方法并没有充分挖掘位场本身的几何信息.此外，构造指数的选取对反演结果有重要的影响.如何选择合适的构造指数，对于欧拉反褶积反演方法来说是一个至关重要的问题.

事实上，谱矩可以对位场数据的表面形貌做较为细致的描述，例如粗糙度、曲率等(黄逸云, 1984, 1985).因此，我们希望借助于谱矩的工具对磁异常数据进行分析，以达到深度反演的目的.孙艳云、杨文采等利用二阶谱矩及统计不变量成功地提取了地壳变形带信息，提出了区域重力场多尺度刻痕分析方法，并将这些方法运用于中国大陆重力场刻痕识别中(孙艳云和杨文采, 2014; 杨文采等, 2015a, b, c; Sun et al., 2016).付丽华等提出了一种基于谱矩的描述表面数据粗糙度的特征因子，将其用于自由空气异常数据中来圈定地貌的类型(付丽华等, 2017).因此，若能找到位场本身的几何信息与场源埋深的关系，则可以对磁源体深度进行更加精准的反演.

本文将提出基于算术平均顶点曲率(四阶谱矩统计不变量)的磁源体埋深估计方法.介绍了谱矩用于表面形貌识别的原理，探讨了磁场的曲率信息与场源埋深的关系.然后通过理论模型展示该方法的反演效果，最后将其用于塔里木盆地航磁异常的反演中.结果表明，这种新的磁源体深度反演方法是有效的.与欧拉反褶积方法相比，无需调节参数，而且快速简单，同时具有较好的反演精度.

1 谱矩方法计算磁源深度的原理 1.1 谱矩

区域位场数据T(x, y)的p+q阶谱矩定义为其功率谱密度函数的p+q阶矩，定义为(Thomas, 1987; Yanagi et al., 2001; 杨叔子和吴雅, 2007)

(1)

其中，G(f_x, f_y)是位场数据T(x, y)的功率谱密度函数.G(f_x, f_y)可以定义为自协方差函数R(s₁, s₂)=ET(x, y)T(x-s₁, y-s₂)的傅里叶变换(Stoica and Moses, 2007)，即

(2)

对(2)式求Fourier逆变换，可以得到

(3)

并且对(3)式两边求r=p+q阶导数，

(4)

当r为偶数时，

(5)

因此，

(6)

当r=4时，m₄₀, m₀₄为其中两个四阶谱矩，文献(杨叔子和吴雅, 2007)给出具体的求解公式：

(7)

即数据的四阶谱矩表示随机表面的曲率方差，分别表示T(x, y)二阶导数的自协方差，描述的是位场数据表面的曲率情况，m₄₀为x方向上的曲率的方差；m₀₄为y方向上的曲率的方差.

由于上述表面的谱矩依赖于坐标系，即它随着坐标系的旋转而改变.因此，需要定义与坐标系旋转无关的统计量来刻画表面的几何特征，尤其是对场源深度比较敏感的曲率特征.

1.2 表面统计不变量与算术平均顶点曲率

任何点的顶点曲率规定为该点上的主曲率的平均值，即定义为

顶点曲率

S_sc表示位场表面评价在区域0≤x≤l_x，0≤y≤l_y内所有点的顶点曲率平均值(arithmetic mean summit curvature of the surface)，即

(8)

Thomas(1982)推出上式可以由四阶谱矩得到：

(9)

事实上，S_sc为描述区域内算术平均顶点曲率，即位场数据区域内的弯曲程度，与场源深度特征关系较为敏感.接下来，将以球状磁源体和板状磁源体为对象，研究磁源体的埋深与曲率之间的关系.该方法只针对垂直磁化模型，在实际应用中，不考虑剩磁或剩磁较弱的情况下，可以采用当地地磁场参数先化极再进行反演.

1.3 球状磁源体埋深理论结果

假设单个球状磁源体中心埋深为R，磁化强度为M，体积为v，磁矩m=Mv，球心坐标为(0, 0, R)，则磁异常公式为

其中，μ₀为真空的磁导率，I为磁化强度倾角，A′为剖面与磁化强度水平投影夹角.

当磁化强度倾角I=90°时，

计算磁异常的曲率：

因此

又因为  

所以，磁源体埋深可以由下面公式计算：

(10)

1.4 板状磁源体埋深理论结果

垂直板状磁源体的磁异常公式为

(11)

其中，μ₀为真空的磁导率，M_s为磁化强度.b为一半板长，h为垂直板上顶埋深.

事实上，公式(11)可以变形为

(12)

其中，

并且，

所以，则当b≪h时，板状磁源体的埋深估计结果为

(13)

2 理论模型数据实验

为了检验上述方法的正确性，分别对球状磁源体和垂直板状磁源体的埋深进行理论试验，并将反演结果与欧拉反褶积方法做对比分析.首先，对磁异常数据ΔT运用公式(9)求取对应的曲率值，然后，利用公式(10)和公式(13)分别对球状磁源体和板状磁源体进行埋深的估计.

2.1 球状磁源体埋深估计实验

采用图 1的球状磁源体模型进行理论试验.球体半径为8 m，磁化强度倾角为90°，x轴和y轴方向点距均为1 km.图 1a，图 2a，图 3a分别是单个球状磁源体在不同埋深情形下产生的磁异常数据，球心坐标分别为(0 m, 0 m, -20 m)、(0 m, 0 m, -30 m)、(0 m, 0 m, -40 m).图 1b，图 2b，图 3b是不同埋深下的球状磁源磁异常对应的曲率计算结果.可以看出，不同的埋深情形下，基于谱矩的曲率结果都是在磁源体埋深垂直上方的位置，即横纵坐标为(0 m, 0 m)处达到最大值.而且，随着磁源体的深度值越大时，磁异常最大值出现了明显的下降：150→40→20，相应的曲率结果也出现了较为明显的下降：2→0.3→0.05.根据公式(10)可以利用磁异常值和曲率值进行磁源体埋深值的估计.当磁源体埋深为20 m时，公式(10)反演出的结果为20.1872 m，当磁源体埋深为30 m时，基于谱矩的反演结果为30.1249 m，当磁源体埋深为40 m时，其反演结果为40.0938 m.

表 1给出了当球状磁源体埋深从20 m到80 m时基于谱矩的方法与基于欧拉反褶积方法估计出的埋深结果及相对误差.其中，构造指数是根据王明等(2012)的建议选为3，总梯度模倍数选为0.6.当磁源体埋深为20 m，30 m时，基于欧拉反褶积反演出的估计值比本文所提方法的相对误差要小，但是当埋深增加时，本文所提的方法估计精度越来越高，而欧拉反褶积方法反演出的结果相对误差确越来越大.当埋深超过60 m时，原有的总梯度模倍数参数已经不能很好的反演，此时，我们需要调整该参数，将该参数调节为0.01.在新的参数下反演出来的结果相对误差仍然超过了2%.从该实验中可以看出，与欧拉反褶积方法相比，本文提出的基于谱矩的磁源体深度反演方法无需调节参数，在理论上能够保持高精度的反演结果(与真实值相比相对误差不超过1%).

2.2 板状磁源体埋深估计实验 2.2.1 单个板状磁源体埋深估计实验

考虑单个垂直板状磁源体的埋深估计.实验中采用的模型是直立垂直向下延伸的板状体，半板长b=4 m，理论埋深h分别取表 2中的理论值，用公式(13)反演出来的埋深估计值见表 2.板状磁源体理论深度值为30 m时，埋深估计为30.51 m，相对误差为1.70%.随着埋深的增加，基于谱矩的方法反演出的结果精度也越来越高，当理论埋深为100 m时，反演估计结果为100.15 m，相对误差为0.15%.理论埋深为200 m时，反演估计结果为200.08 m，相对误差为0.04%.从表 2可以看出，基于谱矩的反演结果与真实的理论值非常接近，这与理论推导是十分吻合的，当理论埋深值h越大于半板长b时，公式(13)的估计精度越高.然而，基于欧拉反褶积的方法受构造指数的影响特别大，不同的构造指数会给结果带来较大的偏差.表 2中给出了不同埋深结果下欧拉反褶积的结果与对应的构造指数值.可以看出，欧拉反褶积的方法应用效果与构造指数有着密切的关系.而本文所提的方法则无此问题.

2.2.2 多个板状磁源体埋深估计实验

本实验利用谱矩的方法确定多个板状磁源体埋藏深度.检验模型由4个直立下延有限的板状体组成，每个板状体的水平位置、上顶埋深、板长均不相等，地磁倾角为90°，测线方位角为0°.四个板状磁源体水平中心位置分别为1000 m，2000 m，3000 m和4000 m.上顶埋深理论值分别为80 m，200 m，300 m和100 m.半板长分别为15 m，30 m，50 m和25 m.磁异常测线图以及反演结果如图 4.图 4上部分子图为磁异常测线值，下部分子图为实际板状体位置和基于谱矩的方法反演的结果，四个板状体埋深反演值分别为84.55 m，215.30 m，316.91 m，108.57 m，与理论模型基本一致.

3 在塔里木盆地地区的应用效果

塔里木盆地地区的航磁数据来自国家地质总局航空物探大队和各大石油公司，不同测区的数据经过标定和同化处理.该数据采用兰勃特坐标系，经纬度范围分别为：x=-451000~900000，y=3250~711550，数据网格间距为：2 km×2 km，点数为678×356.首先将航磁异常利用中纬度常规化极方法进行处理，获得的化极异常图如图 5所示.

为了验证基于谱矩的方法对于磁源体深度反演中的有效性，首先将航磁异常表面数据进行曲率求取，求得的结果见图 6.从图 6中看出，曲率值越大，则显示出该位置处磁源体的埋深值较浅，相反的，曲率值越小则此处磁源体影响较小，即埋深较深.接着，我们将塔里木盆地地区航磁数据按照曲率值的大小分为三类:0~8(曲率值小类)，8.001~38(曲率值中类)，38.001~70(曲率值大类).然后分别采用球状磁源体埋深反演公式或者板状磁源体埋深公式，计算磁源体的埋深值.具体的方法如下：(1)对于曲率值为0~8之间的数据，这时磁源体在沉积层中，故而采用球状磁源体埋深公式进行计算.(2)对于曲率值为8.001~38之间的数据，这时大多数磁源体还在沉积层中，采用球状磁源体埋深公式进行计算.(3)对于曲率值为38.001~70之间的数据，大多数磁源体来自结晶基底，采用板状磁源体埋深公式进行计算.

在谱矩方法计算出的塔里木盆地磁源体深度埋藏图(图 7)中，深度共分为3个等级：0~5 km，5~10 km，10~20 km，分别用红、绿和蓝三种颜色标识.可以看出，比较欧拉反褶积结果(杨文采等，2012)，基于谱矩的磁源体深度反演方法得到的结果并非是一个一个的点状不连续的标记，而是连续的区片.究其原因，区域磁场由大量异常组成，相互连接，计算的关键之一是避免相连异常干扰，而谱矩方法用曲率计算磁源体深度，能最接近磁源体中心, 反演出的深度结果受相连异常干扰较少.尤其是柯坪地区二叠纪火成岩体，在磁性图上的边界清楚，为准确的解释提供了依据.总的来说，塔里木盆地深度为0~5 km深度的磁源体分布面积相对较小，主要出现在柯坪隆起，巴楚隆起西北部，深入至古董山断裂，阿尔金山前带和孔雀河斜坡西段(图 8).塔里木盆地深度为5~10 km主要集中在麦盖提—巴楚—阿瓦提—顺托果勒地区，阿尔金山前带以及孔雀河断裂带上(图 9).塔里木盆地深度为10~20 km表示的磁性体源自本地区的结晶基底(图 10)，分布在全盆地，面积很大.

根据塔里木地区的物性统计结果(杨文采, 1997, 2009; 杨文采等, 2012, 2015c, d, 2015e)，柯坪隆起，巴楚隆起西北部，和田河气田至古董山断裂地区的磁异常主要来自中生代构造隆起地区的二叠纪玄武岩和辉长岩.阿尔金山前带地区的磁性体可能与超基性岩体有关，也与中生代盆地东南缘的造山作用有关联.上述关于塔里木地区的磁性体三维分布的反演结果和地震调查、重力场三维反演结果和大地电磁法电阻率分布的反演结果都是兼容的(侯遵泽等, 2011；于常青等, 2012；瞿辰等, 2013；杨文采等, 2015d, e)，而且增加了磁异常源体分布的信息.比较欧拉反褶积结果可知(杨文采等，2012)，基于谱矩的磁源体深度反演的理论假设更少，计算方法更加简明快捷，取得的结果更加准确精细.因此，我们建议更多使用基于谱矩的磁源体深度反演，完成区域磁异常数据分析和解释.

4 结论

本文基于以随机过程为基础的谱矩理论，研究并提出了磁性体埋藏深度反演的一种方法.该方法运用曲率(四阶谱矩的统计不变量)提取磁异常的几何信息，并将所提取的信息用于深度的反演.具体内容包括：(1)理论上推导出基于谱矩的球状磁源体和板状磁源体的埋深反演公式.(2)理论模型数据试验较好地验证了本文提出的磁源体埋深估计方法的有效性.(3)塔里木盆地地区的磁性体深度反演试验表明，新的方法能很好地反演出磁源体深度，为克拉通沉积盆地磁异常源的深度划分提供了丰富的信息，有助于将它们与形成的地质作用和时代联系起来，为准确解释区域磁异常图提供可靠的依据.上述关于塔里木地区的磁性体三维分布的反演结果和地震调查、重力场三维反演结果和大地电磁法电阻率分布的反演结果都是兼容的，而且增加了磁异常源体分布的信息.与欧拉反褶积方法相比，基于谱矩的磁源体深度反演的理论假设更少，计算方法更加简明快捷，取得的结果更加准确精细，建议推广.