2020, Vol. 63
2. 东华理工大学地球物理与测控技术学院, 南昌 330013;
3. 中国矿业大学资源与地球科学学院, 徐州 221116
, CHEN Xiao1,2
, YANG HaiYan2,3, ZHANG ZhiYong1,2, GUO YiHao1,2, LIU Xing1,2, YE YiXing1,2
2. School of Geophysics and Measurement-Control Technology, East China University of Technology, Nanchang 330013, China;
3. School of Resources and Geosciences, China University of Mining and Technology, Xuzhou 221116, China
地球物理联合反演可以促进不同的地球物理方法互为约束，降低单一方法的多解性，是地球物理反演领域的重要发展方向(杨辉等，2002；于鹏等，2006；陈洁等，2007).现阶段联合反演的实现方式总体可以分为两类：一类是基于岩石物性约束(Heincke et al., 2006；Colombo and De Stefano，2007)；另一类是基于构造几何约束(Haber and Oldenburg, 1997；Gallardo and Meju, 2003；彭淼等，2013).岩石物性约束型联合反演是通过岩石物性参数之间所满足的经验或统计性物性关联来实现岩石物性参数的耦合.例如密度和速度之间的Gardner公式，是常见的经验性岩石物性关联，且已应用于地球物理联合反演中(Liu et al., 2011).由于经验性物性关联的适用性有限，学者们通过分析研究区域实测的岩石物性参数的特征进而建立起统计性物性关系，是常见的岩石物性关联确定方式.如Dell′Aversana(2001)根据测井资料电阻率与速度交绘图，采用双对数的函数关系建立了电阻率和速度之间的物性关联，并将其应用于地震、MT和重力数据的综合解释当中；Jegen等(2009)根据大洋钻探(ODP)的资料统计出速度与电阻率、密度的关系，在地震信息的约束下进行联合反演，获取了玄武岩下伏地层的地质构造信息；Moorkamp等(2011)将电导率与慢度、重力与慢度的函数关系用于重电震三维联合反演.此外，当岩石物性关联统计特征复杂时，通过某些中间参数为桥梁，如:孔隙度、流体饱和度等(Maier et al., 2010；Dell′Aversan et al., 2016)，确定物性参数之间的关联特征也是可行的，然而这种岩石物性关联显然不易建立.近些年来，模糊c均值(FCM)聚类算法也逐渐应用于岩石物性约束型联合反演(Sun and Li, 2011).Carter-McAuslan等(2015)对此相关技术进行了详细介绍，并且针对先验信息精准和不精准的情况，分析了模糊c均值聚类联合反演方法的灵敏度.但模糊聚类方法需要将聚类约束项引入目标函数一起进行极小化处理，无疑增加了联合反演的实现难度.而且利用模糊聚类所确定的岩石物性关联也具有局限性，可能只适用于某种地质条件(Colombo and Rovetta, 2018).需要指出的是，岩石物性关联的有效建立决定了联合反演的可靠性，但是对于某些地球物理属性，如地震速度和电阻率，很难建立起二者之间普遍适用的经验关系(殷长春等，2018).针对岩石物性关联信息不易建立、适用性有限的问题，陈晓等(2016, 2017)以基于模拟退火算法的MT和地震联合反演为例提出了宽范围物性约束技术.模型试验表明，该技术容易实现，具有一定的容错性，可降低联合反演对先验信息的要求.张磊(2016)在实现了基于模拟退火算法的MT与重力宽范围岩石物性约束联合反演的基础上，提出了随机正反比的宽范围岩石物性约束方案，增加了其应用的灵活性.郭曼(2018)将宽范围物性约束技术引入到基于模拟退火算法的MT和重力贝叶斯联合反演当中，进一步扩展了该技术的应用范围.
差分进化算法(DE)由Storn和Price(1997)提出，作为一种全局寻优算法，其核心是利用个体间的差分信息作为扰动量，通过变异和交叉操作，使在得到的新变异个体中又能保留上代个体的优势基因，并经过贪婪算法选择优势个体，最终逐渐使种群向最优解位置靠拢.DE算法采用实数编码，操作简单，鲁棒性好，在地球物理反演领域已有应用.如Wang和Gao(2012)针对高维波形反演问题，将地球物理模型参数分解成若干部分，每个部分引入局部适应度函数来设计新的变异算子进而指导变异方向，并实现了不同部分之间的协同进化；宋维琪等(2013)将DE算法加入到微地震贝叶斯反演，提高了其反演定位的精度；Gao等(2019)在多变异DE算法(MMDE)的基础上结合深度学习技术，提出了一种新的可开展训练学习的深度网络，并将其应用于地震阻抗反演中.但目前来看，关于DE算法的地球物理反演文献依旧较少，只有少量用于二维问题当中.且与经典的全局寻优算法(如模拟退火，遗传算法等)相比，DE算法在地球物理联合反演应用更加少见，也还未见有文献利用DE算法实现MT和重力的联合反演.在DE算法改进方面，种群结构的改进被众多学者所关注.吴亮红等(2007)提出一种双种群DE算法，结合不同变异策略的优势，并用平均熵来初始化种群，结果表明该算法全局搜索能力强、鲁棒性好、收敛速度快.孟红云等(2008)提出一种基于双群体搜索机制的DE算法，用于求解约束多目标的优化问题.谭跃等(2010)采用双种群结构，使其中一个种群所有个体在每代DE算法后，再对这些个体实施交叉和变异策略以增加其寻优能力.暴励和曾建潮(2011)在双种群结构DE算法中，加入人工蜂群算法作为另一种群进化方法，再将两种群进行交互学习，提高算法的收敛速度.王天意(2015)在双种群DE算法中，将其一种群的变异算子和交叉算子取定值，另一种群则随迭代次数自适应调整，实现了大地电磁反演.Wang等(2018)在双种群突变策略基础上，加入个体适应的选择机制和归档技术进行检测保护停滞不前和趋同的个体，提高了多目标全局最优解的寻优能力和精度.双种群结构的改进DE算法具有很强的灵活性，可以结合不同变异策略特点加快算法收敛速度，也可以将差分进化算法结合其他智能非线性算法，进而促使种群间协同进化，提高算法寻优能力.
可以看出，相对于传统的点对点映射式岩石物性约束联合反演(即利用一个经验或统计性物性关联的公式将一种岩石物性模型直接映射到另外一种岩石物性参数模型)，宽范围物性约束技术是一种新型的物性约束技术，易于实现且容错率高.当岩石物性关联复杂甚至难以建立时，该技术在一定程度上可以减小对先验信息的依赖.但目前该技术只是结合模拟退火算法实现的，该技术能否结合其他优化算法实现岩石物性参数在一定“范围”内耦合还未得到验证.并且该技术结合其他算法的实现，可以进一步验证宽范围物性约束思想在联合反演领域的适用性，减小岩石物性约束方式联合反演的局限性，增强其灵活性.此外，DE算法在地球物理联合反演领域应用不足，而双种群设置方案是对DE算法进行改进和完善的一条重要思路.基于此，本文以双种群设置为基本框架，改进DE算法，并以MT和重力联合反演为例，研发基于改进DE算法的宽范围物性约束技术.
1 目标函数根据正则化反演理论，MT和重力正则化联合反演通常由数据拟合项和模型约束项组成:
|
(1) |
式中，Pα(m)表示联合反演目标函数，dobsMT、dobsG和dcalMT、dcalG分别是MT和重力的实际观测数据和理论计算数据，WMT和WG分别表示MT和重力的加权因子，αMT和αG分别是MT和重力的正则化因子，▽mMT和▽mG分别对应MT和重力的模型参数的梯度.本文采取双正则化因子自适应调整方案(张磊，2016)，分别选择初始种群结果中MT与重力的数据拟合项及模型约束项之比作为各自初始正则化因子，当MT与重力拟合效率下降比例低于某一定值则对各自的正则化因子进行自适应衰减.MT正演采用有限差分法，重力正演采用网格累加法.
2 基于改进的DE算法的宽范围物性约束技术 2.1 改进的双种群DE算法本文以吴亮红(2007)提出的双种群DE算法为基本框架，在变异因子的自适应确定以及局部寻优能力的提高方面进行了改进.改进的双种群DE算法由双种群的生成、变异操作、交叉操作、选择操作以及混沌搜索操作组成：
(1) 双种群的生成.根据种群数量，将种群分为种群1和种群2，两个种群随机生成初始个体xi1.
(2) 变异操作.文献(Storn and Price, 1997)显示(2)式的操作全局寻优能力强，但收敛速度慢，(3)式的操作收敛速度快，但全局寻优能力差.为了兼顾收敛速度和寻优能力，在新算法的变异操作环节，种群1采用公式(2)，种群2采用公式(3)：
|
(2) |
|
(3) |
其中，νit+1表示变异后的个体，下标i≠p1≠p2≠p3≠p4表示变异操作需要利用互不相同的个体，F是变异因子.此外，在迭代寻优的过程中，每隔一定代数令两个种群中的个体按照一定概率进行随机交换进而保证两个种群的独立性和交互性.
为了将个体的变异算子与个体的适应度相联系，本文借鉴前人的思路(McGinley et al., 2011)，利用(4)式自适应确定变异算子：
|
(4) |
其中，[Fmin，Fmax]是变异因子的取值范围，在公式(2)、(3)中为[0，1].f(xi)为种群中第i个体的目标函数值，f(xbest)和f(xworst)分别为当代种群中最优个体和最差个体的目标函数值.λ是一小数，是为了保证式(4)分母不会为零，本文取λ=10-13，避免对结果造成较大影响.
(3) 交叉操作.本文按照(5)式将得到的变异个体和上代个体进行交叉操作：
|
(5) |
其中Uidt+1表示交叉后的个体，CR是交叉因子，取值范围为[0, 1].
(4) 选择操作.两个种群根据(6)式所示的贪婪选择机制来判断保留的下一代个体是新更新的个体Uidt+1还是父代个体xidt.最后，在两个种群完成DE算法寻优之后，比较两个种群的最优个体best1和best2对应的适应度，确定两个种群的最优个体best12.
|
(6) |
(5) 混沌搜索操作.为了提高DE算法的局部寻优能力，本文将混沌算法引入DE算法.利用混沌搜索方式在步骤(4)得到的最优个体best12附近进行再搜索，具体如下形式：
|
(7) |
其中
针对目前岩石物性约束方式存在的问题，我们提出了宽范围物性约束技术，该技术的基本思想是将岩石物性先验信息与所采用的优化算法相结合，不再简单地将先验物性关联映射获得的模型直接代入联合反演运算，而是在一定范围内进行再搜索，该技术旨在发挥先验物性信息的导向作用同时，又充分利用优化算法的寻优能力.
陈晓等(2016, 2017)将宽范围物性约束思想融入模拟退火算法的模型扰动环节.结合上述基本思想，并基于2.1节所提出的改进的DE算法，本文尝试在DE算法的变异操作和混沌搜索操作环节融入宽范围物性约束思想.
2.2.2 基于改进DE算法的宽范围物性约束技术本文以MT和重力联合反演为例，通过电阻率和密度的宽范围耦合细节，来介绍基于改进DE算法的宽范围物性约束技术的流程：
(1) 生成初始种群.设置最大种群数量，最大迭代次数，并将初始种群分为种群1和种群2.
(2) 变异操作.种群1的电阻率模型按照(2)式进行变异，种群2的电阻率模型按照(3)式进行变异.不同于传统的岩石物性约束联合反演利用经验性或统计性物性关联直接映射出地球物理模型，种群1的密度模型按照(8)式变异，种群2的密度模型按照(9)式变异.
|
(8) |
|
(9) |
式中，vt+1, iG表示第t+1代种群中第i个密度模型个体.f(xt, iMT)为第t代种群的第i个体的参考密度模型，该模型是由传统的物性模型转换获得，即利用先验物性关联由步骤(2)中的电阻率模型直接映射所得，xavgG表示为截至上次迭代，每一代密度模型的最优解的平均模型，F1为变异算子，F2为一随机数.可以看出，(8)和(9)式右边有三项组成，f(xt, iMT)的引入旨在保证先验物性信息可以发挥导向作用，xt, p1G-xt, p2G和xt, p1G-xt, p2G+xt, p3G-xt, p4G这两个不同搜索项的引入旨在保证两个种群的独立性，xavgG的引入旨在保证寻优的效率.
(3) 交叉操作.电阻率和密度模型按照公式(5)进行交叉.
(4) 选择操作.将交叉之后的模型进行正演计算，按照(6)式判断是否接受.在两个种群均完成上述操作后，选出这两种群中最优的模型best12，进行步骤(5).
(5) 混沌搜索操作.利用公式(10)对电阻率模型进行混沌搜索，利用公式(11)对密度模型进行宽范围混沌搜索：
|
(10) |
|
(11) |
对最优解个体Xbest12的邻域内进行混沌搜索，若有比最优解更好的解则更新最优解，否则不保存进行下一步.
(6) 重复步骤(2)到(5)，直到满足最大迭代次数.
3 模型试验为了验证改进DE算法的效果，以及基于改进DE算法的MT和重力宽范围物性约束联合反演的有效性，下面设计了三个模型试验.模型试验1和模型试验2分别模拟了先验物性信息“精确”和“不精确”的条件，模型试验3则进一步讨论了物性分布部分共界面条件下，宽范围物性约束联合反演的适用性.
3.1 模型试验1图 1为设计的电阻率和密度模型，从左至右，三个异常体的电阻率分别为50、300和1500 Ωm，密度分别为2.58、2.65和2.70 g·cm-3，异常体的上下边界埋深分别为2 km和4 km.MT模拟频率范围为0.001~320 Hz，按对数间隔取38个频点，测点间距1000 m.
|
图 1 理论模型(色标与图 3、5、7一致) (a)电阻率理论模型；(b)密度理论模型. Fig. 1 Theoretical model (color scale is the same as the Figs. 3, 5 and 7) (a) Synthetic resistivity model; (b) Synthetic density model. |
为了验证改进的DE算法和联合反演的有效性，并验证宽范围物性约束技术在先验物性条件“精确”条件下的适用性，本文设计了MT和重力单独反演以及三种联合反演方案试验，四种反演方案具体如下：方案一，基于改进DE算法的MT和重力单独反演；方案二，基于未改进DE算法的传统点对点映射式岩石物性约束联合反演；方案三：基于改进DE算法的传统点对点映射式岩石物性约束联合反演；方案四：基于改进DE算法的宽范围岩石物性约束联合反演.四种方案DE算法最大种群数量均设置为300，最大迭代次数为50次.在联合反演过程中，将异常体的真实埋深浮动±15%，真实电阻率扰动±30%，真实密度扰动±0.02 g·cm-3作为反演参数的解空间.在基于非线性优化算法的联合反演中，解空间一般是根据先验信息设定的，而初始模型则是在解空间内随机生成的(杨辉等，2002；于鹏等，2008；胡祖志等，2020).基于此，本文的初始模型是物性参数在其解空间内随机生成的.模型试验中所模拟的电阻率和密度之间的先验物性关联见图 2，可以看出此时的先验物性关联是“精确”的.
|
图 2 电阻率和密度的岩石物性关联图 Fig. 2 Physical correlation diagram of the resistivity and density |
图 3为四种反演方案的结果.对比方案一单独反演和方案二、三、四联合反演的效果，四种反演方案的重力异常曲线虽然都已拟合，但是单独反演结果(图 3d)不仅界面刻画较差，密度远远偏离真实模型，从图 4的不同方案的反演结果与真实模型的残差分析也可以看出，三种联合反演方案重力效果均要优于单独反演.可见，联合反演可以改善单独重力反演的结果，在界面和物性还原上有更好的效果.
|
图 3 不同方案的结果 (a)电阻率初始扰动模型；(b)密度初始扰动模型；(c)方案一电阻率单独反演结果；(d)方案一密度单独反演结果；(e)方案二电阻率联合反演结果；(f)方案二密度联合反演结果；(g)方案三电阻率联合反演结果；(h)方案三密度联合反演结果；(i)方案四电阻率联合反演结果；(j)方案四密度联合反演结果；(k)均方误差曲线；(l)重力异常拟合曲线. Fig. 3 Results of different schemes (a) Initial perturbed resistivity mode; (b) Initial perturbed density model; (c) Resistivity separate inversion result of the scheme 1; (d) Density separate inversion result of the scheme 1; (e) Resistivity joint inversion result of the scheme 2; (f) Density joint inversion result of the scheme 2; (g) Resistivity joint inversion result of the scheme 3; (h) Density joint inversion result of the scheme 3; (i) Resistivity joint inversion result of the scheme 4; (j) Density joint inversion result of the scheme 4; (k) Mean square error curves; (l) Gravity anomaly fitting curves. |
|
图 4 不同方案的反演结果与理论模型的残差图 (a)方案一电阻率单独反演结果；(b)方案一密度单独反演结果；(c)方案二电阻率联合反演结果；(d)方案二密度联合反演结果；(e)方案三电阻率联合反演结果；(f)方案三密度联合反演结果；(g)方案四电阻率联合反演结果；(h)方案四密度联合反演结果. Fig. 4 Residual plots of different schemes of inversion results and synthetic model (a) Resistivity separate inversion result of the scheme 1; (b) Density separate inversion result of the scheme 1; (c) Resistivity joint inversion result of the scheme 2; (d) Density joint inversion result of the scheme 2; (e) Resistivity joint inversion result of the scheme 3; (f) Density joint inversion result of the scheme 3; (g) Resistivity joint inversion result of the scheme 4; (h) Density joint inversion result of the scheme 4. |
对比三种联合反演方案试验，从界面的还原情况上分析，方案二界面起伏比较严重，方案三和方案四界面还原程度则更高，比较好地刻画出了理论模型的上下边界；从物性的还原情况上分析，方案三和方案四无论是电阻率还是密度的反演结果比方案二都更加连续；从收敛速度上分析(图 3k)可以看出，方案三和方案四的收敛速度都要快于方案二的收敛速度，方案三在15次迭代左右的拟合精度就与方案二在50次迭代左右的拟合精度相当，方案四在20次迭代左右也达到了相同精度；从重力异常的拟合情况上分析，方案三和方案四反演结果对应的重力异常与理论的重力异常更贴合.
3.2 模型试验2模型试验1验证了改进的DE算法的效果，基于3.1节中的地球物理模型，本文还设计了先验物性信息存在一定偏差的地球物理模型(图 5)，从左至右，三个异常体的电阻率分别为50、500和1500 Ωm，密度分别为2.60、2.65和2.72 g·cm-3，异常体的上下边界埋深、MT频点个数、测点距与模型试验1设置相同.模型试验中所模拟的电阻率和密度之间的先验物性关联见图 6，可以看出此时的先验物性关联是“不精确”的.
|
图 5 理论模型(色标与图 7一致) (a)电阻率理论模型；(b)密度理论模型. Fig. 5 Theoretical model (color scale is the same as the Fig. 7) (a) Synthetic resistivity model; (b) Synthetic density model. |
|
图 6 电阻率和密度的岩石物性关联图 Fig. 6 Physical relationship between the resistivity and density |
为了对比宽范围物性约束技术与传统的点对点物性约束方式的特点，并验证宽范围物性约束技术在先验物性条件“不精确”条件下的适用性，本文设计了如下的对比试验：方案一，基于改进DE算法传统的点对点映射式岩石物性约束联合反演；方案二，基于未改进DE算法宽范围岩石物性联合反演；方案三，基于改进DE算法的宽范围岩石物性约束方案.三种方案DE最大种群数量设置为300，最大迭代次数为50次.在联合反演过程中，将异常体的真实埋深浮动±15%，真实电阻率扰动±30%，真实密度扰动±0.02 g·cm-3作为反演参数的解空间.
不同方案的联合反演结果见图 7.对比方案二(图 7e—f)和方案三(图 7g—h)反演结果：方案三即基于改进DE算法较之于方案二未改进DE算法的宽范围约束方案，界面起伏情况不仅有较好改善，而且无论是电阻率还是密度反演结果都更加的连续，对比方案三(图 9e—f)与方案二(图 9c—d)的残差图也可以明显看出，方案三的反演结果更加逼近理论模型.再从图 7i均方误差曲线分析可以看出，方案三从第5次迭代开始，均方误差下降效果就已经开始优于方案二，最后能够收敛到比方案二更高的精度，图 7j中的重力拟合曲线也更加贴合.因此，改进的DE算法方案效果优于未改进DE算法方案；对比方案一(图 7c—d)和方案三(图 7g—h)反演结果:方案一即改进DE算法传统岩石物性约束方案，结合图 6岩石物性关系存在的偏差情况一起分析，可以看出由于“不精确”先验物性信息的引入，第一个和第三个异常体的密度反演结果偏低，第二个异常体的密度反演结果偏高，且图 9b也可以明显看出和理论模型存在巨大差异.而方案三即改进DE算法宽范围岩石物性约束方案，电阻率和密度反演结果都比方案一的反演结果更加连续，第二个异常体密度也能在背景中清楚分辨出来，三个异常体物性反演结果也都更加逼近理论模型.整体分析以上三种方案，可以看出，方案二和方案三采用宽范围物性技术，所获得反演结果均可在一定程度上摆脱“不精确”先验物性信息引入的影响.
|
图 7 不同方案的结果 (a)电阻率初始扰动模型；(b)密度初始扰动模型；(c)方案一电阻率联合反演结果；(d)方案一密度联合反演结果；(e)方案二电阻率联合反演结果；(f)方案二密度联合反演结果；(g)方案三电阻率联合反演结果；(h)方案三密度联合反演结果；(i)均方误差曲线；(j)重力异常拟合曲线. Fig. 7 Results of different schemes (a) Initial perturbed resistivity mode; (b) Initial perturbed density model; (c) Resistivity joint inversion result of the scheme 1; (d) Density joint inversion result of the scheme 1; (e) Resistivity joint inversion result of the scheme 2; (f) Density joint inversion result of the scheme 2; (g) Resistivity joint inversion result of the scheme 3; (h) Density joint inversion result of the scheme 3; (i) Mean square error curves; (j) Gravity anomaly fitting curves. |
|
图 8 不同方案随迭代次数的物性耦合 Fig. 8 Physical coupling diagram of different schemes with iteration number |
|
图 9 不同方案的反演结果与理论模型的残差图 (a)方案一电阻率联合反演结果；(b)方案一密度联合反演结果；(c)方案二电阻率联合反演结果；(d)方案二密度联合反演结果；(e)方案三电阻率联合反演结果；(f)方案三密度联合反演结果. Fig. 9 Residual plots of different schemes of inversion results and synthetic model (a) Resistivity joint inversion result of the scheme 1; (b) Density joint inversion result of the scheme 1; (c) Resistivity joint inversion result of the scheme 2; (d) Density joint inversion result of the scheme 2; (e) Resistivity joint inversion result of the scheme 3; (f) Density joint inversion result of the scheme 3. |
图 8是三种方案随迭代次数变化的物性耦合图.本文选取了初始种群中的最优个体所对应的电阻率和密度制作第0次迭代对应的物性耦合图.由于三种方案采用了同样的初始模型，所以三种方案的物性耦合图是一致的，电阻率和密度随机分布解空间中.随着迭代的进行，由于“不精确”先验信息的影响，采用传统的点对点映射式物性关联的方案一对应的电阻率和密度始终分布在先验的物性关联信息上，而采用宽范围物性约束的方案二和方案三，在发挥先验物性信息导向作用的同时，还能发挥DE算法的寻优能力，逐渐搜索到真实值.此外，方案二基于未改进DE算法的联合反演结果对应的物性耦合分布比较分散，结合图 7的反演结果，可以看出，方案三的寻优能力是优于方案二的.
3.3 模型试验3为了进一步模拟更为复杂的情况，本文设计了电阻率和密度部分共界面的地球物理模型(图 11a—b)，在均匀物性背景下，电阻率模型存在三个异常体，密度模型存在两个异常体.电阻率模型异常块体A上下界面分别为2 km和4 km，电阻率为60 Ωm.块体B埋深与异常体A相同，电阻率为50 Ωm.块体C上下界面分别为0.5 km和2 km，电阻率为10 Ωm.密度模型只有两个异常块体A、B，密度分别为2.65、2.69 g·cm-3.本次模型试验中，MT模拟频率范围、测点间距、种群数量、迭代次数与模型试验2相同，上下界面扰动范围为真实值的±15%.
假设目前的先验信息不足以判定右边的异常块体是否共界面分布，先验信息仅有：左侧异常体电阻率分布的大致范围30~90 Ωm，密度分布的大致范围2.63~2.67 g·cm-3；右侧异常体电阻率分布的大致范围为5~60 Ωm，密度分布的大致范围为2.67~2.71 g·cm-3，同时也将这些物性分布范围设为本次模型试验的解空间.
本文首先根据先验物性信息以及左右异常体的空间位置，建立了两套岩石物性关联.具体而言，根据左边异常块体物性参数分布范围的上限lnln90 Ωm，2.67 g·cm-3，以及下限lnln30 Ωm，2.63 g·cm-3，利用两点公式，建立物性线性关系，密度和电阻率满足：
|
(12) |
对于右边异常块体，同样利用物性参数分布范围的上限lnln60 Ωm，2.71 g·cm-3，下限lnln5 Ωm，2.67 g·cm-3建立物性关联，密度和电阻率满足：
|
(13) |
图 10为基于上述方式建立的电阻率和密度的岩石物性关联图.
|
图 10 电阻率和密度的岩石物性关联图 Fig. 10 Physical relationship between the resistivity and density |
然后，按照均匀物性背景条件下存在三个电阻率和密度完全共界面分布的异常体条件进行初始建模，即三个异常体在解空间内随机生成初始模型(图 11c—d).
|
图 11 联合反演结果 (a)电阻率理论模型；(b)密度理论模型；(c)电阻率初始扰动模型；(d)密度初始扰动模型；(e)电阻率联合反演结果；(f)密度联合反演结果；(g)电阻率与理论模型残差图；(h)密度与理论模型残差图；(i)均方误差曲线；(j)重力异常拟合曲线. Fig. 11 Joint inversion result (a) Synthetic resistivity model; (b) Synthetic density model; (c) Initial perturbed resistivity mode; (d) Initial perturbed density model; (e) Resistivity joint inversion result; (f) Density joint inversion result; (g) Residual plots of resistivity result and synthetic model; (h) Residual plots of density result and synthetic mode; (i) Mean square error curves; (j) Gravity anomaly fitting curves. |
最后，利用先验物性关联信息，即(12)和(13)式进入电阻率宽范围岩石物性约束流程，电阻率、密度联合反演结果见图 11e—f，电阻率、密度联合反演结果与理论模型的残差对比情况见图 11g—h，联合反演迭代误差曲线见图 11i，重力异常拟合情况见图 11j，联合反演结果物性耦合情况见图 12.
|
图 12 联合反演物性耦合图 Fig. 12 Coupling physical properties in joint inversion results |
分析图 11e—j结果可以看出，异常体的界面和物性均得到了较好的还原，重力异常拟合曲线也比较贴合.与异常块体A以及B的密度反演值与背景密度具有较明显的差异相比，异常块体C的密度反演值与背景密度接近，提示异常块体C与均匀背景可能仅存在电性差异，而密度差异微弱.此外，进一步分析物性耦合图(图 12)可以看出，电阻率和密度联合反演结果在先验物性信息周边一定范围内实现了耦合.并且，异常块体C的反演密度值在2.68 g·cm-3附近，与背景密度十分接近.综上，本文推测异常块体C处仅为电性异常体并非为密度异常体.
通过该试验可以看出，宽范围约束技术适用于由于岩石物性部分共界面导致的物性关联特征的“不精确”的条件，在一定“范围内”实现耦合也为还原物性部分共界面的真实模型提供了“空间”.
4 结论本文以双种群设置方案为框架改进了DE算法，并实现了基于DE算法的MT和重力宽范围物性约束联合反演.通过模型试验，可以看出：
(1) 与传统DE算法相比，本文改进的DE算法的寻优效果更强，收敛速度更快，并且成功运用到MT和重力联合反演中.
(2) 本文提出的基于改进DE算法的宽范围物性约束技术具有一定的容错性，适用于先验信息“不精确”的条件，可以降低联合反演对先验信息的要求，在一定程度上适用于岩石物性部分共界面的条件.
(3) 本文算法的实现，进一步验证了宽范围物性约束思想在联合反演领域的适用性，可以预见宽范围物性约束思想具有融入其他非线性优化算法的寻优过程，进而在一定“范围”内实现岩石物性参数耦合的潜质.
需要指出的是，目前宽范围物性约束思想已结合模拟退火算法和DE算法得以实现，如何将宽范围物性约束思想融入线性优化算法是作者今后的研究重点.
Bao L, Zeng J C. 2011. A hi-group differential artificial bee colony algorithm. Control Theory and Applications (in Chinese), 28(2): 266-272. |
Carter-McAuslan A, Lelièvre P G, Farquharson C G. 2015. A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios. Geophysics, 80(1): W1-W15. DOI:10.1190/geo2014-0056.1 |
Chen J, Wen N, Chen B Y. 2007. Joint inversion of gravity-magnetic-electrical-seismic combination survey:progress and prospect. Progress in Geophysics (in Chinese), 22(5): 1427-1438. DOI:10.3969/j.issn.1004-2903.2007.05.013 |
Chen X, Yu P, Deng J Z, et al. 2016. Joint inversion of MT and seismic data based on wide-range petrophysical constraints. Chinese Journal of Geophysics (in Chinese), 59(12): 4690-4700. DOI:10.6038/cjg20161228 |
Chen X, Yu P, Deng J Z, et al. 2017. A new framework for geophysical joint inversion. Oil Geophysical Prospecting (in Chinese), 52(4): 851-858, 883. |
Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data:Application to prestack depth imaging. The Leading Edge, 26(3): 326-331. |
Colombo D, Rovetta D. 2018. Coupling strategies in multiparameter geophysical joint inversion. Geophysical Journal International, 215(2): 1171-1184. DOI:10.1093/gji/ggy341 |
Dell'Aversana P. 2001. Integration of seismic, MT and gravity data in a thrust belt interpretation. First Break, 19(6): 335-341. DOI:10.1046/j.1365-2397.2001.00158.x |
Dell'Aversana P, Bernasconi G, Chiappa F. 2016. A global integration platform for optimizing cooperative modeling and simultaneous joint inversion of multi-domain geophysical data. AIMS Geosciences, 2(1): 1-31. |
Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data. Geophysical Research Letters, 30(13): 1658. |
Gao Z Q, Pan Z B, Zuo C, et al. 2019. An optimized deep network representation of multimutation differential evolution and its application in seismic inversion. IEEE Transactions on Geoscience and Remote Sensing, 57(7): 4720-4734. DOI:10.1109/TGRS.2019.2892567 |
Guo M. 2018. A Bayesian joint inversion of magnetotelluric and gravity based on the wide-range of petrophysical constraints[Master's thesis] (in Chinese). Nanchang: East China University of Technology.
|
Haber E, Oldenburg D. 1997. Joint inversion:a structural approach. Inverse Problems, 13(1): 63. |
Heincke B, Jegen M, Hobbs R, et al. 2006. Joint Inversion of MT, Gravity and Seismic Data Applied to Sub-basalt Imaging. SEG: 784-789. DOI:10.1190/1.2370374 |
Hu Z Z, Shi Y L, Liu Y X, et al. 2020. Nonlinear constrained joint inversion of MT and gravity data. Oil Geophysical Prospecting (in Chinese), 55(1): 226-232. |
Jegen M D, Hobbs R W, Tarits P, et al. 2009. Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints:Preliminary results of sub-basalt imaging off the Faroe Shelf. Earth and Planetary Science Letters, 282(1-4): 47-55. DOI:10.1016/j.epsl.2009.02.018 |
Liu G, Meng X, Guo L. 2011. The gravity and seismic sequential inversion and its GPU implementation. International Workshop on Gravity, Electrical & Magnetic Methods and Their Applications, Beijing, China: 10-13. DOI:10.1190/1.3659038 |
Lu Y L, Zhou J Z, Li Y H, et al. 2008. Adaptive differential evolution algorithm combined with chaotic search. Computer Engineering and Applications (in Chinese), 44(10): 31-33, 39. |
Maier R, Heinson G, Tingay M, et al. 2010. Improved imaging of the subsurface using a gravity and magnetotelluric joint inversion.//21st Geophysical Conference. Sydney, 1-4.
|
McGinley B, Maher J, O'Riordan C, et al. 2011. Maintaining healthy population diversity using adaptive crossover, mutation, and selection. IEEE Transactions on Evolutionary Computation, 15(5): 692-714. DOI:10.1109/TEVC.2010.2046173 |
Meng H Y, Zhang X H, Liu S Y. 2008. A differential evolution based on double populations for constrained multi-objective optimization problem. Chinese Journal of Computers (in Chinese), 31(2): 228-235. |
Moorkamp M, Heincke B, Jegen M, et al. 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data. Geophysical Journal International, 184(1): 477-493. DOI:10.1111/j.1365-246X.2010.04856.x |
Peng M, Tan H D, Jiang M, et al. 2013. Three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints. Chinese Journal of Geophysics (in Chinese), 56(8): 2728-2738. DOI:10.6038/cjg20130821 |
Song W Q, Gao Y K, Zhu H W. 2013. The differential evolution inversion method based on Bayesian theory for micro-seismic data. Chinese Journal of Geophysics (in Chinese), 56(4): 1331-1339. DOI:10.6038/cjg20130427 |
Storn R, Price K. 1997. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4): 341-359. DOI:10.1023/A:1008202821328 |
Sun J J, Li Y G. 2011. Geophysical inversion using petrophysical constraints with application to lithology differentiation.//12th International Congress of the Brazilian Geophysical Society. Rio de Janeiro, Brazil: SEG, 2644-2648.
|
Tan Y, Tan G Z, Wu X D. 2010. Dual population differential evolution algorithm based on crossover and mutation strategy. Computer Engineering and Applications (in Chinese), 46(18): 9-12. |
Wang C, Gao J H. 2012. High-dimensional waveform inversion with cooperative coevolutionary differential evolution algorithm. IEEE Geoscience and Remote Sensing Letters, 9(2): 297-301. DOI:10.1109/LGRS.2011.2166532 |
Wang T Y. 2015. Comprehensive research on the iterative finite element and the improved differential evolution of magnetotelluric[Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences.
|
Wang Z J, Zhan Z H, Lin Y, et al. 2018. Dual-strategy differential evolution with affinity propagation clustering for multimodal optimization problems. IEEE Transactions on Evolutionary Computation, 22(6): 894-908. DOI:10.1109/TEVC.2017.2769108 |
Wu L H, Wang Y A, Zhou S W, et al. 2007. Research and application of pseudo parallel differential evolution algorithm with dual subpopulations. Control Theory and Application (in Chinese), 24(3): 453-458. |
Yang H, Dai S K, Song H B, et al. 2002. Overview of joint inversion of integrated geophysics. Progress in Geophysics (in Chinese), 17(2): 262-271. DOI:10.3969/j.issn.1004-2903.2002.02.011 |
Yin C C, Sun S Y, Gao X H, et al. 2018. 3D joint inversion of magnetotelluric and gravity data based on local correlation constraints. Chinese Journal of Geophysics (in Chinese), 61(1): 358-367. DOI:10.6038/cjg2018K0765 |
Yu P, Wang J L, Wu J S, et al. 2006. Review and discussions on geophysical joint inversion. Progress in Exploration Geophysics (in Chinese), 29(2): 87-93, 134. |
Yu P, Dai M G, Wang J L, et al. 2008. Joint inversion of gravity and seismic data based on common gridded model with random density and velocity distributions. Chinese Journal of Geophysics (in Chinese), 51(3): 845-852. DOI:10.3321/j.issn:0001-5733.2008.03.025 |
Zhang L. 2016. Regularization joint inversion of two-dimensional magnetotelluric and gravity based on petrophysical constraints[Master's thesis] (in Chinese). Nanchang: East China University of Technology.
|
暴励, 曾建潮. 2011. 一种双种群差分蜂群算法. 控制理论与应用, 28(2): 266-272. |
陈洁, 温宁, 陈邦彦. 2007. 重磁电震联合反演研究进展与展望. 地球物理学进展, 22(5): 1427-1438. DOI:10.3969/j.issn.1004-2903.2007.05.013 |
陈晓, 于鹏, 邓居智, 等. 2016. 基于宽范围岩石物性约束的大地电磁和地震联合反演. 地球物理学报, 59(12): 4690-4700. DOI:10.6038/cjg20161228 |
陈晓, 于鹏, 邓居智, 等. 2017. 地球物理联合反演新框架研究. 石油地球物理勘探, 52(4): 851-858, 883. |
郭曼. 2018.基于宽范围岩石物性约束的MT与重力贝叶斯联合反演[硕士论文].南昌: 东华理工大学.
|
胡祖志, 石艳玲, 刘云祥, 等. 2020. 电磁与重力数据非线性约束联合反演. 石油地球物理勘探, 55(1): 226-232. |
卢有麟, 周建中, 李英海, 等. 2008. 基于混沌搜索的自适应差分进化算法. 计算机工程与应用, 44(10): 31-33, 39. |
孟红云, 张小华, 刘三阳. 2008. 用于约束多目标优化问题的双群体差分进化算法. 计算机学报, 31(2): 228-235. |
彭淼, 谭捍东, 姜枚, 等. 2013. 基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演. 地球物理学报, 56(8): 2728-2738. DOI:10.6038/cjg20130821 |
宋维琪, 高艳珂, 朱海伟. 2013. 微地震资料贝叶斯理论差分进化反演方法. 地球物理学报, 56(4): 1331-1339. DOI:10.6038/cjg20130427 |
谭跃, 谭冠政, 伍雪冬. 2010. 基于交叉变异策略的双种群差分进化算法. 计算机工程与应用, 46(18): 9-12. |
王天意. 2015.大地电磁迭代有限元与改进差分进化正反演算法研究[博士论文].北京: 中国地质大学.
|
吴亮红, 王耀南, 周少武, 等. 2007. 双群体伪并行差分进化算法研究及应用. 控制理论与应用, 24(3): 453-458. |
杨辉, 戴世坤, 宋海斌, 等. 2002. 综合地球物理联合反演综述. 地球物理学进展, 17(2): 262-271. DOI:10.3969/j.issn.1004-2903.2002.02.011 |
殷长春, 孙思源, 高秀鹤, 等. 2018. 基于局部相关性约束的三维大地电磁数据和重力数据的联合反演. 地球物理学报, 61(1): 358-367. DOI:10.6038/cjg2018K0765 |
于鹏, 王家林, 吴健生, 等. 2006. 地球物理联合反演的研究现状和分析. 勘探地球物理进展, 29(2): 87-93, 134. |
于鹏, 戴明刚, 王家林, 等. 2008. 密度和速度随机分布共网格模型的重力与地震联合反演. 地球物理学报, 51(3): 845-852. DOI:10.3321/j.issn:0001-5733.2008.03.025 |
张磊. 2016.基于岩石物性约束的MT与重力二维正则化联合反演[硕士论文].南昌: 东华理工大学.
|