基于FCM聚类算法的二维RMT电阻率与介电常数联合反演

2022-06-02易柯张志勇李曼原源郭一豪周峰

地球物理学报 2022年6期

易柯，张志勇,2*，李曼，原源，郭一豪，周峰,3

1 东华理工大学地球物理与测控技术学院，南昌 330013 2 东华理工大学核资源与环境国家重点实验室，南昌 330013 3 中南大学有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083

0 引言

射频大地电磁(Radio-magnetotelluric, RMT)通过接收正交的高频(10kHz～1MHz)电场与磁场分量进行勘探，相较于(音频)大地电磁((A)MT)测深，RMT对百米内的地下介质具有更好的分辨率.该方法在欧洲已应用于地下水(Turberg et al., 1994; Pedersen et al., 2005; Ismail et al., 2011)、环境工程(Tezkan et al., 2000, 2005)、地质构造(Linde and Pedersen, 2004a,b)、地质灾害(Shan et al., 2014, 2016; Wang et al., 2016)等勘探，但在我国RMT方面的研究相对较少(原源等, 2015, 2019).

RMT初期采用与甚低频电阻率(VLF-R)观测方法相似的标量形式(Benson et al., 1997; Turberg et al., 1994)；随后发展了10 kHz～250 kHz频段的张量Enviro-MT仪器(Bastani, 2001)；当前仪器频段得到进一步拓展，研发了10 kHz～1 MHz的张量RMT-F系统(Tezkan and Saraev, 2008)；得益于高频电磁勘探观测时间短的优点，发展了连续观测的RMT系统，方便了水域勘探(Bastani et al., 2015).此外，研发的人工场源射频大地电磁(CSRMT)技术有效提高了观测数据的信噪比(Pedersen et al., 2005; Bastani et al., 2009; Ismail et al., 2011; Saraev et al., 2017).

随RMT观测频段不断提高，位移电流的影响不容忽视.理论计算表明(图1)，对于相对介电常数εr为5的均匀介质，电阻率ρ大于100 Ωm且频率超过100 kHz时，位移电流占比超过1%，电磁波的趋肤深度在高频、高阻时受频率影响变小；随着相对介电常数的增加，位移电流在相对更低频和低阻处即产生较大的影响.高阻介质(ρ≥10000 Ωm)考虑位移电流的二维RMT正演研究表明，受介电常数(εr=5)的影响，视电阻率和相位值从10 kHz开始降低(Kalscheuer et al., 2008)；位移电流受地形影响严重，采用非结构化网格进行带地形正演可有效提高计算精度(原源等, 2015).试算表明对RMT数据同时反演电阻率和介电常数是可行的，但缺少联合约束时，待解参数的增加将直接影响反演结果的精确性(原源等, 2019).此外，RMT与直流联合反演电阻率，提高了资料解释精度(Seher and Tezkan, 2007; Bastani et al., 2012; Yogeshwar et al., 2012; Wang et al., 2018).

图1 趋肤深度和位移电流与传导电流之比左侧为趋肤深度，右侧为位移电流与传导电流之比，(a)、(c)，(b)、(d)对应的相对介电常数εr分别为5、81.Fig.1 Skin depth and the ratio of displacement current to conduction currentLeft two figures represent the skin depth, right two figures represent the ratio of displacement current to conduction current, the relative permittivity of (a) and (c), (b) and (d) are 5 and 81 respectively.

物性相关和结构耦合是构建联合约束的两种主要方法.当不同物性参数存在岩石物理学的耦合关系时可采用物性相关约束(Lees and VanDecar, 1991; Nielsen and Jacobson, 2000)，但该耦合关系通常只适用于局部地区(陈晓等, 2016).采取交叉梯度约束(Gallardo and Meju, 2003, 2004)无需得到参数间的显式解析关系；然而，当物性参数变化无明显方向性时，交叉梯度的约束效果有限(王俊等, 2013).近年来，模糊C均值(Fuzzy C-means, FCM)聚类算法开始应用于联合反演，并取得了良好效果(Paasche and Tronicke, 2007; Carter-Mcauslan et al., 2015; Sun and Li, 2015, 2016).基于正则化理论，FCM聚类算法通过聚类分析得到的隶属度可直接用于反演迭代(Paasche and Tronicke, 2007)；将岩石地球物理先验信息以参考聚类中心的方式融入反演过程(Sun and Li, 2015)，可有效提高反演效果.此外，FCM聚类算法在磁法(Li and Sun, 2016)、重力(Maag and Li, 2018)等数据反演中的应用均取得了良好效果.

目前，尚无基于联合约束的RMT电阻率和介电常数反演研究，为提高联合反演效果，本文在经典最小结构模型正则化约束的基础上，采用FCM聚类算法开展了二维RMT电阻率与介电常数的联合反演研究.研究工作首先验证了二维有限单元正演计算的精度；随后开展了介电常数、电阻率单独反演，进一步验证反演算法各环节的正确性；最后模拟实际工作情况设计了带地形模型，进一步讨论了FCM聚类联合反演的效果以及参数选择对反演结果的影响.

1 全电流二维RMT有限元正演

取时谐因子eiω t时的频率域麦克斯韦方程为

(1)

(2)

二维地电模型见图2，y为地质体走向，z方向垂直向下，磁导率μ取真空磁导率μ0，则电磁场分量与地电参数σ、ε只沿x、z方向变化.

图2 二维地球物理模型Fig.2 Sketch of two dimensional geophysical model

二维条件下，全电流Helmholtz方程为

TE模式：

(3)

TM模式：

(4)

式(3)、(4)可统一表示为偏微分方程：

(5)

综上所述，二维RMT边界条件为

(6)

式(5)、(6)与求解下列变分问题等价：

(7)

对区域进行非结构化三角单元离散后，采用有限单元法对(7)式变分问题进行求解，得以下线性系统：

Bu=0,

(8)

其中B为刚度矩阵，u为节点场值，求解(8)式线性方程，得到各节点主场，并通过主场梯度计算辅助场(徐世浙, 1994)，最后可计算视电阻率和相位响应：

(9)

(10)

Zre、Zim分别为阻抗Z的实部与虚部.

为验证正演计算精度，采用非结构化网格有限单元法对100 Hz～1 MHz频段进行了数值计算，并与一维解析解进行对比分析.均匀半空间和三层水平层状模型对比结果见图3，其中图3a、3b分别为均匀半空间(ρ=10000 Ωm,εr=5)视电阻率和相位响应，一维解析解与二维有限元解曲线吻合，视电阻率平均相对误差为0.4%，相位平均相对误差为0.1%；当频率超过30 kHz后，由于位移电流的影响视电阻率大幅度下降，相位则从10 kHz开始明显下降.图3c、3d分别为三层层状模型视电阻率和相位响应，视电阻率平均相对误差为2.2%，相位平均相对误差为1.1%.综上，正演对比试验验证了开发的二维程序在100 Hz～1 MHz频段计算的正确性，计算精度满足要求.

2 RMT反演理论

2.1 目标函数

根据正则化理论(Tong et al., 2018)构建RMT单独反演介电常数和电阻率的目标函数分别为

P(mε)=φd(mε)+λεφm(mε)

(11)

P(mρ)=φd(mρ)+λρφm(mρ)

(12)

其中，φd为数据拟合泛函，φm为模型稳定泛函，λε、λρ为模型稳定泛函的正则化因子，Wd为数据协方差矩阵，可表示为

其中，s为极小正数，nd为观测数据总个数，i为第i个观测数据的噪音水平.Wε、Wρ为最小结构模型加权矩阵，dobs为观测数据，A(mε,mρ)为正演响应，mε、mρ为待解参数，为参考模型.本文研究中，采用最小二乘法计算非结构化网格模型粗糙度(Lelièvre and Farquharson, 2013)，取参考模型与模型粗糙度之比为1‰(αε 0=αρ 0=0.001)，模型粗糙度可表示为

(13)

其中，αx、αz为比例系数，取αx=αz=1.

图3 100 Hz～1 MHz正演二维有限元解与一维解析解对比(a)、(b) 分别为均匀半空间(ρ=10000 Ωm, εr=5)视电阻率、相位响应; (c)、(d) 分别为三层水平层状模型视电阻率、相位响应.Fig.3 Comparison of 2D FEM solution with 1D analytic solution from 100 Hz to 1 MHz(a) and (b) show the calculated apparent resistivity and phase of homogenous model with ρ=10000 Ωm and εr=5 respectively, (c) and (d) show the calculated apparent resistivity and phase of a three-layered model, respectively.

双参数联合反演目标函数为

(14)

(15)

其中C为聚类中心个数，M为模型单元总个数；mj=(mεj,mρj)T为第j个模型单元物性值，其中mεj、mρj分别为第j个单元的介电常数与电阻率值；vk=(vεk,vρk)T为反演的第k个聚类中心，其中vεk、vρk分别为第k个介电常数与电阻率的聚类中心；ujk为第j个模型单元物性值对第k个聚类中心的隶属度，q为模糊化参数，本文取q=2.0；tk=(tεk,tρk)T为第k个参考聚类中心，tεk、tρk分别为第k个介电常数与电阻率的参考聚类中心；ηk为第k个参考聚类中心tk的权重因子，表示第k个参考聚类中心的置信度.

vk、ujk可表示为(Sun and Li, 2016)

(16)

(17)

聚类项φFCM中第一项可表示为

(18)

结合式(14)、(15)、(18)可将双参数联合反演目标函数表示为

(19)

2.2 目标函数优化求解

对于双参数mε、mρ的反演，采用迭代的高斯-牛顿算法(Zaslavsky et al., 2013)优化求解目标函数(19)的最小值.聚类项中vk和ujk均采用迭代法求解(Sun and Li, 2016)，其第n+1次迭代形式为

(20)

(21)

目标函数(19)的第n+1次迭代形式为

(22)

式(22)中J=(Jε,Jρ)T为灵敏度矩阵，Jε、Jρ分别为介电常数与电阻率的灵敏度矩阵，为减少计算工作量，采用互换定理(Rodi, 1976)进行计算；Δm=(Δmε,Δmρ)T为模型改变量，Δmε、Δmρ分别为介电常数与电阻率的模型改变量.

对目标函数(19)的第n+1次迭代形式求偏导并令其等于0：

(23)

可得高斯-牛顿迭代方程：

(24)

采用双共轭梯度稳定算法(BICGSTAB)求解式(24)可得新的模型参数：

m(n+1)=m(n)+γTΔm(n),

(25)

其中γ=(γε,γρ)T为沿改进量Δm的搜索步长(Zhdanov, 2002)，γε、γρ分别为介电常数与电阻率的搜索步长.

3 模型试算

3.1 电阻率与介电常数反演

为验证算法的正确性，首先进行了相对介电常数异常和电阻率异常单独反演模型试验.设计模型如下：在相对介电常数为5、电阻率为3000 Ωm的均匀介质中，设置长宽分别为400 m、50 m的矩形异常体，异常体上顶埋深20 m.当进行介电常数反演时，取矩形异常体电阻率为3000 Ωm，相对介电常数为80，即只有介电常数异常(图4a)；当进行电阻率反演时，取矩形异常体电阻率为100 Ωm，相对介电常数为5，即只有电阻率异常(图4c).

在10 kHz～400 kHz频段对数等间距取13个频点，点距10 m，共100个测点，采用有限单元法进行正演计算，将视电阻率数据加入5%的随机噪声作为反演数据，进行反演试算.反演取初始正则化因子λε和λρ均为1000，后期根据数据误差、模型误差等信息进行经验选取；反演的初始模型为背景均匀模型，在迭代72次后得到反演的相对介电常数模型图4b与电阻率模型图4d.反演结果表明矩形相对介电常数和电阻率异常体的物性值和轮廓均得到较为准确的恢复，表明本文反演程序正确.

图4 相对介电常数和电阻率单独反演结果(a)、(c) 分别为设计的相对介电常数和电阻率模型; (b)、(d) 分别为反演的相对介电常数和电阻率模型.Fig.4 Separate inverted model of relative permittivity and resistivity(a) and (c) are designed relative permittivity and resistivity model respectively, (b) and (d) are inverted relative permittivity and resistivity model respectively.

3.2 FCM聚类双参数联合反演

为研究FCM聚类对介电常数和电阻率联合约束的效果，设计了带地形的多异常模型见图5，背景电阻率为1000 Ωm，相对介电常数为10，其中加入三个矩形异常体编号(1)、(2)、(3)，其属性见表1.

图5 设计的带地形模型(a) 相对介电常数模型; (b) 电阻率模型.Fig.5 Sketch of designed model with terrain(a) Relative permittivity model; (b) Resistivity model.

表1 地下异常体属性Table 1 Properties of the underground blocks

通过正演计算图5所示模型的响应，测点距、测点数、频点与模型试算3.1节相同.分别讨论参考聚类中心权重因子ηk、聚类项权重因子β、最小结构正则化因子λε与λρ对反演的影响.

首先讨论参考聚类中心权重因子ηk对反演结果的影响，根据位移电流与传导电流的比值至少为1‰，介电常数能产生明显影响，设置λε=1、λρ=1000；聚类项的权重因子β为1固定不变；参考聚类中心为t1[10000 Ωm，3]、t2[10 Ωm，50]、t3[20 Ωm，30]、t4[1000 Ωm，10]，方括号内第1个值为电阻率，第2个值为相对介电常数.在本文算例中认为四个参考聚类中心的置信度相同，后文中将ηk统一表示为η，分别对参考聚类中心权重因子η为1、100、10000的情况进行试算，同时对比不施加FCM聚类联合约束(β=0)的双参数反演结果，试算结果见图6，左侧、右侧分别为相对介电常数与电阻率反演结果，(a)、(e)为无聚类约束的结果，(b)、(f)，(c)、(g)，(d)、(h)分别对应η=1、100、10000.

由图6可见电阻率反演均取得了较好的效果，并且随参考聚类中心权重因子η的增加边界略有改善.不施加联合约束得到的相对介电常数反演模型(图6a)中，(3)号高阻异常体的介电常数得到较好恢复，但两个低阻异常体的介电常数并未恢复，这与低阻异常体内位移电流不足有关；图6b中低阻体的介电常数仍未得到有效恢复；图6c和图6d中低阻体的介电常数得到一定恢复，但边界、物性值与真实模型均有较大差距，总体上图6c和图6d的反演模型类似，通过加大参考聚类中心权重的影响已经到达极限.

根据图6的反演结果可知，基于FCM聚类的电阻率与介电常数联合反演具有可行性.为改善反演效果，进一步讨论聚类项权重因子β、最小结构正则化因子λε与λρ的影响.当η=100时聚类约束对异常的恢复有一定改善，但更大的η对反演结果没有改善，所以取η=100、λε=1、λρ=1000讨论β=10、100时的反演结果，见图7.图7左侧、右侧分别为相对介电常数与电阻率反演结果，(a)、(c)，(b)、(d)分别对应β=10、100.

图6 不同参考聚类中心权重因子η对反演结果的影响(λε=1、λρ=1000，β=1)左侧为相对介电常数反演模型，右侧为电阻率反演模型，其中(a)、(e)为无聚类约束的结果，(b)、(f)，(c)、(g)，(d)、(h)分别对应η=1，η=100，η=104.Fig.6 Test for target cluster centers weight factor η with λε=1, λρ=1000 and β=1Left four figures are inverted relative permittivity models and right four figures are inverted resistivity models, (a) and (e) show the recovered models without FCM clustering constraint, and the test conditions for (b) and (f), (c) and (g), (d) and (h) are η=1, η=100, η=104 respectively.

对比分析图6和图7可以发现，随着β值增大，反演的相对介电常数异常数值有明显改善，当β增大到10，图7a中(1)号低阻异常体的相对介电常数值接近真实值30，但边界未得到明显改善；当β增大到100时(图7b)，相对介电常数异常数值有所改善，但异常体边界变差，而(2)号异常体右下方还出现假异常；电阻率反演模型，β=10时(图7c)高阻异常体电阻率值更接近真值，当β增大到100时(图7d)，(2)号异常体电阻率模型出现了类似于介电常数模型的假异常.

图7 不同FCM聚类约束权重因子β对反演的影响(λε=1、λρ=1000，η=100)左侧为相对介电常数反演模型，右侧为电阻率反演模型，其中(a)、(c)，(b)、(d)分别对应β=10，β=100.Fig.7 Test for FCM clustering constraint weight factor β with λε=1、λρ=1000，η=100Left two figures are inverted relative permittivity models and right two figures are inverted resistivity models, the test conditions for (a) and (c), (b) and (d) are β=10, β=100 respectively.

为分析FCM聚类联合约束效果，绘制了无联合约束反演结果与施加联合约束(β=10，η=100)相对介电常数与电阻率交会图，分别见图8(a、b)，橙色点为真实模型，蓝色点为反演模型.无FCM聚类联合约束(图8a)中反演物性值未形成有效的聚类，而施加联合约束(图8b)反演物性在真实聚类中心附近聚集.

图8 有无FCM聚类约束时反演的相对介电常数和电阻率交会图对比(a) 无聚类约束的反演物性交会图; (b) FCM聚类联合反演物性交会图，橙色点为真实模型物性，蓝色点为反演模型物性.Fig.8 Cross-plot comparison of relative permittivity and resistivity inversion with or without FCM clustering constraint(a) is the cross-plot of inversion without clustering joint constraint, and (b) is the cross-plot of inversion with FCM clustering joint constraint, the orange dots represent true model values and the blue dots represent inverted model values.

为了进一步改进反演效果，在β=10、η=100条件下，继续讨论λε、λρ对反演效果的影响.分别采用λε=10、λρ=104和λε=100、λρ=105两组参数进行反演，反演结果见图9.图9(a、b)为相对介电常数反演模型，图9(c、d)为对应的电阻率反演模型，其中(a、c)，(b、d)分别对应λε=10、λρ=104和λε=100、λρ=105.由图9可见当λε=10、λρ=104时无论电阻率还是介电常数都较λε=100、λρ=105时更好.

图9 不同最小结构正则化因子λε、λρ对反演的影响(β=10，η=100)左侧为相对介电常数反演模型，右侧为电阻率反演模型，其中(a)、(c)，(b)、(d)分别对应λε=10、λρ=104，λε=100、λρ=105.Fig.9 Test for minimum structure regularized factors λε and λρ with β=10, η=100Left two figures are inverted relative permittivity models and right two figures are inverted resistivity models, the test conditions for (a) and (c), (b) and (d) are λε=10、λρ=104, λε=100、λρ=105 respectively.

图6、图7、图9综合对比分析表明，10～400 kHz频段RMT数据基于FCM聚类联合约束的反演算法可以实现介电常数与电阻率双参数反演，尤其是在参考聚类中心tk为真值的前提下，合理选择参考聚类中心权重因子η、聚类项权重因子β、最小结构正则化因子λε与λρ等参数可取得理想反演效果.