基于模糊C均值聚类约束的井-地磁法联合数据空间反演

2021-07-16徐凯军季春晖梁舒瑗

中国石油大学学报（自然科学版） 2021年3期

徐凯军, 李猛, 季春晖, 梁舒瑗

(中国石油大学(华东)地球科学与技术学院,山东青岛 266580)

目前深部矿产资源勘探已成为中国勘探重点[1]。井中磁测对深部矿体的磁场变化非常敏感,在深部金属矿勘探中发挥着重要作用[2-4]。考虑到井的数量有限,因此结合地面磁测数据,进行井-地磁异常联合反演,可以有效提高反演结果的空间分辨率[5-7]。三维磁异常反演计算量大,为了提高地球物理数据反演计算效率,一些学者开展了数据空间地球物理反演算法,极大地减少反演时间[8-11]。加入先验信息约束可以提高反演可靠性[12],基于模糊C均值聚类(FCM)的约束反演,可以将先验地球物理信息加入到反演中,对反演结果进行聚类,使反演的矿体边界清晰,反演结果更符合地质解释[13-19]。井中磁测不仅可以观测不同深度的三分量磁异常,同时钻井岩心、磁化率测井曲线等信息能够提供较精确的地下物性信息,更适合进行FCM约束反演,提高反演可靠性。笔者将数据空间反演方法和FCM约束反演方法相结合,实现基于模糊C均值聚类约束的井地磁异常数据空间联合反演。

1 井地磁异常正反演理论

1.1 井中磁异常正演方法

磁法正演通常将地下剖分为一系列长方体网格，采用不规则网格，矿体区域网格剖分较为细密。井中观测点接收到的磁场等于所有地下网格在观测点所引起的磁场总合。观测点处产生的磁场三分量异常表达式如下：

(1)

其中

式中，Bx、By和Bz分别为x、y、z三个方向的磁异常；μ0为真空磁导率；Mx、My和Mz分别为x、y、z三个方向的磁化强度。

1.2 井-地磁异常联合反演方法

1.2.1 数据空间反演算法

三维磁法模型空间反演通常计算量大,数据空间反演是基于模型空间反演算法的一种改进,减少反演时间,提高了反演效率。

磁法数据反演的目标函数通常为

(2)

式中,φ(m)为目标函数;φd(m)为数据拟合目标函数;φm(m)为模型目标函数;d为观测数据向量;F(m)为正演数据向量;Wd为数据对角加权矩阵;m为模型向量;Wm为模型加权矩阵;m0为参考模型向量。反演计算就是对目标函数求解极小值的过程。因此对目标函数求偏导,并令其等于0,可得

(3)

式中,J为雅克比矩阵。基于式(3)求解,可以得到反演迭代公式,常规磁法数据反演都是基于模型空间,其反演迭代公式为

m=m0+(JTWdJ+Wm)-1JTWd(d-F(m)+J(m-m0)).

(4)

式中,JTWdJ+Wm为M×M阶矩阵;M为模型空间个数。

在模型空间反演迭代公式基础上进行转换,得到数据空间反演迭代公式[8]:

(5)

1.2.2 模糊C均值聚类约束反演

模糊C均值聚类算法通过指定聚类中心的个数,自动计算出每个样本数据相对于聚类中心的“亲昵程度”(隶属度),从而确定样本点的类别,完成数据分类的目的。FCM的数学表达式[13]为

(6)

其中

式中,M为需要聚类数据的个数;C为聚类中心的个数;pj为第j个数据;vk为第k个聚类中心值;μjk为第j个数据pj相对第k个聚类中心vk的隶属度;q为隶属度μjk的模糊系数,取值范围为[1,∞)。当q=1时,聚类结果会产生类似“硬聚类”的效果,即每个数据对象仅分配给一个聚类中心,q值越大时聚类效果越模糊。

Sun等[13]在FCM算法的基础上加入了聚类中心参考项,改进后的FCM算法为

(7)

其中

式中,tk为第k个聚类中心的参考值;ηk为评价第k个聚类中心是否有已知参考值的系数。ηk的取值范围为0或1,当第k个聚类中心有已知参考值时,ηk=1;当第k个聚类中心没有已知参考值时,ηk=0。

mc=vμ.

(8)

加入FCM约束的反演方法在拟合观测数据的同时对模型的物性做出了约束限制,因此反演结果即符合观测到的地球物理场又满足先验的统计岩石物理信息,同时FCM聚类方法使反演异常体边界更加清晰。

1.2.3 井-地磁异常联合反演

常规地面磁测数据具有较高的水平分辨率,但纵向分辨率较差,而井中磁测数据具有较高的纵向分辨率,但受井数量限制,控制范围有限。充分的利用地面磁测和井中磁测的优点,进行井-地磁异常联合反演,能够提高磁异常反演结果的空间分辨率。

对于井-地联合反演,式(2)中φd可写为

(9)

式中,γ为权重因子;γ1后的部分为地面磁测数据的残差;γ2→γn后的部分为不同井中磁测数据的残差。

联合反演中要使井中磁测和地面磁测的最大振幅保持在同一个数量级内,当井中磁测的振幅远大于地面磁测振幅时,适当降低井中磁测数据的权重,使反演结果不会因为过于拟合井中磁测数据而造成地面磁测数据的欠拟合的问题。

2 井-地磁异常特征分析

2.1 井中三分量磁异常特征

地下为单一矿体,长、宽均为200 m,高为100 m,矿体顶面埋深为400 m。布设4个井,其中井L1穿过矿体,井L2～L4未穿过矿体,模型与井位分布如图1(a)所示。矿体的磁化强度为1 A/m,磁倾角为90°,磁偏角为0°,地磁场倾角为90°,地磁偏角为0°。井中测点从z=0起,相隔20 m设一个测点,正演计算井中测点的三分量磁异常。图1(b)为不同井的三分量磁异常曲线,图中灰色部分显示了矿体的深度范围,井位于矿体中心线上(x=500 m),由于对称性,磁异常Bx分量为0。可以看出,当井穿过矿体时,磁异常Bz分量在矿体位置出现显著变化,变化极值点可以确定磁性体的上下边界。当井没有穿过矿体时,磁异常垂直分量Bz的峰值基本对应矿体中心位置,磁异常水平分量By正负峰值基本对应矿体的上下边界,随着井离矿体越远,磁异常幅值越小。

图1 模型及井中三分量磁异常曲线

2.2 井-地磁异常对比分析

设置地下同时存在3个矿体的模型,模型和井分布如图2(a)所示。矿体的磁化强度为1 A/m,磁倾角为90°,磁偏角为0°,地磁场倾角为90°,地磁偏角为0°。图2(b)为井中三分量磁异常曲线图。可以看出,磁异常具有很好地纵向分辨率,能够识别出离井较近的两个纵向叠置矿体不同位置,而矿体3由于离井较远,异常较弱,井中三分量磁异常没有明显响应。图2(c)为该模型的地面磁异常等值线图,可以看到两个明显的异常区域,但无法识别纵向叠置的异常体。因此井中磁测具有较高的纵向分辨率,但受限于测井位置,横向控制范围有限,而地面磁异常具有较好的横向分辨率,但纵向分辨率差,因此结合两者的优势开展井-地磁异常联合反演是提高反演分辨率的有效方式。

图2 模型及井-地磁异常

3 模型反演

3.1 数据空间反演与模型空间反演对比

建立4个顶面埋深在300 m,尺寸为500 m×500 m×480 m的磁性体,如图3(a)所示,磁化强度为2 A/m,磁化倾角为90°,磁偏角为0°,背景区域无磁性,地磁倾角为90°,地磁偏角为0°。异常体区域剖分成50×50×40个网格单元,每个网格的尺寸为50 m×50 m×40 m。地面上从坐标原点起,在x和y方向上隔50 m设置一个测点,观测点总数为2 500个。图3(b)为正演计算得到的地面磁异常等值线图。基于普通微机分别开展数据空间反演和模型空间反演,分析反演的计算效率。基于模型空间反演耗时为50 min,而基于数据空间反演耗时为15 min,可以看到基于数据空间反演明显减少了反演所用的时间。

图3 模型及地面磁异常

图4分别展示了模型空间和数据空间三维反演结果示意图,可以看出模型空间反演和数据空间反演都能反演出矿体的空间位置,但反演的物性值与真实值2 A/m都有较大偏差。

图4 反演结果三维示意图

3.2 地面磁测反演与井-地磁测联合反演对比

设计两个倾斜脉状体共同组合成的模型,模型中设置两口井,井位坐标分别为L1(650,1 000)、L2(1 500,1 000),如图5(a)所示。矿体在深度剖面上呈现“y”形状,背景区域的磁性为零,矿体的磁化强度为2 A/m,磁化倾角75°,磁偏角25°。图5(b)为正演得到的地面磁异常等值线图,地面磁异常场的总体分布特征为中间高四周低,磁异常存在着两个高值区间。基于正演获得井中三分量磁异常,从z=0起间隔50 m设置一个测点,图5(c)显示了井中磁测三分量曲线。井L1和L2的磁异常Bz分量表现了明显的穿矿特征。

图5 模型及井-地磁异常图

针对上述模型仅使用地面磁异常数据进行反演,反演结果如图6(a)所示,图中仅展示磁化强度高于0.4 A/m的反演结果,以展示反演结果的整体特征。反演结果显示深部异常体形态与实际矿体有较大差异,说明地面磁测纵向分辨率差。反演结果的磁化强度极大值为1.4 A/m,与实际矿体的磁化强度值2 A/m也有较大的差异。

图6 模型反演结果

联合地面磁测和井中三分量磁异常进行联合反演,反演结果如图6(b)所示。可以看出,相比于单独地面磁异常反演,井-地磁异常联合反演结果更接近于真实矿体模型,“y”型矿体得到很好地恢复,说明井-地磁异常联合反演能够结合地面磁测横向分辨率高和井中磁测纵向分辨率高的优点,提高反演结果的空间分辨率,但反演结果的磁化强度极值为2.4 A/m,与实际矿体的物性值仍存在一定的差异。

3.3 FCM约束反演

针对图5(a)模型,设置不同的聚类参数,分析不同聚类参数时的反演效果。真实模型的聚类中心为2,地质背景磁化强度为0,矿体磁化强度为2 A/m。设置二组不同聚类参数,第一组聚类中心数C为2,聚类中心参考值仅输入一个0值(即反演中背景层的磁性已知,而矿体的磁性未知);第二组聚类中心数C为4,聚类中心参考值仅输入0和2 A/m(反演中猜测地下存在4种磁性体,但仅已知两种岩石的物性值)。

分别使用这二组参数开展基于FCM约束的井-地磁异常联合反演。第一组的反演结果如图7所示。可以看出,基于FCM约束反演的物性分布直方图更接近于真实模型。反演异常体的边界清晰,而且反演物性值与真实模型较为一致。

图7 第一组参数的FCM约束反演结果

第二组参数的反演结果如图8所示,由于聚类中心数C设为4,反演结果物性分布直方图也显示了磁化强度主要集中分布在0、0.4、1和2 A/m 4个区间内,反演异常体边界变得模糊,与真实模型有了明显差异。

图8 第二组参数的FCM约束反演结果

通过不同聚类参数的反演结果对比分析可以看出,FCM约束反演的效果与聚类中心数C的选取有关,C取值越接近真实情况,反演效果越好。在使用FCM约束反演方法时,聚类中心数C的取值应尽量与地下岩石的物性种类数目一致。钻井和测井资料能够提供较为精确的先验物性信息,这也为FCM约束反演提供了好的基础。

3.4 冬瓜山铜矿建模及反演

冬瓜山铜矿位于安徽铜陵地区,已做过不同程度的矿产调查工作,该区的地质构造情况已有较为详细的认知。以冬瓜山铜矿的地质模型(图9)为基础,调研了相关矿物的磁化率分布区间[20-21],然后建立冬瓜山地区的三维磁化强度模型,如图10所示。矿体的磁化角度按当地的地磁场角度设置,即磁化倾角47.4°,磁化偏角-5.54°。

图9 冬瓜山铜矿典型地质模型[20]

图10 冬瓜山铜矿三维磁化强度模型

在该模型中设置3口井,井位坐标分别为L1(700,500),L2(300,600),L3(1 000,1 000),井深1 000 m。地面测点间隔100 m,井中测点间隔40 m,通过正演计算分别得到地面磁力异常和井中三分量磁力异常,然后加入5%的高斯噪音作为观测数据。图11显示了地面磁异常和井中三分量磁异常。

图11 地面和井中磁异常

分别进行地面磁异常反演、井-地磁异常联合反演和基于FCM约束的井-地磁异常联合反演,分析不同反演的效果。图12显示了3种方法的反演效果。

图12 不同反演方法结果对比

可以看出,仅使用地面磁异常反演效果较差,铜矿脉的埋藏较深且上方存在磁性岩体遮挡,地面磁异常反演无法从中分辨出不同的磁性矿体,更无法确定出埋藏最深的铜矿脉的位置形态。井地磁异常联合反演提高了矿体的空间分辨率,上、中、下3部分的矿体在反演结果里都得以体现,但矿脉的连续性较差,而且反演的磁化强度数值和真实值有一定的误差。基于FCM约束的井地磁异常反演效果得到了进一步提升,反演的矿体边界较为清晰,磁化强度值更接近于真实的矿体物性值。