基于震电效应的矿化度的反演方法

1.本发明属于地质探测领域，具体涉及一种基于震电效应的矿化度的反演方法。


背景技术：

2.在土壤、岩石等天然孔隙介质中，孔隙流体(通常为水)中通常会含有带电离子，形成孔隙电解质溶液。孔隙地层总体上是电中性的，但是由于固相骨架选择性吸附孔隙溶液中的某种离子，使得原本中性的溶液中存在净剩电荷，在孔壁附近形成了双电层。由于置换和破键作用，使得固相基质表面带负电荷，溶液中的部分阳离子就被吸附到矿物的表面，形成吸附层，同时，溶液中的离子重新分布，形成扩散层，吸附层和扩散层一起构成了双电层。当弹性波在孔隙介质中传播时，孔隙流体相对固相骨架运动，孔隙电解质溶液中的净剩带电离子随之运动形成电流和电磁场。因此，在孔隙介质中，不仅固相运动与流体运动是耦合的，而且弹性波场与电磁场也是耦合的。frenkel最先从理论上研究了渗流运动与电磁场的耦合，在他的工作基础上，biot建立一套完整的流体饱和孔隙介质声学理论，pride在frenkel和biot的工作基础上，建立宏观震电耦合的控制方程组。目前震电效应相关的研究主要是集中在野外和实验室观测，以及部分学者对其进行数值模拟研究，包括动电效应，地球发电机效应，压电效应，压磁效应等不同机制。(pride&haartsen,1996；haartsen&pride,1997；hu&wang,1999；garambois&dietrich,2002；guan&hu,2008；gao&hu,2010；hu&gao,2011).
3.油气资源的储层主要为含流体多孔介质。单一的地震勘探无法区分地下储层流体性质，单一的电法勘探不能准确区分地下储层位置。震电勘探集中了地震勘探和电法勘探的优点，可以较为准确的分析出相应地层的未知地层参数。因此在地球物理勘探中，针对含流体多孔介质，基于双电层动电效应的震电勘探方法由于对多孔介质的参数敏感，使它们成为探测地下介质的潜在手段，将在油气和水资源勘探领域具有广阔的应用前景。
4.但是目前，很少有学者对利用震电信号进行反演研究，guan et al.(2013)提出了一种在流体饱和的多孔地层中反演震电测井渗透率的方法。macchioli-grande et al.(2020)研究了利用贝叶斯反演从一组含噪声的合成数据推断冰川介质参数的可能性。然而，这些反演研究都是基于均匀半空间的假设，对于地下结构进行成像还是存在一定缺陷。


技术实现要素：

5.鉴于上述，本发明的目的是提供一种基于震电效应的矿化度的反演方法，能够实现利用震电数据进行地下矿化度参数的探测。
6.为实现上述发明目的，本发明提供以下技术方案：
7.一种基于震电效应的矿化度的反演方法，包括以下步骤：
8.(1)根据反演地质构建初始模型，其中，模型参数为矿化度；
9.(2)获取能够用于反演的观测地震电磁场数据；
10.(3)采用pride方程组对模型进行正演计算，以得到作为正演结果的预测地震电磁
场数据；当第一次迭代计算时，模型为初始模型，其他次迭代计算时，模型为更新模型；
11.(4)根据观测地震电磁场数据、预测地震电磁场数据以及模型粗糙度构建震电效应一维正则化反演的目标函数；
12.(5)以目标函数求极小等价于其导数为0为目标，采用高斯-牛顿方法进行迭代反演计算模型参数的更新量，以得到更新模型；
13.(6)迭代重复步骤(3)-(5)，直到迭代终止，得到的最新更新模型作为最终的反演模型。
14.在一个实施例中，在对目标函数进行求解时，计算观测地震电磁场数据的灵敏度，具体计算为：
[0015][0016]
其中，mj表示第j个模型的模型参数，δmj表示模型参数的更新量，fs(mj)表示对第j个模型的正演响应数据，即预测地震电磁场数据，s表示正演响应数据的索引，j
sj
表示第j个模型关于第s个正演响应数据构建的目标函数；
[0017]
依据观测地震电磁场数据的灵敏度筛选观测地震电磁场数据用于模型参数的更新量的计算。
[0018]
在一个实施例中，所述依据观测地震电磁场数据的灵敏度筛选观测地震电磁场数据用于模型参数的更新量的计算，包括：将超过灵敏度阈值的观测地震电磁场数据对应的噪声误差设置为极大，以使这些数据在反演过程中不发挥作用。
[0019]
在一个实施例中，在迭代反演计算过程中，为反演过程中施加电导率上下限约束以减少虚假异常数据，具体包括：
[0020]
基于对数参数的电导率转换函数将模型参数变化到转换域内，对反演结果进行电导率上下限约束，以减少了反演的多解性，即利用下述公式将模型参数mk变换到转换域内模型参数hk，迭代之后再将hk利用下述公式变换为mk，其中ak为下限，取值范围为小于0.001，bk为上限，取值范围为大于等于1，
[0021]hk
＝log(m
k-ak)-log(b
k-mk),ak＜mk＜bk,k＝1,2,...,m.
[0022][0023]
在一个实施例中，在迭代反演计算过程中，通过对目标函数进行求解以获得预测地震电磁场数据和观测地震电磁场数据的拟合差，将拟合差小于设置阈值作为反演迭代终止条件，在不满足反演迭代终止条件时，依据求解目标函数获得的模型参数的更新量来更新模型。
[0024]
与现有技术相比，本发明具有的有益效果至少包括：
[0025]
本发明实施例提供的基于震电效应的矿化度的反演方法中，在为探测地质结构构建初始模型和获取观测地震电磁场数据的基础上，依据观测地震电磁场数据、正演得到的预测地震电磁场数据以及模型粗糙度来构建震电效应一维正则化反演的目标函数，采用高斯-牛顿方法进行迭代反演计算模型参数的更新量，以得到更新模型，这样可以获取符合探测地质结构的最终反演模型，大大提升了震电反演的效率和准确度。
附图说明
[0026]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。
[0027]
图1是本发明实施例提供的基于震电效应的矿化度的反演方法的流程图；
[0028]
图2是本发明实施例提供的反演模型结果与真实模型对比图，中间异常层的位置分别位于(a)2-4km；(b)4-6km；(c)6-8km；(d)11-13km.蓝色虚线是震源所在深度；
[0029]
图3是本发明实施例提供的观测地震电磁场数据与预测地震电磁场数据以及残差的对比图，中间异常层的位置分别位于(a)2-4km；(b)4-6km；(c)6-8km；(d)11-13km，残差放大了10倍绘制；
[0030]
图4是本发明实施例提供的反演模型结果的协方差矩阵，中间异常层的位置分别位于(a)2-4km；(b)4-6km；(c)6-8km；(d)11-13km。
具体实施方式
[0031]
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
[0032]
图1是本发明实施例提供的基于震电效应的矿化度的反演方法的流程图。如图1所示，实施例提供的基于震电效应的矿化度的反演方法，包括以下步骤：
[0033]
步骤1，为探测地质结构构建初始模型。
[0034]
在构建初始模型时，设置几何模型信息、传递数据信息以及控制参数。其中，几何模型信息包括弹性性质、电性性质等，具体参数见表1，传递数据信息包括震源及接收点个数与位置等。控制参数包括最大迭代次数和正则化因子λ等。在迭代过程中需要更新的模型参数即为矿化度。
[0035]
表1
[0036]
[0037][0038]
步骤2，获取能够用于反演的观测地震电磁场数据。
[0039]
实施例中，用于反演的观测地震电磁场数据可以来自有实测数据，当为实测数据时，在被应用之前，还需要进行去噪预处理。
[0040]
步骤3，采用pride方程组对模型进行正演计算，以得到作为正演结果的预测地震电磁场数据。
[0041]
实施例中，采用pride方程组对模型进行正演计算时，pride方程组表示为：
[0042][0043][0044][0045][0046][0047][0048]
其中，表示拉普拉斯算符，σ(ω)表示动态电导率，ω表示圆频率，u是固相位移，w＝φ(u-uf)是渗流位移，uf是平均流体位移，φ是孔隙度，p是孔隙流体压强，τ是应力张量，i是单位张量，ρf是孔隙流体密度，ρ是孔隙介质的等效密度，η是孔隙流体粘滞系数，κ(ω)是动态渗透率，l(ω)表示动电耦合系数，h,c,m和g是四个单独的弹性模量，h是磁场，e是电场，f和f是分别作用在流相和整个孔隙介质上的体力密度，ε和μ是孔隙介质的介电常数和磁导率；
[0049]
将模型参数和中心频率带入pride方程组，经计算求解得到正演结果；其中，在计算固相位移u、渗流位移w以及平均流体位移uf时需要用到中心频率。其中，模型参数为表1中所示的所有参数。
[0050]
由于步骤3是迭代计算过程，在第一次迭代时，模型为初始模型，除第一次迭代外的其他次迭代计算时，模型为更新模型。
[0051]
步骤4，根据观测地震电磁场数据、预测地震电磁场数据以及模型粗糙度构建震电
效应一维正则化反演的目标函数。
[0052]
实施例中，基于l2范数的目标函数，该目标函数用于正则化反演过程来优化反演更新模型，具体根据观测地震电磁场数据、正演得到的预测地震电磁场数据以及模型粗糙度构建正则化反演的目标函数为：
[0053][0054]
其中为数据拟合项，为模型粗糙度约束项，β为正则化因子，u和k的具体形式为：
[0055]
u＝wd(d
obs-d
prd
)，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0056]
k＝wm(m-m
ref
).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0057]
上式中，d
obs
为观测数据向量，d
prd
是由模型m计算的电磁响应，m
ref
为参考模型，如果初始给定参考模型，即参考模型不变；如果未给定，则第一次迭代为初始模型，之后为每一次迭代为更新的模型，wd为一个对角阵，其元素是观测响应中噪声的倒数。和的一般形式可表示为：
[0058][0059]
其中，x表示上式中的u或k。第n次迭代的目标函数可表示为
[0060][0061]
而
[0062]
u＝wd(d
obs-d
n-1-jδm)，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0063]
k＝wm(m
n-1
δm-m
ref
)，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0064]
式中，d
n-1
是由上一次迭代得到的模型m
n-1
计算得到的响应向量，δm＝m
n-m
n-1
，j是灵敏度矩阵。
[0065]
步骤5，以目标函数求极小等价于其导数为0为目标，采用高斯-牛顿方法进行迭代反演计算模型参数的更新量，以得到更新模型。
[0066]
实施例中，在对正则化反演的目标函数进行求解时，采用一阶向前差分的方法近似计算灵敏度矩阵(雅可比矩阵)，即计算观测地震电磁场数据的灵敏度。具体计算为：
[0067][0068]
其中，mj表示第j个模型的模型参数，δmj表示模型参数的更新量，fs(mj)表示对第j个模型的正演响应数据，即预测地震电磁场数据，s表示正演响应数据的索引，j
sj
表示第j个模型关于第s个正演响应数据构建的目标函数。
[0069]
实施例中，灵敏度用作观测地震电磁场数据的筛选基准，依据观测地震电磁场数据的灵敏度筛选观测地震电磁场数据用于模型参数的更新量的计算。具体地，将超过灵敏度阈值的观测地震电磁场数据对应的噪声误差设置为极大，以使这些观测地震电磁场数据在反演过程中不发挥作用。
[0070]
实施例中，还需要计算更新模型粗糙度，即在得到当前迭代的模型参数的更新量m时，根据参考模型m0，按照公式(7)和(9)计算模型粗糙度约束项φm(k)，以更新目标函数，用于下一次反演。
[0071]
实施例中，为了使反演结果更加接近真实情况，在迭代反演计算过程中，为反演过程中施加电导率上下限约束以减少虚假异常数据，具体包括：基于对数参数的电导率转换函数将模型参数变化到转换域内(称为对数变换)，对反演结果进行电导率上下限约束，以减少了反演的多解性。
[0072]
即利用公式(15)将模型参数mk变换到转换域内模型参数hk，迭代之后再将hk利用公式(16)变换为mk，其中ak为下限，取值范围为小于0.001，优选为1e-5，bk为上限，取值范围为大于等于1，优选为1，
[0073]hk
＝log(m
k-ak)-log(b
k-mk),ak＜mk＜bk,k＝1,2,...,m.
ꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0074][0075]
步骤6，迭代重复步骤3-步骤5，直到迭代终止，得到的最新更新模型作为最终的反演模型。
[0076]
实施例中，每次迭代计算时，根据上次迭代更新的更新模型继续执行步骤3-步骤5，直到到达迭代终止，满足反演终止条件或达到最大迭代次数，获取符合探测地质结构的最终反演模型，同时输出正演响应。
[0077]
具体地，在迭代反演计算过程中，通过对目标函数进行求解以获得预测地震电磁场数据和观测地震电磁场数据的拟合差，将拟合差小于设置阈值作为反演迭代终止条件，在不满足反演迭代终止条件时，依据求解目标函数获得的模型参数的更新量来更新模型。
[0078]
实施例中，反演迭代终止条件为拟合差小于1，其中，拟合差rms的计算为：
[0079][0080]
其中，ndata是观测地震电磁场数据个数，s为观测地震电磁场数据的索引，d
pre，s
和d
obs,s
是正演计算得到的正演响应(即预测地震电磁场数据)和观测地震电磁场数据，wd为数据方差矩阵。
[0081]
实施例还提供利用上述基于震电效应的矿化度的反演方法的反演算例，以验证反演算法的正确性。
[0082]
图2是本发明实施例提供的反演模型结果与真实模型对比图，中间异常层的位置分别位于(a)2-4km；(b)4-6km；(c)6-8km；(d)11-13km。横着虚线是震源所在深度；异常层厚度为2km，震源深度为10km，震源为双力偶点源m
xz
＝m
zx
＝m0＝4.37
×
10
17n·
m，相当于矩震级为5.7级的地震强度，接受点位置为(100km,5km,1m)，震源时间函数为：
[0083][0084]
其中，中心频率f0为0.5hz，脉冲持续时间为4s，频率范围0-1hz。
[0085]
图3是本发明实施例提供的观测地震电磁场数据与预测地震电磁场数据以及残差的对比图，中间异常层的位置分别位于(a)2-4km；(b)4-6km；(c)6-8km；(d)11-13km，残差放大了10倍绘制。从图2和3可以看出反演模型与真实模型吻合较好，且响应曲线也能够较好的拟合，虽然图2中的(d)存在一定差距，这是由于地震电磁场数据主要对震源深度以上的
地层参数灵敏度较大。这些都可以说明本发明的基于震电效应的矿化度反演算法是可靠的，具有一定的实用性。图4是本发明实施例提供的反演模型结果的协方差矩阵。从图4可以看出，反演是稳定收敛的且反演结果与真实模型的差距较小，进一步证明反演算法的正确性。
[0086]
以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。