基于双磁性界面和变磁化率模型的居里面深度反演计算

2019-10-23刘卓曾昭发李静何荣钦霍祉君林景宜

世界地质 2019年3期

刘卓, 曾昭发, 李静，何荣钦，霍祉君，林景宜

1.吉林大学地球探测科学与技术学院，长春 130026;2.国土资源部应用地球物理重点实验室，长春 130026

0 引言

居里面是地壳中岩石铁磁性矿物在温度升高到居里点变为顺磁性物质时的深度界面，为磁法勘探的磁性层底界面，表征地下温度场的分布，对地热田评价开发、油气资源的预测、地震火山灾害防治及原生热液矿产勘查等具有重要指导意义[1]。居里面深度计算方法包括谱分析方法和界面反演方法等。谱分析方法[2--3]是通过计算磁异常数据对数功率谱来确定磁源深度，但对数功率谱的复杂性和混乱性难以有效地划分磁源平均深度，具有较大的误差。Parker在位场界面正演计算中引入快速傅里叶变换(FFT)[4]。Oldenburg根据Parker公式，提出频率域深度界面迭代反演方法[5]，改变了过去采用棱柱或长方体模型反演界面深度不连续性，同时计算速度快[6]。Parker--Oldenburg反演算法存在3个问题：①该方法属于单磁性界面反演方法，磁性体通常是由双界面或多界面组成的,利用单界面计算与实际情况不符。王万银[7]和相鹏[8]采用双磁性界面模型来提高反演结果精度。②在反演过程中假设磁化率参数为常数，而实际地质介质磁化率参数在横向和纵向上均具有较复杂的变化。冯娟等[9]采用了垂向变物性参数模型来计算基底界面，提高了反演精度。③在Parker--Oldenburg公式中包含向下延拓因子中的不收敛项，影响计算的稳定性。尽管可以通过添加滤波器来确保迭代的收敛，但将会使信号的高频信息损失，降低反演结果的分辨率和精度。肖鹏飞[10]、张冲等[11]采用多项式迭代方法来改进传统Parker--Oldenburg算法，提高计算稳定性。

针对上述问题，前人研究中分别针对不同问题提出了对应的解决方案，没有考虑综合影响。笔者综合了上述3个问题提出了一种改进Parker--Oldenburg界面的反演方法。通过引入垂向变磁化率因子和双界面模型来提高模型的准确性，同时采用多项式迭代算法来提高算法稳定性，提高反演精度。通过理论模型测试和实测数据应用验证了该算法的可靠性以及结果的准确性，计算的居里面深度更加具有理论和实际应用价值。

1 基于双磁性界面和垂向磁性变化的Parker--Oldenburg反演方法

1.1 基于双磁性界面和垂向磁性变化的磁异常计算

垂向磁化率和双磁性界面模型条件下频率域正演磁异常公式为：(关于公式(1)推导过程详见附录)

(1)

对公式(1)进行重组，得到如下公式：

(2)

频率域位场的向上延拓公式为：

F[U0]=e-ωh0F[U]

(3)

可以看出公式(2)和公式(3)在形式上具有相似性。

1.2 磁异常--深度变换的迭代公式

为迭代计算得出界面反演公式，首先假定地下某一深度处的磁异常初始值设为Δz(x,y,z)(1)=Δz(x,y,0)。利用公式(3)得出地表的异常值计算初值如下：

Δz(x,y,0)(1)=F-1[e-ωh0F[Δz(x,y,z)(1)]]

(4)

结合位场迭代计算方法[12]，可以利用Δz(x,y,z)(1)、Δz(x,y,0)、Δz(x,y,0)(1)获得Δz(x,y,z)(2)，即：

Δz(x,y,z)(2)=Δz(x,y,z)(1)+s(Δz(x,y,0)-Δz(x,y,0)(1))

(5)

以此类推，经过n次迭代计算后的结果：

Δz(x,y,z)(n+1)=Δz(x,y,z)(n)+s(Δz(x,y,0)-Δz(x,y,0)(n))

(6)

当|Δz(x,y,0)-Δz(x,y,0)(n)|≤ε时，认为Δz(x,y,z)(n+1)≈Δz(x,y,z)(n)，对公式(6)推导可得：

(7)

提取深度结果得到：

(8)

公式(8)即为双界面模型约束下的下界面反演公式。

2 数值模型模拟计算

图1和图2为设计的模型。从图1可见模型由位于地下30 km和10 km处的两个磁性界面组成，每个磁性界面均由两个大小不同、方向相反的半球组成，30 km深度处底界面的较大凹的半球直径为40 km，方向向下，最大埋深为40 km；较小隆起半球的直径为30 km，方向向上，最小埋深22.5 km；10 km深度处顶界面的较大凹半球直径也为40 km，最大埋深为20 km；较小隆起半球的直径为30 km，最小埋深2.5 km。目标体磁化率设为M0=1 A/m。为了对比分析，本文另设模型体的变磁化率设为M=M0e-0.178z。分别采用本文提出算法和传统Parker--Oldenburg方法对图1所示模型开展界面深度反演。提取模型俯视图的界面对角线AB(图2)剖面反演结果进行分析。

图1 双磁性界面模型体分布图Fig.1 Dual magnetic interface model

图2 模型体俯视图对角线AB位置图Fig.2 Diagonal AB top view of interface model

图3为传统Parker--Oldenburg方法界面深度反演结果，图4为本文方法的磁异常界面深度反演结果。图中黑实线为磁异常模型值，点划线和黑虚线分别表示变磁化率和常磁化率情况下计算深度值。表1和表2为传统算法与本文提出算法的计算参数对比。

2.1 变磁化率对界面深度反演的影响

对比图3和图4中点划线和虚线的反演结果，可以看出应用同一算法变磁化率模型下的反演结果要优于常磁化率条件下的反演结果。

对比表1的数据可以看出，在应用传统Parker--Oldenburg方法条件下，变磁化率模型反演结果迭代20次剩余RMS误差为0.019 2；常磁化率模型反演结果迭代30次剩余RMS误差为0.019 3，明显提高收敛速度。

图3 传统Parker--Oldenburg算法反演模型下界面深度结果Fig.3 Bottom interface depth inverted results obtained by conventional Parker--Oldenburg algorithm

图4 新算法反演模型下界面深度结果Fig.4 Bottom interface depth inverted results obtained by new algorithm

数值模拟最大埋深/km最大深度位置/km最大埋深/km最小深度位置/km(87,87)点差值/km(36,36)点差值/km收敛准则/km剩余RMS/km迭代次数P--O算法反演40.00(44,44)22.50(90,90)-----常磁化率模型39.47(44,44)25.82(93,93)3.481.300.020.01930030变磁化率模型40.09(44,44)24.83(91,91)2.481.500.020.01920020新算法反演40.00(44,44)22.50(90,90)-----常磁化率模型40.08(43,43)25.46(90,90)2.901.460.020.0052003变磁化率模型40.05(44,44)22.48(90,90)0.182.120.020.0000655

在应用本文提出的算法变磁化率模型反演结果迭代5次剩余RMS误差为6.5×10-5，常磁化率模型反演结果迭代3次剩余RMS误差为0.005 2；变磁化率模型的RMS值要远小于常磁化率模型的RMS值，表明变磁化率因子可以显著地提高新法的反演精度和收敛的速度。

上述对比分析验证了将变磁化率因素引入到界面反演公式中的确可以提高反演结果的精度和可靠性。

2.2 迭代算法对界面反演结果的影响

对比表2的数据可以看出，在应用变磁化率模型条件下，新算法的反演结果迭代5次剩余RMS误差为6.5×10-5；Parker--Oldenburg方法的反演结果迭代20次剩余RMS误差为0.019 2；在应用常磁化率模型条件下，新算法的反演结果迭代3次剩余RMS误差为0.005 2；Parker--Oldenburg方法的反演结果迭代30次剩余RMS误差为0.019 3。

上述对比分析验证了新算法相较于Parker--Oldenburg方法具有更高的精度和效率。

表2 常规Parker--Oldenburg方法与本文所述方法反演参数对比

3 松辽盆地航磁数据居里面计算

松辽盆地是中国重要的含油气盆地，也是中国的一个重要热盆，地热资源丰富。根据磁异常计算区域的居里面深度，为分析区域油气和地热资源，特别是深部干热岩资源具有重要的意义。为了对松辽盆地地热资源进行有效的评估，将本文所述方法应用于松辽盆地居里面深度计算。

3.1 异常分离

图5为研究区化极磁异常。航磁异常的幅值范围在-200～300 nT之间，异常值较弱。磁异常方向呈现NE向和SN向展布规律，局部磁异常呈EW向磁异常带被NE向异常干扰、错断。磁异常方向变化反映了板块构造作用的应力方向及方式的改变，反映了研究区后期主要受太平洋板块俯冲作用的影响[13]。化极磁异常特征分布的复杂性反映了松辽盆地漫长的地质历史中经历了多期次复杂构造运动[14]，磁异常形态差异表明不同刚性程度各微板块在外部应力作用下产生了多样的地质构造形态，这种构造形态具有继承性和叠加性[15]。

图6为前人研究资料获得的松辽盆地基底深度等值线图[16]。为去除松辽盆地上覆地层的磁性影响，其磁性主要来源于泉头组和登娄库组的地层磁性，故选取泉头组和登娄库组地层的平均磁化率75×10-5SI[17]，并结合公式(1)计算得到松辽盆地基底及上覆地层的正演磁异常(图7)，从总磁异常去除松辽盆地基底及上覆地层的正演磁异常可得到如图8所示的剩余磁异常，剩余异常可以看作是以居里面和基底为双界面的模型体所产生的磁异常，为了引入随深度变化的磁化率因子来约束居里面的反演结果，笔者参考前人研究资料中有关松辽盆地的井温--深度数据[13]，通过数据拟合得到了松辽盆地地温随深度的变化。地温随深度变化关系为T=19.334e0.451 6z[13]，其中T为井温，z为深度。大多数磁性物质的磁性随温度呈反比例变化，故以泉头组和登娄库组地层的平均磁化率75×10-5SI为起算点，得到基底与居里面的之间的磁化率随深度的变化关系为κ=3.88×10-5e-0.451 6z。

3.2 松辽盆地居里面计算结果及解释

采用图8的剩余磁异常利用变磁化率模型计算居里面深度如图9所示。反演得到的居里面深度范围为18～25 km，松辽盆地居里面反演结果显示盆地中北部居里面深度最浅，向外依次变深。

图10是松辽盆地的地温梯度分布规律等值线图，松辽盆地中北部地温梯度的变化规律与居里面正好相反，二者的对应关系满足理论上的居里面深度与地温梯度的反比例关系，地温梯度高值区域与居里面深度浅的区域有很好的一致性；此外利用居里点温度和反演得到的居里面深度与计算得到的地温梯度和实测的地温梯度在具体数值精度上也有良好的对应关系。由此推测居里面上隆是松辽盆地地温梯度高的深部原因。

图5 松辽盆地化极磁异常Fig.5 Polarization magnetic anomaly of Songliao Basin

图6 松辽盆地基底深度Fig.6 Basement depth of Songliao Basin

图7 松辽盆地基底及上覆地层正演磁异常Fig.7 Forward magnetic anomalies of basement and overlying strata in Songliao Basin

图8 松辽盆地去除基底及上覆地层原始磁异常后的剩余磁异常Fig.8 Residual magnetic anomalies obtained by removing forward anomalies of basement and overlying strata from original magnetic anomalies

图9 松辽盆地居里面深度反演结果Fig.9 Inverted Curie depth of Songliao Basin

图10 松辽盆地地温梯度图[14]Fig.10 Geothermal gradient map of Songliao Basin

居里面反演结果与通过测井资料获得的地温梯度资料无论是在分布形态还是在数值精度上均有良好的对应关系，通过对居里面反演结果与实测地温梯度之间的对比研究，验证了本文获得的居里面反演结果具有较高的反演精度。

在去除了基底磁性的影响后，剩余磁异常在EW向上的构造特征相比于原始化极磁异常EW向上的构造特征更加明显，表明松辽盆地深部EW向构造被NS向构造所干扰、错断[15,18--19]。

He[20]等根据深部地震Vp/Vs比值和接收函数成像结果，发现在松辽盆地中部地下410 km和660 km深度不连续带存在凹陷，推断深部存在地幔热柱上涌现象。地幔柱的位置与本文计算居里面隆起(图9)吻合，认为松辽盆地中北部区域居里面隆起与地幔热柱上涌相关[21--22]。

松辽盆地中北部基底存在深大断裂构造，为深部热源提供良好的导热通道；松辽盆地深部基底分布有大量花岗岩，存在背斜构造，为热量的传导与储存提供良好的地质条件。居里面隆起的区域是地热包括干热岩重要的有利区，推断图9中居里面深度较浅的松辽盆地中北部区域地热资源丰富，结合松辽盆地地温梯度图(图10)，可以将地热勘探远景区缩小至地温梯度最高值附近，获得图9所示的1、2、3指示的区域。结合前人的研究资料[23--24]，远景区2区域存在大量花岗岩分布，同时处于背斜构造，利于热量的传导和储存；远景区3的花岗岩分布较少，处于向斜构造；远景区1，没有明显的褶皱构造，花岗岩分布最少。通过分析1、2、3区域的花岗岩分布特征和褶皱构造，验证了2区域的地热开发远景最优，3区域其次，1区域最次。因此在进一步勘探中建议优先在2区域中寻找。