一种基于陆地与水下直流电联合测量的地质反演方法

1.本发明属于陆地与水下直流电探测技术领域，尤其涉及一种基于陆地与水 下直流电联合测量的地质反演方法。


背景技术：

2.迄今为止，直流电的正反演方法与理论较为成熟，地面直流电法应用广， 相对稳定。但是，纵观目前的进展，将陆地与水下直流电联合测量与反演的方 法研究很少，技术还不完备。水库、水电工程坝体安全检测受限于常规方法的 局限亟需突破，开展考虑动水环境的陆地与水下直流电联合测量与反演方法研 究尤为重要。


技术实现要素：

3.本发明提供一种基于陆地与水下直流电联合测量的地质反演方法，旨在解 决上述存在的问题。
4.本发明是这样实现的，一种基于陆地与水下直流电联合测量的地质反演方 法，包括以下步骤：
5.s1：在目标地域中排列布置电极，所述目标地域包括水域，所述水域两侧 分别为陆地，电极排列从一侧陆地逐步向水域底部布置，直至布置到另一侧陆 地上；
6.s2：采用声呐探测装置，确定水域中水体的空间分布，并采集水域不同深 度及不同点位的水样进行分析，确定水体的电阻率空间分布结构，建立水域中 水体及水底地形起伏三维模型；
7.s3：根据陆地及水域底部的电极布置，采用温纳装置或偶极-偶极装置变换 供电电极与接收电极位置，建立不同供电-接收观测系统，获取观测数据；
8.s4：对观测数据进行预处理，去除观测数据中去掉不正常的数据点；
9.s5：根据预处理后的观测数据，获取视电阻率拟断面图，分析视电阻率拟 断面图，初步判断目标地域中地下地质情况；
10.s6：根据目标地域中水体及地质特点建立不同地质模型，进行二、三维正 演模拟计算分析，其中，整个正演模拟计算区域采用非结构化三角单元或四面 体单元进行二、三维网格剖分；
11.s7：根据水域中水体及水底地形起伏三维模型，建立包含三维水体电阻率 结构的初始模型，并使用非结构化三角单元或四面体单元进行网格剖分，以实 现二、三维反演。
12.进一步的，在步骤s3中，采用温纳装置或偶极-偶极装置变换供电电极与 接收电极位置具体包括：采用电极a、b、m、n四个电极，电极a为正极， 电极b为负极，电极m和电极n为测点电极；将电极a、电极b、电极m和 电极n不断地变换位置，从而向地下不同深度输入电流，并获得当前测点电极 m、n的电位差；记录每次变换的电极a、电极b、电极m和电极n位置、电 流及电压，最终视电阻率按以下方式计算：
[0013][0014]
其中，δu
mn
为m、n测量电极对之间的电位差，i为测量电流，k为装置 系数，温纳装置、偶极-偶极装置的装置系数k具体计算法则分别为：
[0015][0016][0017]
其中，am、an、bm、bn、mn分别为各电极之间的距离。
[0018]
进一步的，在步骤s4中，对观测数据进行预处理具体包括：首先对观测数 据进行傅里叶变换，分析动水对电极位置的扰动以及数据的影响，并对同一个 观测点随时间的变化数据进行滤波处理；之后做傅里叶反变换，叠加观测点随 时间变化的数据，并求均值与方差，进一步整理为视电阻率数据，对各点的视 电阻率数据进行编辑，去掉超出正常范围的数据。
[0019]
进一步的，在步骤s5中，获取视电阻率拟断面图具体包括：对于温纳装置 而言，以电极a和电极b之间距离的一半为纵轴，电极m和电极n之间中心 点为测点横向位置，获取视电阻率拟断面图；对偶极-偶极装置而言，电极a、 b之间中点与电极mn之间中点的距离一半为纵轴，测量装置中心点为测点横 向位置，获取视电阻率拟断面图。
[0020]
进一步的，在步骤s6中，进行二、三维正演模拟计算分析具体包括：设置 水体均匀电阻率，水底以下不规则异常体电阻率，其余部分为特定的背景电阻 率；
[0021]
二维与三维电位计算方程分别如下：
[0022]
‑▽
[σ(x,z)
▽
φ(x,y,z)]＝i(x,y,z)
[0023]
‑▽
[σ(x,y,z)
▽
φ(x,y,z)]＝i(x,y,z)
[0024]
其中，σ(x,z)为电导率的二维分布，σ(x,y,z)为电导率的三维分布，
▽
φ(x,y,z) 为电位的三维分布，i(x,y,z)为一个电流源项；
[0025]
利用有限单元法离散上式，得到一个线性方程组：
[0026]
kv＝s
[0027]
其中，k为稀疏矩阵，包含网格单元几何信息以及电阻率结构信息，v为 待求解向量，s为右端项；
[0028]
进一步利用向量v以及插值矩阵，获得计算空间内任意测点位置的视电阻 率值。
[0029]
进一步的，在步骤s7中，具体包括：
[0030]
设定反演目标函数：
[0031][0032]
其中，m为所有电阻率参数组成的列向量，λ为正则化因子，为数据 目标函数，为模型目标函数；
[0033][0034]
其中，δd为观测数据与正演数据的残差向量，wd为数据加权矩阵；
[0035]
wd＝diag{1/σ1,1/σ2,...1/σj...1/σm}
[0036]
其中，数据加权矩阵为一个对角矩阵，对角元素为每个数据方差的倒数；
[0037][0038]
其中，wm为模型协方差矩阵，用于控制模型在横向和纵向的光滑程度， 该矩阵使用梯度算子或拉普拉斯算子对局部模型约束。
[0039]
进一步的，在步骤s7中，反演方法采用高斯牛顿法，具体包括：
[0040]
s71：选择初始电阻率模型m0，设定反演最大迭代次数，以及最小拟合差；
[0041]
s72：求解线性方程组获得迭代的更新方向pk；
[0042]
s73：设定一个初始步长，并开始一维线性搜索，寻优当前迭代的模型更新 步长α，并按如下公式更新模型：
[0043]mk 1
＝mk αpk[0044]
s74：若当前拟合差小于最小阀值，或迭代次数已经大于最大迭代次数，则 终止反演；若上述条件都未满足，则跳至第二步，开始新的迭代；
[0045]
在线性方程组中，包含如下灵敏度矩阵：
[0046][0047]
其中，灵敏度矩阵的元素为不同观测数据关于模型参数的偏导数，该矩阵 为一个满秩矩阵；
[0048]
反演流程中，拟合差按如下计算：
[0049][0050]
进一步的，在反演过程中，λ通过以下机制调整：
[0051]
当反演拟合差难以下降时，λ缩小为原来的十分之一；
[0052]
当线性搜索难以找到合适的更新步长时，λ缩小为原来的十分之一。
[0053]
与现有技术相比，本发明的有益效果是：本发明基于直流电联合测量，从 实际角度出发，设计了多个技术要点，包括：从陆地到水中的递进式测量方法， 有效控制岸边到深水区的地质情况；消除动水对电极扰动的影响；建立水库模 型，消除水体对地下介质反演的影响。
附图说明
[0054]
图1为本发明提供的陆地-水下-陆地电极整体布置示意图；
[0055]
图2为本发明提供的实施方案流程图；
[0056]
图3为本发明实施例提供的电极布置示意图；
[0057]
图4为本发明实施例提供的声呐探测示意图；
[0058]
图5为本发明实施例提供的供电电极与测量电极移动示意图；
[0059]
图6为本发明实施例提供的电阻率正演模型；
[0060]
图7为本发明实施例提供的正演网格剖分示意图；
[0061]
图8为本发明实施例提供的视电阻率拟断面图；
[0062]
图9为本发明实施例提供的反演结果一；
[0063]
图10为本发明实施例提供的反演网格剖分；
[0064]
图11为本发明实施例提供的反演结果二；
[0065]
图12为本发明实施例提供的反演结果三。
具体实施方式
[0066]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实 施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅 仅用以解释本发明，并不用于限定本发明。
[0067]
在本发明的描述中，需要理解的是，术语“长度”、“宽度”、“上”、
ꢀ“
下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、
ꢀ“
底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置 关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或 元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发 明的限制。此外，在本发明的描述中，“多个”的含义是两个或两个以上，除 非另有明确具体的限定。
[0068]
实施例1
[0069]
如图1所示，为实现从岸边逐渐向深水处过度，从而达到有效控制水体下 方及周围相关地质体的目的，本发明首先要求电极排列从陆地逐步向水中布置， 并再次布置到陆地上；
[0070]
为降低水体对反演的影响，本发明要求在设计观测系统及反演前，需要采 用声呐装置，控制当前整个工区的水体及水底地形起伏三维模型，并采集不同 深度及不同点位的水样品，分析水体的电阻率三维分布，以及是否存在明显的 电阻率宏观分层。
[0071]
按照图1所示的陆地及水下电阻率布置，可以采用温纳装置或偶极-偶极装 置等变换供电电极与接收电极位置，建立不同供电-接收观测系统，从而达到更 精准的地质识别能力。
[0072]
由于水流及地下电流的影响对数据的影响，必须在观测数据中去掉不正常 的数据点。首先进行傅里叶变换，分析动水对电极位置的扰动以及数据的影响， 并对同一个观测点随时间的变化数据进行滤波处理。之后做傅里叶反变换，叠 加观测点随时间变化的数据，并求均值与方差，进一步整理为视电阻率数据， 对各点的视电阻率数据进行编辑，去掉超出正常范围的数据。
[0073]
对温纳装置而言，以电极a和电极b之间距离的一半为纵轴，电极m和 电极n之间中心点为测点横向位置，获取视电阻率拟断面图；对偶极-偶极装置 而言，电极a、b之间中点与电极m、n之间中点的距离一半为纵轴，测量装 置中心点为测点横向位置，获取视电阻率拟断面图。根据成像结果，切出不同 位置的水平切片与垂直切片，可初步判断渗漏点位置。
并为正演模拟分析与反 演计算，提供初始模型，减小反演多解性。
[0074]
为适应复杂地形，本发明要求使用三角形或四面体非结构化单元进行网格 剖分。此外，根据当前工区水体及地质特点建立不同地质模型，进行二、三维 正演模拟分析，初步判断当前观测系统的实际分辨能力。
[0075]
为消除水体对反演的影响，本发明要求根据水体及水底三维结构，建立包 含三维水体电阻率结构的初始模型，并使用三角形或四面体非结构化单元进行 网格剖分，以实现二、三维反演。反演前标记水体电阻率结构及地形起伏，在 反演中作为固定参数，不参与模型修改。
[0076]
本发明基于直流电联合测量，从实际角度出发，设计了多个技术要点，包 括：从陆地到水中的递进式测量方法，有效控制岸边到深水区的地质情况；消 除动水对电极扰动的影响；建立水库模型，消除水体对地下介质反演的影响。
[0077]
实施例2
[0078]
如图2所示，本发明公开一种基于陆地与水下直流电联合测量的地质反演 方法，包括以下步骤：
[0079]
1、陆地-水下-陆地电极布置
[0080]
如图3所示，图中纵向z代表高程，横向x代表水平距离。需要将电极从 陆地布置到水下，再布置到陆地上，水下电极与水底面接触，水深约1.5米， 电极从左岸度过水下，一直布设到右岸。
[0081]
2、建立水体模型
[0082]
如图4所示，根据声呐探测结果，确立水体空间分布；根据gps测定，获 取实际水面、水体的高程；根据水区不同点位和深度的水样分析，确定水体总 体的电阻率空间分布结构。
[0083]
若水体存在明显的宏观电阻率分层现象，还需建立具有分层结构的水体电 阻率结构，以便更精确的在正反演中去除水的影响。
[0084]
3、建立观测系统
[0085]
按照图1所示的陆地及水下电阻率布置，可以采用温纳装置、偶极-偶极装 置等变换供电电极与接收电极位置，建立不同供电-接收观测系统，从而达到更 精准的地质识别能力。
[0086]
如图5所示，将温纳装置应用于地质探测，温纳装置的用法是将正极a和 负极b，以及测点电极m、n不断的变换位置，从而向地下不同深度输入电流， 并获得当前mn电位差。每次变换的a、b、m、n位置、电流及电压需要记 录，最终视电阻率按以下方式计算：
[0087][0088]
其中，δu
mn
为m、n测量电极对之间的电位差，i为测量电流，k为装置 系数，具体计算法则为：
[0089][0090]
其中，am、an、mn分别为各电极之间的距离。
[0091]
4、数据预处理
[0092]
去掉灵敏度较大的正演数据，以稳定下一步的反演计算。若处理数据为实 测数据，必须分析实测数据中不正常的飞点，并对实测的电压、电流进行叠加 处理。通过频谱分析及滤波，可去掉流动水以及其它外来电流对观测数据的影 响。
[0093]
此外，通过多次观测，还需要计算当前数据点的方差，方差越大的数据在 反演中权重越小，以此降低质量差的数据在反演中的作用。
[0094]
5、定性分析
[0095]
布极方式为温纳装置，以正极a( )和负极b(-)之间的距离一半为纵 轴，测量电极对m和n之间中心点为测点横向位置，获取视电阻率拟断面图， 视电阻率以上述公式为计算依据。
[0096]
通过分析视电阻率拟断面图，可以初步判断地下地质情况，尤其是潜在的 渗漏点横向位置。
[0097]
6、网格剖分与正演模拟计算
[0098]
正演模拟既是反演计算中重要的组成部分，同时也是实际问题的关键分析 过程。因此步骤6中需要实现网格剖分与正演模拟计算。
[0099]
图6是实施例提供的正演模型图，设置水体均匀电阻率为10ωm，水底以 下不规则异常体电阻率为20ωm，其余部分背景电阻率为100ωm，电极分布与 图3一致。
[0100]
如图7所示，整个正演计算区域采用非结构化三角单元剖分，三角单元具 有较强的适应复杂地形的能力，可提高正反演的计算精度。
[0101]
电位计算方程如下：
[0102]
‑▽
[σ(x,z)
▽
φ(x,y,z)]＝i(x,y,z)
[0103]
其中，σ(x,z)为电导率的二维分布，
▽
φ(x,y,z)为电位的三维分布，i(x,y,z)为 一个电流源项；
[0104]
利用有限单元法离散上式，得到一个线性方程组：
[0105]
kv＝s
[0106]
其中，k为稀疏矩阵，包含网格单元几何信息以及电阻率结构信息，v为 待求解向量，s为右端项；
[0107]
线性方程组可以用直接解或迭代解求解，直接解精度较高，当待求问题较 小时，直接解更快速。迭代解适用于大规模计算问题，内存消耗小。
[0108]
进一步利用最终的解向量v以及插值矩阵，即可以获得计算空间内任意测 点位置的视电阻率值。
[0109]
图8为图6模型最终的视电阻率拟断面图，从图中可以看到视电阻率的分 布可以在一定程度上反应异常体的存在。
[0110]
7、反演计算
[0111]
采用非结构化三角单元剖分整个计算区域，对水体及陆地结构离散化，并 可以设置不同模型区域的电阻率值，以及电极总个数、电极之间的几何距离、 电极入水深度、电极装置类型、视电阻率总个数、模型尺寸、反演最大迭代次 数。
[0112]
本发明的反演计算首先需要设定目标函数，实施例中的反演目标函数为：
[0113]
[0114]
其中，m为所有电阻率参数组成的列向量，λ为正则化因子，为数据 目标函数，为模型目标函数；
[0115][0116]
其中，δd为观测数据与正演数据的残差向量，观测数据与正演数据都为对 数化后的视电阻率值，wd为数据加权矩阵；
[0117]
wd＝diag{1/σ1,1/σ2,...1/σj...1/σm}
[0118]
其中，数据加权矩阵为一个对角矩阵，对角元素为每个数据方差的倒数； 其作用在于控制不同数据在反演中的权重，提高高质量数据的影响，降低低质 量数据的作用。此外，数据加权矩阵还可以起到归一化的作用，对稳定反演有 非常大的帮助。
[0119][0120]
其中，wm为模型协方差矩阵，用于控制模型在横向和纵向的光滑程度， 该矩阵使用梯度算子或拉普拉斯算子对局部模型约束。
[0121]
λ为正则化因子，其作用为控制与之间的关系。λ值越大，则反 演电阻率结构越光滑，分辨率越低；λ越小，反演结果分辨率越高，但反演中 数据项作用强于模型项，可能导致整体反演的多解性增强，反演出现冗余的错 误结构。针对实测数据的反演，λ可选用经验值0.01至10，在反演前可以通过 人为设定来验证反演结果的可靠性。
[0122]
在反演过程中，λ通过以下机制调整：
[0123]
当反演拟合差难以下降时，λ缩小为原来的十分之一；
[0124]
当线性搜索难以找到合适的更新步长时，λ缩小为原来的十分之一。
[0125]
在反演开始前，还需要设定初始模型，越靠近真实模型的初始模型，越能 使反演获得更精确的结果。若反演数据为实测数据，则需要根据测井，地质钻 探等，获取水体以及围岩电阻率值，从而控制已知结构的电阻率分布，达到精 准反演的要求。
[0126]
本发明选用的反演方法具体为高斯牛顿法，高斯牛顿法流程如下：：
[0127]
s71：选择初始电阻率模型m0，设定反演最大迭代次数，以及最小拟合差；
[0128]
s72：求解线性方程组获得迭代的更新方向pk；
[0129]
s73：设定一个初始步长，并开始一维线性搜索，寻优当前迭代的模型更新 步长α，并按如下公式更新模型：
[0130]mk 1
＝mk αpk[0131]
s74：若当前拟合差小于最小阀值，或迭代次数已经大于最大迭代次数，则 终止反演；若上述条件都未满足，则跳至第二步，开始新的迭代；
[0132]
需要说明的是，步长寻优将当前反演退化为单变量寻优问题。此外，在线 性方程组中，包含如下灵敏度矩阵：
[0133][0134]
其中，灵敏度矩阵的元素为不同观测数据关于模型参数的偏导数，该矩阵 为一个满秩矩阵；当求解问题较大时，该矩阵求解慢。为此，本发明使用互易 定理求解灵敏度矩阵，以加快反演计算。
[0135]
反演流程中，拟合差按如下计算：
[0136][0137]
基于以上反演过程及步骤，在图6模型的正演数据基础上加入2％随机噪 音作为观测数据，验证本发明的反演效果。
[0138]
首先将水体作为自由参数进行反演，且水体与陆地结构都从100ωm开始 反演，反演结果如图9所示，反演网格剖分如图10所示。从图9的结果可以看 到，即使是水下探测，若将水体作为自由参数反演且初始电阻率与陆地一致， 难以将水底下的低阻异常体探明。这是因为，水体相对于地下异常体而言，吸 引人工电流场的作用非常强烈。因此，将水体作为自由参数参与反演，会压制 水下异常，无法完成探查水底地质结构的目的。
[0139]
为了避免上述问题，将水体作为已知结构参与反演，即根据声呐探测结果 获取水体结构分布，但依然将水体作为自由参数，在反演中不固定，最终的反 演结果如图11所示。从图中的结果可以看到，尽管水体已是已知值，但作为自 由参数参与计算，仍然难以识别水底低阻异常体。
[0140]
为进一步提高反演效果，第三个反演计算，将水体作为已知结构并全程作 为固定参数不进行修改。图12为再次反演的结果，从结果可以看到，将水体作 为已知且固定的参数参与反演，对水底下方的低阻识别有了明显改善，反演结 果显示在横向10m至20m位置有较大的异常体，这与图6所示水下异常体位置 基本对应。
[0141]
以上仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的 精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保 护范围之内。