利用迭代反演方法解决变磁性磁界面反演问题

张向宇，吴健生

（1.广州海洋地质调查局 ，广州 510760；2.同济大学 海洋地质国家重点实验室，上海 200092）

南海是环太平洋西部的最大边缘海，面积约350×104km2。深入开展南海北部中生代沉积盆地研究，对于开拓我国海洋油气勘探新领域具有十分重要的意义。近年来，随着勘探工作的开展，由于地震方法受射线覆盖强度、反射能量有限及分辨尺度等因素的限制，使地震资料难以清晰地反映盆地凹陷下埋深较大的基础层的构造，从而影响对基底结构的了解和构造演化的推断，因此使得重磁勘探的手段得到了进一步的重视。

南海地区经历了多次构造运动，使得地层岩石磁性在横向上呈现较大差异，目前常用的磁界面反演方法如欧拉算法、频谱分析法、维纳反褶积法以及特征点法和切线法等，都是针对常磁性磁界面的计算，若直接将其应用于南海北部磁性界面的反演中必定会带进很大误差，这就给常用反演方法的应用提出了挑战。通过以往的研究，对目前较常用的迭代反演方法进行改进，获得适用于变磁性磁界面反演的方法，为类似南海地区这种磁性变化较大区域的磁性基底反演提供了一种解决方法。其次磁性界面可能与基底界面有关联，这样通过反演磁性界面可以得到基底的情况，这也为进一步研究地质—地球物理解释工作提供了依据。

国内、外有关学者对磁界面迭代反演方法有了较多研究，国外方面，Parker［1］在1972年提出了界面正演方法；Oldenburg D W［2］提出了起伏界面的迭代反演方法，此方法一经提出，便得到了众多地球物理学者的重视，并得到了广泛运用；Guspi［3］在Oldenberg提出方法的基础上，提出了一种非线性迭代反演法；Barbosa［4］对Guspi提出的方法进行了改进，解决了不宜编程实现的缺点。国内方面，赵百民等［5］对用于磁界面反演的一般方法如Parker迭代反演方法、欧拉算法等进行了简要介绍；相鹏等［6］进行了基于Parker算法的磁性双界面正、反演研究，并应用理论模型进行研究验证，且运用到实测数据中；董焕成等［7］、王家林等［8］、管志宁等［9］也在相关书籍中分别对界面迭代反演、直接反演法、统计反演法、线性反演法等反演方法的常用方法原理进行了介绍和应用。基于前人对磁界面反演方法的研究，为本文变磁性磁界面迭代反演提供了一定的依据。

1 磁界面迭代反演方法原理

1.1 常磁性磁界面迭代反演方法

设界面以下磁化强度恒定，以上无磁性（M ＝0），取过平均深度的水平面为起算平面，如图1所示。根据引力位公式及场位关系式，考虑界面起伏远小于平均深度，得到垂直磁化垂直磁异常公式（在xoy平面上）为式（1）。

图1 磁界面反演计算示意图Fig.1 Map of inversion for magnetic interface

其中 h0为界面平均深度；Δh＝h－h0为界面相对平均深度的起伏；Mz为磁化强度M 在z方向上的分量；μ0为真空磁导率。

式（1）写成近似积分形式为式（2）。

式（2）为2 m×2 n个垂直磁化有限延伸直角棱柱体的垂直磁异常之和，即垂直磁化磁界面垂直磁异常的近似公式。当已知磁异常值和磁化强度值时，通过解线性方程组就可得到界面起伏值，实际计算时应采用迭代法计算。其步骤为：①将实际的Za值作为迭代的值，并给出Δhij的初值Δ；②将Zadie（k，l）带入式（2）解出迭代中的 Δhij的改正值δΔhij，将δΔhij与相加后再代入式（2）计算的理论值；③计算理论值与实际值的均方误差，当满足迭代终止条件时即完成了反演计算［7］。

1.2 变磁性磁界面迭代反演方法及改进

式（2）表示可将磁界面近似看作由2m×2n个磁化强度相等的垂直磁化直立长方体组合而成。对于变磁性磁界面，每个长方体的磁化强度值不等，则垂直磁化垂直磁异常公式变为式（3）。

其中 Mij为第（i，j）个直立长方体的磁化强度。

同样采用上述提到的迭代法进行反演即可得到界面深度。式（3）可理解为某一观测点（x，y，z）处的磁异常值是组成磁界面的所有长方体在点（x，y，z）处产生的磁异常的和。将式（3）改写成矩阵形式为式（4）。

其中 H为界面深度向量；M 为磁化强度向量；Z为垂直磁化垂直磁异常向量。

在实际计算时，已知Z和M 矩阵，通过矩阵求逆解出H矩阵的结果，便实现了对界面深度的反演。从这个反演过程不难发现，若要达到可以将每一个长方体的磁化强度值和高度值看做不变的近似条件，则划分的长方体数N必定要足够大才可行，而随着N的增大，矩阵计算的工作量也大大增加，这样就不适宜实际计算的实现，考虑到这个因素，我们对变磁性磁界面反演进行改进，获得了一种适宜实际计算的方法。

针对这种情况，需要对变磁性磁界面迭代反演方法进行改进。考虑到以往研究相类似问题时，运用滑动时窗的方法较为便利，故我们也由此入手。若界面的磁化强度值和界面起伏都变化不大的时候，我们可以做进一步近似。取适宜大小的窗口，将窗口处的界面深度值和磁化强度值看做是恒定的，为窗口内所有长方体的深度和磁化强度的平均值，在每一个窗口处通过解矩阵求得界面深度值，当窗口滑动遍历整个计算区域，就得到了整个界面的反演结果。这样通过引入滑动窗口的方式，将大大减小反演计算量。

作者采用的变磁性磁界面迭代反演方法，就是建立在这个思想上对常磁性磁界面迭代反演方法改进得到的。首先输入研究区化极后的Za磁异常数据和事先通过研究区地质资料得到的平均界面深度和磁化强度值；为了减小边界影响，需要将异常下延到平均界面深度处，计算每一个窗口处对应平均界面埋深的磁异常线性正演矩阵的解，得到每个滑动窗口处对应的界面深度，这样通过滑动窗口遍历整个研究区域就可以得到整个区域的界面深度；通过对得到的界面深度值进行磁异常正演与输入的实际磁异常进行比较，判断均方误差是否满足迭代结束条件，若不满足迭代结束条件，则求取界面深度修正值，返回继续进行迭代计算，直到满足迭代精度或迭代次数要求后即得到界面深度反演结果。整个流程见图2。

在反演过程中，为了减小边界影响，在计算窗口处矩阵解前，需要先将异常下延到平均界面深度处，这里采用的是加入正则化因子进行空间域向下延拓的方法。

异常场下延可表示为矩阵形式［8－9］：

其中 Z0为已知场值矩阵；A为延拓系数矩阵；ZH为下延到H深度处平面上场值。

则式（5）可以写成式（6）的形式。

图2 变磁性磁界面深度反演流程图Fig.2 Flowchart of inversion for varying magnetic interface depth

引入一个正则化因子α，有

其中 I为单位矩阵。

由式（7）选取合理α值便可计算得到下延异常，若α的取值过小，则无法较好地控制下延发散的问题；若α的取值过大，则将导致反演误差增大。这样，通过参数α的引入，较好地完成了磁异常的下延。

2 模型正演验证

为了验证方法的可靠性和适用条件，选取模型进行正演计算。作者将选择规则球体模型和模拟实际地形的仿真模型分别进行正演计算。

2.1 单一规则球冠模型正演计算

取球体顶点高度为－4km，平均深度为－4.68 km的单一球冠，测线间距和测点间距均为5km，测线数和测点数均为64，界面磁化强度取如图3（a）所示渐变式，其平均值为10A／m，这里分别取界面磁化强度为真值（变磁性）和为平均值（常磁性）时计算反演结果，得到磁异常正演结果、变磁性反演结果及由常磁性反演结果见图3。

为了更清晰地对比反演结果，取过球冠顶点处的测线数据进行误差分析，剖面反演结果对比图见图4。统计测线数据的反演误差，得到变磁性反演结果相对误差为0.2%～2.4%；而常磁性反演结果相对误差为0.6%～4.5%。根据图3（d）等值线图显示，常磁性反演结果使球缺走势随界面磁化强度渐变方向倾斜，而变磁性反演则没有出现这个问题，同时变磁性磁界面反演的结果比常磁性磁界面反演结果相对误差更小，由此可说明变磁性反演的优势。

图3 单一球冠模型反演结果Fig.3 Inversion of single spherical cap model

图4 球冠模型剖面数据对比图Fig.4 Profile comparison of spherical cap model profile

2.2 仿真模型正演计算

仿真模型模拟的是海盆和海山共存时的海底地形，测线和测点间距均为2km，地形平均高度为－4.25km，其地形等值线图如图5（a）所示，区域磁化强度呈线性渐变，平均值为50A／m，其正演磁异常见图5（b），变磁性反演结果见图5（c）。

同样为了方便分析方法的效果，取如图5（c）黑色线所示中心测线上的数据分析反演误差，因测线经过测区地形起伏变化最大的地方，故该测线较具代表性，剖面数据对比图见图6。分析数据得到相对误差为0.6%～7.4%，其中误差最大点出现在地形起伏顶点处，可以说该反演方法对地形有一定的“抹平”作用，同时整体测线相对误差大小较合理，说明针对接近实际情形的仿真模型，该方法也取得了较好的效果。

2.3 影响因素分析

在对几组模型进行反演验证的过程中发现，观测面高度的选取、测线测点间距的选取和计算范围的大小也会对反演结果产生一定的影响，下面通过单一球冠模型对这几方面影响因素进行简要说明。

2.3.1 观测面高度影响。

对上述单一球冠模型固定其形状，取磁化强度为常值10A／m.，取观测面高度分别为3km、5km和8km，计算反演结果。

将得到的三个反演结果分别与模型地形数据进行对比：①当观测面高度为3km时，反演误差最大为7.1%；②当观测面高度为5km时，反演误差最大为2.2%；③当观测面高度为8km时，反演误差最大为1.2%。由此可以看到，随着观测面高度的增大，反演误差逐渐减小。所以观测面高度取值不宜过小也不宜过大，要根据研究区的地质情况进行合理选择。

2.3.2 测线间距影响

当测线间距取值不合理时，也会影响反演结果。为了考察测线间距取值的不同对反演结果的影响，仍采取单一球缺模型，球缺顶点高度仍为－4km，底界面深度为－5km，磁化强度为10A／m，取测线间距和测点间距相同，分别取测线间距为1km、3 km、5km、10km进行变磁性磁界面反演。

分析反演结果的相对误差：①当测线间距为1km时，反演结果误差最大为10.2%；②当测线间距为3km时，反演误差最大为5.8%；③当测线间距为5km时，反演误差最大为2.2%；④当测线间距为10km时，反演误差为1.1%。

图6 模拟地形剖面数据对比图Fig.6 Profile comparison of simulative topography

由此得出，测线间距越大，反演误差越小，但考虑到计算的数据量，测线间距也不宜取的过大，要根据计算区域的大小，选择合理的测线间距进行计算。

2.3.3 边界影响

在计算时应将根据计算数据的情况调整计算区域大小，不同的计算区域，其边界效应对反演结果的影响是不同的。下面仍然固定球缺模型深度、大小、磁化强度值及测线间距值，变化测点测线数，分别取计算区域为220km×220km、270km×270km、320km×320km和420km×420km计算反演结果。

统计反演结果相对误差，当计算范围为220km×220时，反演误差最大为5.4%；当计算范围为270km×270km时，反演误差最大为3.1%；当计算范围为320km×320km时，反演误差最大为2.2%；计算范围增大到420km×420km时，反演误差最大为1.2%。由此看到，随着计算范围的增大，边界对结果的影响减小，则反演误差也逐渐减小。

根据上面的模型分析及影响因素分析，可以对本文采用的变磁性磁界面迭代反演方法进行以下总结：

（1）通过规则形体模型和模拟实际地形地质模型的验证，得到较理想的反演结果，说明了该方法的可靠性。

（2）采用该方法反演时，对地形存在一定的“抹平”作用，因此对地形起伏较陡峭的区域不适宜应用该方法。

（3）运用变磁性磁界面迭代反演方法的前提是，界面起伏远小于界面平均深度，并且界面磁化强度不能有正负相差较大的情况，否则会导致反演结果无法收敛，从而无法达到迭代反演的目的。

（4）在运用该方法时，存在观测面高度、测线测点间距、计算范围等几个主要因素影响反演结果。因此在对实测数据进行反演计算时，应充分注意这几个因素。

（5）观测面高度取值越大，反演误差越小，因此应该根据研究区实际情况，在合理范围内选取最大可取观测面高度值进行反演。

（6）随着测线、测点间距取值增大，反演误差减小，但为了保证计算的数据量，测线、测点间距不宜取值过大，应根据计算范围的大小选择合理的点距、线距进行反演。

（7）计算范围越大，边界影响越小，则反演误差越小，但为了避免计算量得增大，影响反演计算的速度，也不宜将计算区域扩的过大，所以在对实测数据计算前，应该根据边界影响大小适当调整计算范围。

3 实测数据磁性基底反演

通过对模型进行反演结果对比后，对变磁性磁界面迭代反演方法有了一定的了解，将该方法运用于南海东北部某研究区实测数据中，计算该区域磁性基底深度。研究区地形没有起伏较尖锐的点，并且各测点间磁化强度无变化较剧烈的点，适合应用变磁性磁界面迭代反演方法进行磁性基底深度反演。

在计算前，需要预先得到研究区由磁性基底引起的磁异常及区域磁化强度值，这里先通过化极及小波分解获得由磁性基底引起的Za磁异常（图7），再由研究区地质地球物理资料及相关方法计算得到区域磁化强度值（图8）。将上述数据代入变磁性磁界面迭代反演程序中进行计算，得到研究区磁性基底界面深度分布图，反演结果见图9。图7～图9中棕色线段标示区域断裂，紫色线段标示二级构造单元，各构造单元名称见图例。

由图7－图9进行对比总结：

图7 研究区化极磁异常小波三阶细节场Fig.7 Wavelet third－order details for Magnetic anomaly of research area

图8 研究区基底磁异常Fig.8 Magnetic anomaly for basement of research area

图9 研究区磁性基底深度反演结果Fig.9 Inversion for magnetic basement depth of research area

（1）在图7中，研究区Za异常呈条带状，近东北－西南方向延展，并出现正磁异常、负磁异常条带间隔出现的现象，在东沙隆起板缘区域磁异常呈现最大值约为100nT，在白云坳陷区域则出现磁异常负值，为－20nT～－100nT。

（2）在图8中，研究区磁化强度形态也呈条带状展布，与磁异常延伸方向相近，在图7中对应正磁异常条带的位置（东沙隆起板缘），区域磁化强度近乎呈常值，约为2A／m；从磁化强度的这一高值带的分布特征和强度来看，推论可能是侵入基底的火成岩的反映，其分布位置和展布方向与阳江—统暗沙断裂东南段非常一致，认为其发育与这一断裂有密切关系；而对应磁异常负值的区域（白云坳陷）磁化强度则出现负值，最小达到－10A／m，并向研究区东南方向逐渐减小。

（3）图9中磁性基底在东沙隆起板缘区域深度近乎不变，约为6.7km，延展区域也与磁异常形态相近，向近东北－西南方向延伸，同时对应区域磁异常和磁化强度负值的区域（白云坳陷）磁性基底深度值则较大，最大约为10km左右；两部分中间区域磁性基底深度则较小，为4.5km左右。

从整体反演的磁性基底深度结果看到，南海东北部地区除局部地区外，总体呈现大面积宽缓的等值线特征，一直延伸到洋陆分界处。从磁性基底分布的宏观特征来看，东北部陆缘的构造性质更倾向于非火山型［10－15］。

4 结论

（1）作者使用滑动时窗，对常规磁界面迭代反演方法进行改进得到适用于变磁性磁界面的迭代反演方法，这是与以往磁界面反演方法不同之处，而应用模型进行验证取得了较好的效果，说明了该方法的可靠性较好。

（2）在运用本文提到的变磁性磁界面迭代反演方法前，需要预先得到研究区化极磁异常Za数据、磁化强度数据和界面平均埋深值。

（3）在对实测数据应用本文提及的变磁性磁界面反演方法反演时，要注意观测面高度、边界影响等因素的影响。观测面高度取值要合理，计算范围要根据实测数据正演磁异常情况进行一定调整，防止边界效应影响反演结果。

（4）作者研究变磁性磁界面迭代反演方法具有地形“抹平”作用，对于起伏较大的地区，在地形最值处反演结果有一定误差，在应用实测数据进行计算时要注意这个特点，对于这一缺陷的改进在今后的研究工作中有待解决。

（5）通过实测数据计算，得到研究区磁异常及磁化强度均呈近东北—西南方向条带状特征展布，通过变磁性磁界面反演得到的磁性基底分布特征与磁异常特征相似，平均深度为6km左右，最大深度约10km。这为今后地质－地球物理进一步的研究提供了一定的依据。

［1］PARKER R L.The rapid calculation of potential anomalies［J］.G.J.R，1972，37（3）：447－455.

［2］OLDENBURG D M，MOGILIVILAY P R.Generalized subspace methods for 1arge－scale inverse problems［J］.Geophys J Int，1993，114：12－20.

［3］FERNANDO GRUSPI.Nonilerative nonlinear gravity inversion［J］.Geophysics，1993，58（7）：935－940.

［4］BARBOSA V C F，SILVA B C，MEDEIROS W E.Gravity inversion of basement relidf using approximate equality comstrains on depths［J］.Geophysics，1997，62（6）：1745－1757.

［5］赵百民，郝天珧.反演磁性地质界面的意义和方法［J］.地球物理学进展，2006，21（2）：353－359.

［6］相鹏，刘展.双界面模式Parker算法磁性界面正反演方法［J］.中国石油大学学报：自然科学版，2009，33（1）：37－43.

［7］董焕成.重磁勘探教程［M］.北京：地质出版社，1993.

［8］王家林，王一新，万明浩.石油重磁解释［M］.北京：石油工业出版社，1991.

［9］管志宁.地磁场与磁力勘探［M］.北京：地质出版社，2005.

［10］姚伯初，万玲.中国南海海域岩石圈三维结构及演化［M］.北京：地质出版社，2006.

［11］刘昭蜀，赵焕庭，范时清，等.南海地质［M］.北京：科学出版社，2002.

［12］郝天珧，徐亚，赵百民，等.南海磁性基底分布特征的地球物理研究［J］.地球物理学报，2009，52：2763－2774.

［13］陈洁，高德章，温宁，等.南海磁场特征研究［J］.地球物理学进展，2010，25（2）：376－388.

［14］江凡.南海东北部中生界分布综合地球物理研究［D］.上海：同济大学，2010.

［15］郝天珧，黄松，徐亚，等.南海东北部及临区深部结构的综合地球物理研究［J］.地球物理学报，2008，51（6）：1785－1796.