CN112799140B - 一种基于自然电位反演的渗透率估计方法 - Google Patents

Abstract

本发明属于地热、地球物理电法勘探和工程与环境地球物理研究技术领域，是涉及一种基于自然电位信号反演的渗透率估计方法。包括：S1通过建立水力与电性耦合模型模拟地下自然电位信号的正反演响应获取地下电流密度分布，利用采集的电导率约束还原地下电势场，并根据电势场求取渗透率分布；通过利用现场获得的高密度电阻率数据进行倒转与处理获得电导率分布，结合自电位信号进行正反演计算恢复地下自然电势场，估计地下渗透率分布信息。相比于传统的侵入式方法，利用无损探测技术观测有效地热靶区的地下水多维(时间‑空间)多响应参数信息。将多个电性参数耦合计算恢复渗透率的分布，跳过传统问题中的采样难，测量难等步骤。

Description

一种基于自然电位反演的渗透率估计方法

技术领域

本发明属于地热、地球物理电法勘探和工程与环境地球物理研究技术领域，是涉及一种基于自然电位信号反演的渗透率估计方法。

背景技术

近些年，对地下热能源的探测与利用开发是持续发展的新型技术领域。通过常规的地球物理勘探对地下热流体的流动方向和迁移路径的信息获取是人们所关注的重点。

传统的地球物理方法对地下水流动信息不够敏感，无法检测到渗流的微小变化，相比于电阻率法是对水流变化进行观测的方法之一，可以相对准确的获得地电结构。其中，电法勘探中的自然电位法对水流的变化与路径信息的刻画是非常敏感的一种有效方法，对应于地下自然电流源的产生而对电场的被动测量探测技术，应用此方法可以很好的恢复有关地下水流量的非侵入式信息。

自然电位方法作为一种广泛应用的电磁场地球物理勘探方法，通过地下水的流动运移导致的电荷堆积，形成电离子固液双界面，进而产生微弱的电流密度源发散于地下空间被地表的采集仪器观测到。但在这类方法的研究早期，人们无法获得其观测信号与地下流体位置的相互响应关系，且方法论述只停留在定性解释方面。随着科学技术与理论计算的发展，在上个世纪自然电位法已经拥有了完整的理论体系与量化计算能力。它在工程环境中用于检测和监测大坝和堤防中的异常渗漏，地热区域的热变化与断裂引起水流流变等具有非常高的敏感性。通过观察地表的自然电位数据可以有效的计算出地下空间的电位分布与电流密度分布从而定量出异常分布特征。相比于其他传统的地下水水文观测方法，它的成本较低，具有非侵入式测量等优点，在地球物理和非侵入式流量传感器中拥有很好的前景。目前，自然电位法有待于深入探讨其对地下流体的多种响应机制，其中地层渗透率对流体的运移能力起到关键作用，这也是本专利的出发点，通过电性地球物理方法的无损探测技术还原与水文参数有关的丰富地下信息。

在热能源丰富的地下空间中，由于较高的温度引起的水流上涌或迁移运动是一种常见的异常变化情况。同时，在地热开发区域进行注水压裂，使得地下水流的流动会引起孔隙固液双界面形成吸附双电层，从而发生电离，形成电荷迁移而产生地下自然电位差，这种流体流动产生流动电势可以在地表被观察到。另外，这种热异常的温度差异也会产生热异常电动势。所以在地热区观测的自然电位信号通常是多种地下响应引起的。

可以结合多种响应的自然电位信号和靶区电导率分布信息进行反演来恢复地下电势分布，再通过已知的地热区域压裂信息与其他地层参考项加入到模型的运算中，可获得达西速度，渗透率等所关注的参数分布，对解释地热异常引起的流体流动和路径信息提供更为可靠的技术支持。

发明内容

本发明所要解决的技术问题在于提供一种基于自然电位反演的渗透率估计方法，相比于传统的侵入式方法，利用无损探测技术观测有效地热靶区的地下水多维(时间-空间)多响应参数信息。将多个电性参数耦合计算恢复渗透率的分布，跳过传统问题中的采样难，测量难等步骤，真正意义上实现“透过现象看本质”。

本发明是这样实现的，一种基于自然电位反演的渗透率估计方法，该方法包括：

S1通过建立水力与电性耦合模型模拟地下自然电位信号的正反演响应：

其中正演通过模拟地下水随压力变化的运移程度，再计算耦合电性的连接参数与描述水流流动的速度参数，最后通过电性控制方程计算地下的电流密度源与电导率约束条件还原电势分布；利用第二类边界条件在水力与电性耦合模型的顶部获取自然电位信号，模拟地下水流动激发的电性响应；正演过程的输入参数包括:电导率，渗透率与边界条件，输出参数为电流密度以及其表层接收的自然电位信号；

反演通过建立真实采集信号与正演模拟信号的差值的目标函数，利用最小二乘牛顿迭代或最速下降法使得目标解趋于最优值，其中，正演计算的结果不断迭代来拟合采集信号，来获得最优解状态下的模型响应；通过地表采集的自然电位信号与正演模拟的输出信号进行观测差最小化，获得地下电流密度分布；

S2根据步骤S1获得的地下电流密度分布，利用采集的电导率约束还原地下电势场，并根据电势场求取渗透率分布；

S3根据步骤S1和步骤S2利用多项电性参数与地下流体流动信息对水文观测区域或地热靶区进行评价，通过利用现场获得的高密度电阻率数据进行倒转与处理获得电导率分布，结合自电位信号进行正反演计算恢复地下自然电势场，估计地下渗透率分布信息。

进一步地，步骤S1中，正演包括利用公式(1)的本构方程和连续性来解释非均匀的和各向异性的多孔材料的运移程度：

其中，u表示达西速度，a为比存量，K为渗透率值，Q_s为外部源，h为液压头；

公式(1)遵循第一类边界条件Γ_D和第二类边界条件Γ_N，同时满足初始条件：

initial h＝h₀att＝0

在各向异性介质中，导电电流密度

与源电流密度j_s之和为总电流密度j：

j_s＝Q_Vu (3)

其中，σ表示为电导率值，ψ为电势，Q_V是每单位空隙体积的有效过剩电荷密度，通过经验公式(4)表达有效过剩电荷密度Q_V的渗透率依赖性：

log₁₀Q_V＝-9.2-0.82log₁₀K (4)

根据电荷的连续性方程

具有以下边界条件:

其中，σ表示为电导率值，ψ为电势。

进一步地，步骤S1中，

反演包括：用一组点M上的j_s

ψ(N)＝∫_ΩG(N，M)j_s(M)dV (7)

其中j_s(M)为一组点M上的电流密度源，ψ(N)为一组N电极处的电势， G(N，M)表示在点N处测量自然电位数据和当前位于点M处的内核矩阵，通过利用有限元法求解内核矩阵，

其中|| ||₂为L2范数，G表示N×2M的核矩阵，ω为与电势和内核矩阵G有关的参数项，即为电流密度值j_s，ψ^obs表示一组自电位数据相对应的N个元素的向量，W_d为N×N的对角线加权平方矩阵，且δ₁为偏差的平方数，加入正则化项与深度加权函数φ_m，建立全局目标函数项φ_T：

其中

为N×2M的深度加权矩阵，z₀表示模型单元的观测高度，L是基于ω的一阶导数的光滑度，β为常数项，介于1到4之间，λ表示正则化的约束项0＜λ＜∞，目标函数用标准形式给出：

其中

L_P是通过基于QR分解的Elédn算法推导出的矩阵，求目标函数φ_T最小值的通过公式(11)求取：

在地下渗透率参数已知的情况下计算出原始电流密度分布用于目标函数，在没有先验模型，即

的情况下得到公式(12)：

为减少核矩阵的计算量使用奇异值分解方法，将

H＝min(N，M)，并重新定义

其中，Λ为奇异对角矩阵，其元素Λ_i表示矩阵对角线上的奇异值分量，uv^T为左右奇异向量，

为反演计算后的更新电流密度参数，通过迭代误差的减少获得计算稳定值分布。

进一步地，步骤S2利用步骤S1中反演获得地下真实电流密度分布j_s，计算公式(14)的电导率σ约束还原地下电势场ψ结果：

通过达西速度方程

将公式(15)中的两种势场计算联立得到：

其中，压力p的空间变化等于液压头的变化

将其带入公式(16) 中获得：

认为压力梯度

的变化是与渗透无关的，在估计渗透率K时设为均匀半空间的理想压力分布，或根据具体情况在地形表面或地层中施加不同的压力分布，对公式(17)中的电势项与压力项移至左端，分离常数项得到：

其中，c(k)为关于渗透率的沉积耦合系数，与电势和压力之比有关且恒定等于

μ为地下介质粘度，Q_V(K)对渗透率的依赖关系为：

log₁₀Q_V＝-9.2-0.82log₁₀K (19)

将公式(18)与公式(19)进行联立，并对等式两边取对数得公式(20)：

将公式(20)进一步化简得到公式(21)：

其中等式左端为渗透率K的耦合系数L(k)，通过地表测得的自然电位信号进行反演恢复电势所得，通过公式(20)得到公式(21)求得渗透率分布：

K＝-1.7148×ln[L(k)] (22)。

进一步地，步骤S3具体包括：

S31利用已有的地质资料、地球物理资料、钻井资料研究地下的形态信息，作为地下压力场的初始分布状态，或在信号反演中设为理想地下压力模型；

S32利用高密度或各类电法勘探方法进行地下电导率监测，并作为自然电位信号反演的约束条件；

S33监测自然电位信号，对获取的数据进行预处理分析及实验室内的正反演计算(S1步骤)，最终恢复地下自然电势分布；

S34利用S2步骤计算渗透率耦合系数L(k)；

S35通过化简式估计地下渗透率分布情况。

本发明与现有技术相比，有益效果在于：

(1)本发明对自电位信号进行反演来估计地下渗透率系数，从而获得地下水流的优先流动路径。通过系统的理论推导与模型建立，利用正演模型进行估计方法的可靠性研究，再将反演恢复的渗透信息与正演结果进行对比。

(2)本发明将改进后的推导形式进行数值计算，提高了反演结果的精度，减少了由于计算产生的数量误差。同时在反演过程中保留了原始数据有效信息，提高了反演结果的准确性和完整性。且本发明中提出的方法相对于传统的技术方式在计算效率上有所提高。

附图说明

图1是本发明的方法流程图；

图2为本发明初始模型渗透率分布图；

图3为本发明正反演模型自然电势分布图(a)正演自然电势的分布结果、(b)反演自然电势的分布结果；

图4为本发明模型地表自电位信号；

图5为本发明正反演渗透率分布结果(a)正演模型的渗透率分布结果、(b) 反演的渗透率估计结果；

图6为本发明实测自然电位信号；

图7为本发明实测信号恢复的自然电势结果；

图8为本发明实测电导率分布图；

图9为本发明反演实测信号恢复的地下渗透率分布图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参见图1，一种基于自然电位反演的渗透率估计方法，该方法包括：S1通过建立水力与电性耦合模型模拟地下自然电位信号的正反演响应：

其中正演通过模拟地下水随压力变化的运移程度，再计算耦合电性的连接参数与描述水流流动的速度参数，最后通过电性控制方程计算地下的电流密度源与电导率约束条件还原电势分布；利用第二类边界条件在水力与电性耦合模型的顶部获取自然电位信号，模拟地下水流动激发的电性响应；正演过程的输入参数包括:电导率，渗透率与边界条件，输出参数为电流密度以及其表层接收的自然电位信号；

反演通过建立真实采集信号与正演模拟信号的差值的目标函数，利用最小二乘牛顿迭代或最速下降法使得目标解趋于最优值，其中，正演计算的结果不断迭代来拟合采集信号，来获得最优解状态下的模型响应；通过地表采集的自然电位信号与正演模拟的输出信号进行观测差最小化，获得地下电流密度分布；

S2根据步骤S1获得的地下电流密度分布，利用采集的电导率约束还原地下电势场，并根据电势场求取渗透率分布；

S3根据步骤S1和步骤S2利用多项电性参数与地下流体流动信息对水文观测区域或地热靶区进行评价，通过利用现场获得的高密度电阻率数据进行倒转与处理获得电导率分布，结合自电位信号进行正反演计算恢复地下自然电势场，估计地下渗透率分布信息。

自然电位信号的正演问题是一个耦合场的问题，通过模拟地下水随压力变化的运移程度(水力问题)，再计算耦合电性的连接参数与描述水流流动的速度参数，即为达西速度，最后通过电场的控制方程计算地下的电势分布(电性问题)。这一系列的复杂计算名为正演计算。目的是通过模拟地下的水文渗透率与电阻率分布来尽可能的还原自然电位信号的输出响应。这也是还原问题本质的关键一步。

利用本构方程和连续性来解释非均匀的和各向异性的多孔材料的水力问题：

其中，u表示达西速度(m/s)，a为比存量(m-1)，K为渗透率值(m/s)， Q_s为外部源(s^-1),孔隙水的压力可以表示成水的质量密度ρ_f(kg/m³)，重力加速度g(m/s²)，液压头h(m)的关系表达式p＝ρ_fg(h-z)给出。其中z(m)为以基准点为准的高度。上式要遵循第一类边界条件Γ_D和第二类边界条件Γ_N，同时满足初始条件：

initial h＝h₀att＝0

在各向异性介质中，导电电流密度(根据欧姆定律)

与源电流密度j_s之和为总电流密度j(A/m²)。

j_s＝Q_Vu (3)

其中，σ表示为电导率值(S/m²)，ψ为电势(V)，Q_V是每单位空隙体积的有效过剩电荷密度，通过下述经验公式表达了Q_V的渗透率依赖性：

log₁₀Q_V＝-9.2-0.82log₁₀K (4)

根据电荷的连续性方程

具有以下边界条件:

其中，在上界面(空气-地面界面)施加了第二类边界条件Γ_NV并且在其他边界施加了第一类边界条件Г_DV，这样可以确保地下电势的响应在地表能接收，不向周围底层扩散。上述边界问题表明，异常包含与液压流动路径相关的某些信息。

反演问题与正演问题相反，在知道地表观测的自然电位信号的情况下，建立一个均匀渗透率分布，尽可能的计算其输出后的模拟信号，将模拟信号与真实的自然电位信号之间的误差达到最小，以至于恢复地下的渗透率与电导率，其中通过高密度等方法可以约束电导率的计算，对于反演问题获得更加准确的结果提供帮助。

获得地下源电流密度矢量引起的电势异常空间分布。利用一组点M上的j_s

ψ(N)＝∫_ΩG(N，M)j_s(M)dV (7)

其中j_s(M)为一组点M上的电流密度源，ψ(N)为一组N电极处的电势， G(N，M)表示在点N处测量自然电位数据和当前位于点M处的内核矩阵。 G(N，M)内核函数的计算取决于几种参数的选取：①介质的电导率值②离散元素M的数量。通过利用有限元法可以求解内核矩阵。

其中|| ||₂为L2范数，G表示N×2M的核矩阵(相对应于每个站点测得的自电位值)，ω为与电势和内核矩阵G有关的电流密度值j_s。ψ^obs表示一组自电位数据相对应的N个元素的向量，W_d为N×N的对角线加权平方矩阵，且δ₁为偏差的平方数。因反演会导致求解的不适定性，对分析过程中加入正则化项与深度加权函数φ_m，减少过拟合问题。建立全局目标函数项φ_T：

其中

为N×2M的深度加权矩阵，z₀表示模型单元的观测高度，L是基于ω的一阶导数的光滑度，β为常数项，介于1到4之间。λ表示正则化的约束项0＜λ＜∞，上式目标函数可以用标准形式给出：

其中

L_P是通过基于QR分解的Elédn算法推导出的矩阵。求目标函数φ_T最小值的解决方案：

在地下渗透率参数已知的情况下可计算出原始电流密度分布用于目标函数的先验模型，同样在没有先验模型

的情况下：

为减少核矩阵的计算量使用奇异值分解方法，将

H＝min(N，M)，并重新定义

其中，Λ为奇异对角矩阵，其元素Λ_i表示矩阵对角线上的奇异值分量，uv^T为左右奇异向量。在计算效率上，SVD后的等式更能体现出主要成分的信息，有效的减少了计算量。除此之外，正则化参数λ的选取也是至关重要的，采用了L 曲线方法(φ_d与φ_m的对数-对数交会图)来定义正则化参数的最优值。

S2对渗透率参数的估计计算：

在S1中反演获得地下真实电流密度分布j_s，计算公式(14)的电导率σ约束还原地下电势场ψ结果：

利用本构方程和电性方程得到水力问题与电性问题的解决之后，通过达西速度方程

将公式(15)中的两种势场计算联立：

其中，压力p的空间变化等于液压头的变化

将其带入公式(16) 中获得：

认为压力梯度

的变化是与渗透无关的，在估计渗透率K时设为均匀半空间的理想压力分布，也可根据具体情况在地形表面或地层中施加不同的压力分布(施加不同的边界条件)。对公式(17)中的电势项与压力项移至左端，分离常数项：

其中，c(k)为关于渗透率的沉积耦合系数，与电势和压力之比有关且恒定等于

μ为地下介质粘度，Q_V(K)对渗透率的依赖关系为：

log₁₀Q_V＝-9.2-0.82log₁₀K (19)

将公式(18)与公式(19)进行联立，并对等式两边取对数得公式(20)：

由于计算过程中出现的精度不足等问题，将公式(19)进一步化简得到公式(20)：

其中等式左端为渗透率K的耦合系数L(k)，通过地表测得的自然电位信号进行反演恢复电势所得，通过公式(20)得到公式(21)求得渗透率分布：

K＝-1.7148×ln[L(k)] (22)。

S3利用多项电性参数与地下流体流动信息对水文观测区域或地热靶区进行评价：

通过利用现场获得的高密度电阻率数据进行倒转与处理获得电导率分布，结合自电位信号进行正反演计算恢复地下自然电势场(步骤S1)，估计地下渗透率分布信息(步骤S2).

步骤S3包括：S31利用已有的地质资料、地球物理资料、钻井资料研究地下的形态信息，作为地下压力场的初始分布状态，或在信号反演中设为理想地下压力模型；

S32利用高密度或各类电法勘探方法进行地下电导率监测，并作为自然电位信号反演的约束条件；

S33监测自然电位信号，对获取的数据进行预处理分析及实验室内的正反演计算(S1步骤)，最终恢复地下自然电势分布；

S34利用S2步骤计算渗透率耦合系数L(k)(公式21)；

S35通过化简式估计地下渗透率分布情况(公式22)。

实施例

如图1所示，一种基于自电位信号反演的地下渗透率估计方法，包括如下步骤：

(S1)通过建立水力与电性耦合模型模拟地下自然电位信号的正反演响应：

通过建立水力与电性传输正演模型模拟地下空间自然电位信号对地下水流动的变化响应，之后通过最小二乘反演问题求解计算获得地下电势分布场。

(S2)计算沉积耦合系数，估计地下渗透率分布信息：

首先通过建立模型进行水利与电性问题的解决，在设定好各个边界条件后获得自然电势分布。在考虑无初始压力分布的条件下利用电导率参数进行约束，求取关于渗透率的耦合系数。

log₁₀Q_V＝-9.2-0.82log₁₀K

其中，Δψ为通过泊松方程计算的地下电势分布，Δh为地下压力场分布，ρ为地下空间的密度，g为加速度参数，μ为地层物质的粘度，σ为电导率分布信息 Q_V(K)K为渗透率的变量参数，其中电荷密度Q_V(K)与渗透率K具有线性依赖性 (Jardani and Revil,2009)。

(S3)利用多项电性参数与地下流体流动信息对水文观测区域或地热靶区进行评价：

通过利用现场获得的高密度电阻率数据进行倒转与处理最终获得电导率分布，再结合自电位信号进行反演恢复地下自然电势场，准确估计地下渗透率分布信息。

(S31)利用已有的地质资料、地球物理资料、钻井资料研究地下的形态信息，可以作为地下压力场的初始分布状态，也可以在信号反演中设为理想地下压力模型；

(S32)利用高密度或各类电法勘探方法进行地下电导率监测，并作为SP信号反演的约束条件；

(S33)监测SP信号，对数据进行预处理分析及反演，最终恢复地下自然电势分布；

(S34)计算渗透率耦合系数及适当的带通滤波，以确保数据的有效成分信息的保留；

(S35)通过化简式估计地下渗透率分布情况。

通过利用自电位信号在地热异常区域或深层地热注入和生产井的运作过程中进行测量，可以对地下因热异常而引起的水流流动进行有效的观测。通常，在SP信号幅值高的区域会显示出高渗透率信息。

为了测试本发明实际应用效果，首先通过正演与反演计算分别估计模型渗透率，并验证其可用性，再通过提出的反演算法对实测SP信号进行渗透率估计，进行综合评价。

图2是初始模型渗透率分布图，通过电性与水利方程的计算，在地表与地底施加第一类边界条件，使得水流向更深出渗透，最终获得地下正演自然电势的分布结果(图3(a))，之后提取地表自然电位信号进行反演计算获得由信号(图4)恢复的自然电势分布(图3(b))。分别将两种电势结果与电导率分布作为约束条件，计算出耦合系数L(K)，最终获得正演模型的渗透率分布结果(图5(a))与自电位信号反演的渗透率估计结果(图5(b))。可以看出，通过正演理论计算估计出的渗透率与初始模型给出的渗透率在形状分布上保持一致，并且与周围地层拥有明显的区分性。通过反演自电位信号估计的渗透率与正演模型和初始渗透率结果相一致，数值误差较小，可以有效的刻画出高渗透路径。

图6为在实验测区采集的自然电位数据，通过噪声处理与校正后用于反演来恢复地下自然电势分布(图7)，同时，在测区进行了高密度电法勘测并记录电阻率剖面作为渗透率估算的约束条件(图8)。计算耦合系数L(K)时，考虑到地层的压力场变化，给出多层介质的理想压力分布情况，最终估计出实验测区的渗透率分布结果(图9)。可以看出，在测量剖面上存在多处高渗透性区域，并呈柱状向地底延伸，与实际测量过程中所发现的钻孔井位相对应，拥有良好的估计结果。

通过上述方法，本发明能很好地利用自然电位信号对地下渗透路径进行估计与计算，提高模型精度，以此更准确地圈定地下水流信息与地热流体异常情况。因此，本发明具有突出的实质性的特点和进步。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于自然电位反演的渗透率估计方法，其特征在于，该方法包括：

S1通过建立水力与电性耦合模型模拟地下自然电位信号的正反演响应：

其中正演通过模拟地下水随压力变化的运移程度，再计算电性耦合的连接参数与描述水流流动的速度参数，最后通过电性控制方程计算地下的电流密度源与电导率约束条件还原电势分布；利用第二类边界条件在水力与电性耦合模型的顶部获取自然电位信号，模拟地下水流动激发的电性响应；正演过程的输入参数包括:电导率，渗透率与边界条件，输出参数为电流密度以及其表层接收的自然电位信号；

反演通过建立真实采集信号与正演模拟信号的差值的目标函数，利用最小二乘牛顿迭代或最速下降法使得目标解趋于最优值，其中，正演计算的结果不断迭代来拟合采集信号，来获得最优解状态下的模型响应；通过地表采集的自然电位信号与正演模拟的输出信号进行观测差最小化，获得地下电流密度分布；

S2根据步骤S1获得的地下电流密度分布，利用采集的电导率约束还原地下电势场，并根据电势场求取渗透率分布；

S3根据步骤S1和步骤S2利用多项电性参数与地下水流动信息对水文观测区域或地热靶区进行评价，通过利用现场获得的高密度电阻率数据进行倒转与处理获得电导率分布，结合自然电位信号进行正反演计算恢复地下自然电势场，估计地下渗透率分布信息；

步骤S2利用步骤S1中反演获得地下真实电流密度分布j_s，计算下式(14)的电导率σ约束还原地下电势场ψ结果：

其中，a为比存量，K为渗透率值，Q_V是每单位空隙体积的有效过剩电荷密度，

通过达西速度方程

将公式(15)中的两种势场计算联立得到：

其中，压力p的空间变化等于液压头的变化

将其带入公式(16)中获得：

认为压力梯度

的变化是与渗透无关的，在估计渗透率K时设为均匀半空间的理想压力分布，或根据具体情况在地形表面或地层中施加不同的压力分布，对公式(17)中的电势项与压力项移至左端，分离常数项得到：

其中，c(k)为关于渗透率的沉积耦合系数，与电势和压力之比有关且恒定等于

μ为地下介质粘度，Q_V(K)对渗透率的依赖关系为：

log₁₀Q_V＝-9.2-0.82log₁₀K (19)

将公式(18)与公式(19)进行联立，并对等式两边取对数得公式(20)：

将公式(20)进一步化简得到公式(21)：

其中等式左端为渗透率K的耦合系数L(k)，通过地表测得的自然电位信号进行反演恢复电势所得，通过公式(20)得到公式(21)求得渗透率分布：

K＝-1.7148×ln[L(k)] (22)。

2.按照权利要求1所述的方法，其特征在于，步骤S1中，正演包括利用公式(1)的本构方程和连续性来解释非均匀的和各向异性的多孔材料的运移程度：

其中，u表示达西速度，a为比存量，K为渗透率值，Q_s为外部源，h为液压头；

公式(1)遵循第一类边界条件Γ_D和第二类边界条件Γ_N，同时满足初始条件：

initial h＝h₀ at＝0

在各向异性介质中，导电电流密度

与源电流密度j_s之和为总电流密度j：

j_s＝Q_Vu (3)

其中，σ表示为电导率值，ψ为电势，Q_V是每单位空隙体积的有效过剩电荷密度，通过经验公式(4)表达有效过剩电荷密度Q_V的渗透率依赖性：

log₁₀Q_V＝-9.2-0.82log₁₀K (4)

根据电荷的连续性方程

具有以下边界条件:

其中，σ表示为电导率值，ψ为电势。

3.按照权利要求1所述的方法，其特征在于，步骤S1中，

反演包括：用一组点M上的j_s

ψ(N)＝∫_ΩG(N,M)j_s(M)dV (7)

其中j_s(M)为一组点M上的电流密度源，ψ(N)为一组N电极处的电势，G(N,M)表示在点N处测量自然电位数据和当前位于点M处的内核矩阵，通过利用有限元法求解内核矩阵，

其中‖ ‖₂为L2范数，G表示N×2M的核矩阵，ω为与电势和内核矩阵G有关的参数项，即为电流密度值j_s，ψ^obs表示一组自然电位数据相对应的N个元素的向量，W_d为N×N的对角线加权平方矩阵，且δ_N为偏差的平方数，加入正则化项与深度加权函数φ_m，建立全局目标函数项φ_T：

其中

为N×2M的深度加权矩阵，z₀表示模型单元的观测高度，L是基于ω的一阶导数的光滑度，β为常数项，介于1到4之间，λ表示正则化的约束项0<λ<∞，目标函数用标准形式给出：

其中

L_P是通过基于QR分解的Elédn算法推导出的矩阵，求目标函数φ_T最小值的通过公式(11)求取：

在地下渗透率参数已知的情况下计算出原始电流密度分布用于目标函数，在没有先验模型，即

的情况下得到公式(12)：

为减少核矩阵的计算量使用奇异值分解方法，将

并重新定义

其中，Λ为奇异对角矩阵，其元素Λ_i表示矩阵对角线上的奇异值分量，u为左奇异向量，v^T右奇异向量，

为反演计算后的更新电流密度参数，通过迭代误差的减少获得计算稳定值分布。

4.按照权利要求1所述的方法，其特征在于，步骤S3具体包括：

S31利用已有的地质资料、地球物理资料、钻井资料研究地下的形态信息，作为地下压力场的初始分布状态，或在信号反演中设为理想地下压力模型；

S32利用高密度方法进行地下电导率监测，并作为自然电位信号反演的约束条件；

S33监测自然电位信号，对获取的数据进行预处理分析及实验室内的正反演计算步骤S1，最终恢复地下自然电势分布；

S34利用S2步骤计算渗透率耦合系数L(k)；

S35通过化简式估计地下渗透率分布情况。

Images