基于模糊C均值聚类的二维直流电阻率与射频大地电磁联合反演

张志勇， 易 柯， 谢尚平， 周 峰， 郭一豪， 程 三

（1.核资源与环境国家重点实验室，江西 南昌 330013；2.东华理工大学 地球物理与测控技术学院，江西 南昌 330013）

直流电阻率（direct current resistivity， DCR）与射频大地电磁（radio-magnetotelluric， RMT）是两种浅层地球物理勘探方法。DCR是较早应用的地球物理方法之一，具有较好的浅层勘探效果，理论与应用的发展均较为成熟［1-8］。RMT是近年发展起来的频域电磁方法，其工作频段为10 kHz～1 MHz，适用于地下水［9-11］、环境与工程［12-14］、地质灾害［15-17］等浅地表勘探任务。

研究表明，联合反演相比单一方法可得到更准确的地下模型［18-21］。由于DCR和RMT勘探深度存在重叠，两种方法的联合反演，适用于实测数据处理［22-24］。RMT采用与大地电磁相似的观测方式，易受近地表小尺度异常影响而产生静态位移，通过与DCR数据的联合反演能有效减小该效应影响［25］。此外，在高阻、高频条件下RMT需考虑位移电流，且受地形影响严重［26］，采用非结构网格进行RMT与DCR数据的联合反演有效提高了反演精度与计算效率［27］。相较于单一的DCR或RMT反演，联合反演结合了DCR与RMT对低阻体和高阻体的分辨能力、浅部分辨能力和勘探深度方面的优势，可得到更可靠的地下电阻率模型［28］。

当前DCR与RMT的联合反演研究均采用平滑模型约束，本文为进一步改善电阻率反演效果，将模糊C均值（fuzzy C-means， FCM）聚类方法引入到反演模型约束中。FCM聚类是一种机器学习算法，其原理是根据样本点对所有类的隶属度进行样本的自动分类［29-30］。将具有样本分类能力的FCM聚类应用于模型约束，利用地下地质单元有限分类的实质，将大大提高反演结果与实际地质模型的对应程度。基于FCM聚类的正则化反演，是为了利用不同数据集中得到的模型信息开展地震与重力数据的联合反演［31-32］；另外，在可取得有效地球物理先验信息的情况下，进一步发展了带引导项的FCM聚类联合反演算法［33］，并应用于地震与重力数据的联合反演［34］。此外，FCM聚类算法在磁法［35］、海洋直流电阻率［36］、激发极化［37］、重力［38］、大地电磁［39］等数据反演中的应用均取得了良好效果。然而，在目标函数中引入FCM聚类项势必引起反演目标函数中待确定权重系数的增多，而这些参数的选择对于反演又非常关键［40-41］。

为提高DCR与RMT联合反演在浅层勘探中的应用效果，在经典最小结构模型正则化约束［42］的基础上，引入FCM聚类进行电阻率模型约束，以取得与实际地质模型更接近的反演结果。分别对DCR、RMT单一方法以及联合反演进行了模型试算，检验了算法与程序的正确性，讨论了联合反演的优势；最后分析了FCM聚类约束对提高DCR与RMT联合反演效果的作用。结果表明，联合反演结合了DCR与RMT法的优势，对低、高阻体均实现了高精度反演；FCM聚类约束的应用使得反演的异常体物性值更准确，边界更清晰。

1 DCR与RMT联合反演优势分析

1.1 DCR与RMT正则化反演

DCR与RMT单独反演的正则化目标函数［43］可表示为

式中：m为待解参数；φd(m)为数据拟合项，可统一表示为为数据协方差矩阵，A(m)为正演响应，dobs为反演数据；φm(m)为模型稳定项，可统一表示为为模型加权矩阵，mapr为参考模型；λ为模型稳定项的正则化因子。

本文研究中，DCR与RMT的反演数据分别为视电阻率和阻抗。

Wd的标准形式可表示为

式中：χi为第i个数据的方差；e为极小正数，确保式（2）有效；n为观测数据总个数。

本文采用最小结构模型约束，φm(m)可表示为

采用最小二乘法计算非结构化网格模型粗糙度［44］，取参考模型与模型粗糙度之比为1‰（α0=0.001）。模型粗糙度可表示为

式中：αx、αz为比例系数，通常取αx=αz=1。

1.2 理论反演分析

为讨论DCR与RMT的单独反演能力，设计了图1所示的地下模型。在电阻率为1 000 Ω·m，相对介电常数为40的背景模型中设置编号为（1）、（2）、（3）、（4）的电阻率异常体块，异常体块属性见表1。DCR法采用二极装置，共布设电极124根，电极距为5 m，数据点总数为2 080。RMT法的计算频段为16 384～506 429 Hz，按2的对数等间距取16个频点，测点个数为60，测点距为5 m，参与反演的数据个数为3 840，其中包含阻抗实部与虚部。DCR与RMT测线均位于水平地表，测点位置均关于原点（0，0）对称布设。

表1 异常体块属性Tab.1 Properties of four blocks

图1 设计的计算模型Fig.1 Diagram of synthetic model

采用非结构网格有限单元法进行了DCR与RMT的正演模拟，两种正演数据中各自加入5%的随机噪声，对DCR和RMT数据分别采用高斯-牛顿法［45-46］进行反演，两种方法的正则化因子λ初始值均取1 000，随后每一次迭代的正则化因子为上一次的95%，两种方法均迭代36次，初始模型均为背景电阻率，反演的电阻率模型见图2。其中图2a为DCR数据反演结果，图2b为RMT法横电场（TE）模式数据反演结果，图2c为RMT法横磁场（TM）模式数据反演结果，图2d为RMT法TE、TM模式数据共同反演结果。

图2 DCR和RMT数据单独反演的电阻率模型对比Fig.2 Comparison of resistivity models for DCR and RMT inversions

DCR数据反演（图2a）结果表明，（1）号低阻体对应位置处物性值恢复较好，接近真值，但规模小于实际情况；（2）号高阻体与（3）号低阻体对应位置处物性值与边界均有一定恢复，但（3）号低阻体下边界过于发散，分析是由于当前数据反演深度不足造成；（4）号高阻体对应位置处物性值与边界恢复接近真实模型。RMT法TE模式数据反演（图2b）表明，（1）、（3）号低阻体对应位置处物性值与边界均恢复较好，而（2）、（4）号高阻体对应位置处物性值与边界无明显恢复。RMT法TM模式数据反演（图2c）结果表明，（1）、（3）号低阻体对应位置处物性值与边界均有一定恢复，较TE模式反演得到异常体规模略小；（2）、（4）号高阻体对应位置处恢复仍不理想，但略优于TE模式。RMT法TE与TM模式数据共同反演（图2d）结果表明，（1）、（3）低阻体对应位置处物性值与边界均较单独的TE、TM模式数据反演结果好，（2）、（4）号高阻体对应位置处物性值与边界恢复均略优于单独的TE、TM模式。图3为4次单独反演对应的正则化因子、均方根误差RMS以及模型误差变化曲线，RMS=，两种方法的RMS均呈下降趋势，模型误差均呈上升趋势，反演过程稳定，验证了本文程序设计的正确性。

图3 DCR和RMT数据单独反演的正则化因子、数据均方根误差（RMS）、模型误差变化曲线Fig.3 Curves of regularization factor, RMS, and model object function in single inversions of DCR and RMT data

DCR与RMT单独反演结果表明，DCR反演对高、低阻体均有较好恢复能力；而RMT对低阻体具有高灵敏度，高阻体反演能力则较差。开展DCR与RMT联合反演，一方面加强DCR对低阻体的分辨率；另一方面可避免RMT反演结果遗漏高阻体的风险，实现更精确的浅地表勘探。

2 DCR与RMT的联合反演

2.1 数据比重对联合反演的影响

联合反演的目标函数形式及其优化求解与1.1节单独反演一致，φd中反演数据dobs和正演响应A(m)均由DCR与RMT共同组成，dobs=(ρs，Zyx，Zxy)T，ρs为DCR法视电阻率数据，Zyx、Zxy分别为RMT法的TE、TM模式阻抗数据；A(m)=[ADCR(m)，ARMT(m)]T，ADCR(m)、ARMT(m)分别为DCR与RMT正演响应。对于联合反演，可在标准数据加权矩阵Wd中引入平衡算子调节数据比重［28］，本文采用的形式为

式中：x代表DCR与RMT方法；Wdx代表各类型数据的加权矩阵包含代表各类型数据的平衡算子，δx用于表征数据集对反演模型参数的贡献比例，由于DCR与RMT数据均只反演电阻率，因此取δRMT=δDCR=1.0；Nx代表各类型数据个数，NRMT与NDCR分别代表RMT与DCR数据个数。

对图1所示模型数据，在正则化因子、迭代次数等参数与1.2节单独反演一致的条件下，对基于标准数据加权的Wd以及引入平衡算子且δDCR=δRMT=1.0的Wd分别进行联合反演试算，结果分别见图4a、4b；图5为对应的RMSDCR、RMSRMT变化曲线，图5a对应图4a，图5b对应图4b，虚线表示RMSDCR或RMSRMT的值为1.0。

对比图4a与图2d可知，图4a联合反演结果中高、低阻体对应位置处物性值与边界恢复情况较图2d RMT反演结果均较好，其中（4）号高阻体的物性值由图2d的3 005 Ω·m上升至3 965 Ω·m；对比图4a与图2a可知，图4a联合反演的结果中高阻体对应位置处物性值与边界恢复较图2a DCR反演稍差，但低阻体对应位置处的恢复明显好于图2a DCR反演结果，图4a对应的RMSDCR=1.05，RMSRMT=0.88；当在数据加权矩阵Wd中引入平衡算子后（图4b），（4）号高阻体对应位置的物性值与边界恢复明显提升，接近图2a DCR反演结果，对应的RMSDCR=0.99，RMSRMT=0.90。分析图5的RMS曲线可知，图4a中高阻体反演不理想与DCR数据拟合不足有关，在第20次迭代后图5b中RMSDCR较图5a更好收敛。

图5 联合反演的RMSDCR、RMSRMT变化曲线Fig.5 RMSDCR versus RMSRMT in joint inversions

对比图2与图4，DCR与RMT的联合反演结合了两种方法的优势，对高、低阻体均有较好的反演能力，相较单一方法得到了更准确的地下电阻率模型。在实际反演中对不同数据的方差假设往往不一致，可能导致不同数据的权重差距过大，可通过引入平衡算子，改变式（5）中的δx以调节不同数据的比重，进而改善联合反演效果。

图4 基于不同数据加权矩阵的联合反演电阻率模型对比Fig.4 Comparison of resistivity models for joint inversions based on different data weightings

2.2 基于FCM聚类的联合反演

2.2.1 基于FCM聚类模型约束的DCR与RMT联合反演

研究表明，在反演中引入FCM聚类约束可得到更好的地质分异信息［31，34］，有利于提高地质解释的精度。为进一步提高DCR与RMT数据联合反演效果，引入FCM聚类进行电阻率模型约束。

基于FCM聚类的DCR与RMT联合反演目标函数可表示为

式中：φFCM(m)为FCM聚类项［34］；β为聚类项的权重因子。

聚类项可表示为

式中：M为总模型单元个数；C为聚类中心个数；mj为第j个模型单元的物性值；vk为反演的第k个聚类中心；uqjk为第j个模型单元物性值对第k个聚类中心的隶属度，其中q为模糊化参数，本文取q=2；在获得岩石物理先验信息时，可引入参考聚类中心tk；ηk为第k个参考聚类中心tk的权重因子，表示第k个参考聚类中心的置信度。

基于FCM聚类的DCR与RMT联合反演目标函数可进一步表示为

聚类项中第一项可表示为

采用高斯-牛顿法优化求解式（8）目标函数的最小值，第n次迭代的高斯-牛顿方程为

式中：J为雅可比矩阵。根据式（10）可得新的模型参数为

式中：γ为沿改进量Δm的搜索步长［43］。

2.2.2 FCM聚类模型约束权重的自动选择

聚类项在目标函数中的权重可直接影响聚类反演效果。理论上，在反演初期，当采用均匀模型反演时，无法施行有效的聚类，φFCM(m)的权重β过大将影响反演的正常进行；而随反演的进行，异常信息将逐渐清晰，分类特征将越来越明显，应当有较大的聚类权重。根据以上分析，设计随迭代过程自动调整β值的方法为

2.2.3 理论模型反演分析

为了分析聚类约束对联合反演效果的改进，进行了联合反演与基于FCM聚类的联合反演试算。正则化因子初始值均为10 000，后期根据数据误差以及模型误差等信息进行经验选取，迭代次数均为60次，数据加权方式均采取2.1节引入平衡因子的形式。对于聚类联合反演，取L=10、100、500、 1 000，根据式（12）的方法选取β值；为简化参数讨论，文本根据真实模型的物性值设置3个参考聚类中心，t1=100、t2=1 000、t3=10 000，并将其置信度η均设置为50；将RMS=1.0或达到最大迭代次数作为反演终止条件，确保2.2.2节中β自动选择方案有效。基于以上条件进行模型反演试算，结果见图6与图7。图6为不带聚类约束得到的联合反演结果，图7为FCM聚类联合反演结果，图7a、7b、7c、7d分别对应L=10、100、500、1000。图8为FCM聚类联合反演对应的迭代曲线，图8a为RMS变化情况以及根据RMS计算的β值，RMS值均呈下降趋势并趋于稳定。图8b为FCM聚类项的模型稳定函数φFCM变化值，随着FCM聚类项权重β的增大而增大。

图6 不带聚类约束的联合反演电阻率模型Fig.6 Resistivity model of joint inversion without FCM clustering

图8 不同L取值时聚类联合反演的RMS、β、φFCM的变化曲线Fig.8 Curves of RMS,β,φFCM at iteration of different L values based on FCM

对比图6与图7，随着L的增大即聚类约束的增强，联合反演得到的4个异常体的物性值均更接近真实模型，边界均更清晰，但并不是L值越大异常体边界与实际模型越吻合，反而L=500、1 000时边界差异更大；L=100时反演的电阻率模型物性值与边界最接近真实模型。当前，尚未开展β与L的自适应选择方案研究，L值表征聚类约束强弱，通过试算选择。

图7 不同L取值的FCM聚类联合反演电阻率模型对比Fig.7 Resistivity models for joint inversion based on FCM at different L values

为分析反演早期聚类约束过强对联合反演的影响，对β取固定值与自动选取β两种情况的反演过程进行对比。L=100时反演结果最好，将该反演最后一次迭代的β值5.817作为β的固定值，将其反演结果与L=100自动选择β时的反演结果进行对比，见图9。其中，图9a、9b、9c为β取固定值5.817时的反演结果，图9d、9e、9f为L=100自动选择β时的反演结果，9a、9d，9b、9e，9c、9f对应的迭代次数分别为第18次，第36次与第60次。

对比图9a、9b、9c与图9d、9e、9f可以看出，当β取固定值5.817时，反演初期过大的聚类约束减小了反演向真实模型改进的动力，导致反演的电阻率模型较差，进而影响最终的反演结果，特别是高阻体；而当采取自动选择聚类权重β的策略时，随着迭代次数的增大，反演的电阻率模型越来越接近真实模型。

图9 β取定值5.817时的FCM聚类联合反演与L=100自动选择β的反演结果对比Fig.9 Resistivity models for joint inversions based on FCM at two different β values (one being a fixed β value of 5.817, the other being an automatic selection value of β), and an L of 100

对比图6、7、9可以看出，FCM聚类约束可一定程度上提高DCR与RMT的联合反演效果，对于聚类项权重因子β的选择，可根据数据均方根误差在反演初期较大、后期较小的性质，通过数据均方根误差自动计算得到β值，以避免反演初期聚类约束过强对反演结果造成的不良影响。

3 结论与展望

通过对二维DCR、RMT数据的正则化反演与基于FCM聚类模型约束的联合反演研究，取得以下成果：

（1） 对于二维DCR与RMT数据的单独反演，DCR方法在其勘探深度内对高、低阻体均较敏感，而RMT方法对低阻体分辨能力明显优于高阻体。

（2） DCR与RMT方法的联合反演较单一方法可以得到更准确的地下电阻率模型；联合反演中，DCR与RMT数据在反演中的权重对反演结果影响较大，通过平衡两个数据集在反演中的比例，可有效提高联合反演效果。

（3） FCM聚类约束可改善联合反演效果；反演初期聚类项权重过大将影响反演目标函数优化过程中的下降动力，导致最终的反演结果较差；以数据均方根误差为依据的聚类项权重因子自动选取算法适用于DCR与RMT数据的联合反演。

RMT数据受位移电流的影响，可对介电常数进行反演，后续将进行DCR与RMT数据的电阻率、介电常数联合反演研究。鉴于我国RMT方法研究处于理论研究阶段，课题组正与国外团队合作，进行RMT实测数据反演研究。

作者贡献声明：

张志勇：方法提出，算法设计与代码撰写，论文撰写与修改。

易 柯：算法改进，论文撰写与修改。

谢尚平：试算，论文修改。

周 峰：方法讨论，论文修改。

郭一豪：程序改进，论文修改。

程 三：数值模拟，论文修改。