基于数据动态控制时间光滑约束的四维电阻率反演方法

摘要

本发明公开了一种基于数据动态控制时间光滑约束的四维电阻率反演方法，将采集过程得到的电阻构成数据集合，定义采集一组完整数据集用时和一次反演解译所需时间；将采集的完整数据集均分为多份，对拆分后的小份数据进行重组，得到重组数据集；定义数据响应因子，构造动态调整时间光滑度的四维电阻率反演的目标函数；根据数据响应因子，选取合适的时间光滑因子，构造动态调整时间光滑度的四维电阻率反演方程，求解模型参数增量向量，迭代反演后，求得最终反演结果。

权利要求

1.一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：包括以下步骤：(1)将采集过程得到的电阻构成数据集合，定义采集一组完整数据集用时和一次反演解译所需时间；

(2)将采集的完整数据集均分为多份，对拆分后的小份数据进行重组，得到重组数据集；

(3)定义数据响应因子，定量的表示相邻两次观测数据变化的大小，构造动态调整时间光滑度的四维电阻率反演的目标函数；

(4)根据数据响应因子，选取合适的时间光滑因子，对反演目标函数求最小值，构造动态调整时间光滑度的四维电阻率反演方程，求解模型参数增量向量，迭代反演后，求得最终反演结果。

2.如权利要求1所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(1)中，对某区域进行监测，连续进行数据采集，使用一定的电极排列，对AB两点供电，测量MN的电位差，除以此时的供电电流，即可获得电阻数据。

3.如权利要求1所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(2)中，对拆分后的小份数据进行重组，将任意相邻的多份数据重新组合为新的完整数据集，称为重组数据集。

4.如权利要求1所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(2)中，取重组数据集用于反演，用DRk表示，其中，k＝1，2，…，Nk，则新的数据集合为DR＝(DR1,DR2,…,DRNk)T，Nk表示用于同时反演的重组数据集的个数。

5.如权利要求1所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(3)中，定义 为当前时刻采集的数据，即采集到第in组完整数据集中第jn个元素；此时的数据响应因子 拆分后每一份数据集的个数为Nr＝dnum/NCS， 为第in组、in-1组采集完整数据集中的部分元素，其中，k＝0，1，…，Nr-1，dnum表示一份完整数据集中所包含的元素个数；数据响应因子为相邻两次相同电极排列的数据的2范数，定量的表示相邻两次观测数据变化的大小；

Random为符合正态分布的随机数，数学期望为0，标准差为σ,定义采集一组完整数据集用时为TC,一次反演解译所需时间为TS,NCS＝[TC/TS],[*]是取整运算,采集的完整数据集均分为NCS小份。

6.如权利要求4所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(4)中，目标函数为：Φ＝||eTe||2+λψ+μΓ

等式右端第一项是数据误差项，e是实际观测数据与理论观测数据的差值向量，第二项是模型光滑约束项，第三项即时间光滑约束项，其中ψ、Γ分别为空间和时间上相应的目标函数，λ和μ用于调节空间和时间模型光滑约束对目标函数影响的权重，时间光滑约束项具体如下：其中Mk是一组模型参数，与DRk一一对应，M＝[M1,...,MNk]T，ΔM是M在迭代中的模型参数增量向量，ωk称为时间光滑因子，用来控制时间光滑度的大小，tk表示不同的时刻，k＝1，

2，…，Nk，tk+1-tk为采集拆分后的1小份数据所用的时间间隔，CT为时间光滑约束矩阵。

7.如权利要求6所述的一种基于数据动态控制时间光滑约束的四维电阻率反演方法，其特征是：所述步骤(4)中，动态调整时间光滑度的四维电阻率反演方程如下：(JTJ+λCTC+μCTTCT)ΔM＝{JT(DR-G(M))-μCTTCTM}其中J为偏导数矩阵，C表示光滑度矩阵，G表示正演过程；

求解反演方程得到ΔM，通过迭代得到反演结果。

说明书

基于数据动态控制时间光滑约束的四维电阻率反演方法

技术领域

[0001] 本发明涉及一种基于数据动态控制时间光滑约束的四维电阻率反演方法。

背景技术

[0002] 四维电阻率层析成像技术可以用来表征在地球表面发生的动态过程。该技术通过捕获与时间有关的电导率变化，来反演和解析所监测过程相关的空间属性变化。近年来，已经广泛的应用于环境监测，工程建设和水文地球物理研究等领域，是一种精细化监测的好
手段。

[0003] 目前的监测系统，可以通过安装永久性的现场传感器和遥测技术，远程启动和控制数据采集，高度自动化的对数据传输和处理。与人工重复测量和手动控制永久性设施相
比，可以显著的提高时间分辨率。而自动化处理和反演需要测量数据集具有一致性和完整
性，故仪器一方面考虑了测量数据与特定数据集相关的固有关系，例如电极排列；另一方面提供了时间推移处理的框架和后续数据集之间相互依赖的分析。它为四维电阻率层析成像
技术提供了一个先决条件，可以使用以前收集的数据来自动地改进获取和解释新数据。因
此优化采集方案和协调连续数据集之间的依赖关系是未来的一个发展重点。

[0004] 目前，电法监测的研究取得了一定的效果，但仍存在以下关键问题，尚未解决：

[0005] ①传统的四维电阻率反演方法应用于监测过程时，通常需要等待一组或多组完整排列的数据采集完毕后再进行反演解译。现行的计算机或服务器计算效率已经得到大幅提
高，一般而言，监测过程中采集数据的耗时大于反演计算耗时，反演计算快速的优势并没有得到最大程度的发挥。

[0006] ②现有的电法监测通常是在时间维度上，建立连续数据集之间的相互依赖关系。这些时间正则化和空间正则化结合的技术，减少了反演伪像产生，但时间正则化可能会使
一个模型的时间域过于平稳，不利于发现瞬时变化。被监测区域地质条件往往是复杂的，地下过程是快速多变的。

发明内容

[0007] 本发明为了解决上述问题，提出了一种基于数据动态控制时间光滑约束的四维电阻率反演方法，本发明更加高效的利用反演计算快速的优势，进一步提高监测质量和捕获
潜在的小规模瞬态事件的能力。通过改进了观测数据的提取方式，以及根据观测数据的变
化大小来调整时间域上的光滑度，大大减少了反演过程中的假异常，提高捕获潜在的小规
模瞬态事件的能力，可以更早的、更加准确的发现异常体的形态和定位。

[0008] 为了实现上述目的，本发明采用如下技术方案：

[0009] 一种基于数据动态控制时间光滑约束的四维电阻率反演方法，包括以下步骤：

[0010] (1)将采集过程得到的电阻构成数据集合，定义采集一组完整数据集用时和一次反演解译所需时间；

[0011] (2)将采集的完整数据集均分为多份，对拆分后的小份数据进行重组，得到重组数据集；

[0012] (3)定义数据响应因子，构造动态调整时间光滑度的四维电阻率反演的目标函数；

[0013] (4)根据数据响应因子，选取合适的时间光滑因子，构造动态调整时间光滑度的四维电阻率反演方程，求解模型参数增量向量，迭代反演后，求得最终反演结果。

[0014] 所述步骤(1)中，对某区域进行监测，连续进行数据采集，使用某种特定的电极排列，如温纳装置、施伦贝尔装置、bipole-bipole(跨孔装置)或某几种装置的组合，对AB两点供电，测量MN的电位差，除以此时的供电电流，即可获得电阻数据。

[0015] 所述步骤(2)中，对拆分后的小份数据进行重组，将任意相邻的多份数据重新组合为新的完整数据集，称为重组数据集。

[0016] 所述步骤(2)中，取重组数据集用于反演，用DRk(k＝1,2,...,Nk)表示，则新的数据集合为DR＝(DR1,DR2,…,DRNk)T，Nk表示用于同时反演的重组数据集的个数。

[0017] 所述步骤(3)中，定义数据响应因子 定量的表示相邻两次观测数据变化的大小。

[0018] 所述步骤(3)中，定义 为当前时刻采集的数据，即采集到第in组完整数据集中第jn个元素。则此时的数据响应因子 拆分后每
一份数据集的个数为Nr＝dnum/NCS， 为第in组、in-1组采集完整
数据集中的部分元素。数据响应因子为相邻两次相同电极排列的数据的2范数，定量的表示相邻两次观测数据变化的大小。这里Random为符合正态分布的随机数，数学期望为0，标准差为σ。

[0019] 所述步骤(4)中，目标函数为：

[0020] Φ＝||eTe||2+λψ+μΓ

[0021] 等式右端第一项是数据误差项，Wd为数据权重矩阵，e是实际观测数据与理论观测数据的差值向量。第二项是模型光滑约束项，第三项即时间光滑约束项。其中ψ、Γ分别为时间和空间上相应的目标函数，λ和μ用于调节空间和时间模型光滑约束对目标函数影响的权重。

[0022] 时间光滑约束项具体如下：

[0023]

[0024] 其中Mk是一组模型参数，与DRk一一对应，M＝[M1,...,MNk]T，ΔM是M在迭代中的模型参数增量向量，ωk称为时间光滑因子，用来控制时间光滑度的大小，tk(k＝1,2,...,Nk)表示不同的时刻，tk+1-tk为采集拆分后的1小份数据所用的时间间隔。CT为时间光滑度矩阵。

[0025] 所述步骤(4)中，动态调整时间光滑度的四维电阻率反演方程如下：

[0026] (JTJ+λCTC+μCTTCT)ΔM＝{JT(DR-G(M))-μCTTCTM}

[0027] 其中J为偏导数矩阵，C表示空间光滑度矩阵，G表示正演过程。

[0028] 求解反演方程得到ΔM，通过迭代得到反演结果。

[0029] 与现有技术相比，本发明的有益效果为：

[0030] 1.本发明提出了一种新型数据提取策略，将监测过程中连续采集的多个完整数据集按一定规则和采样时刻进行拆分和重组，构建与反演计算时长相匹配的新的数据集合，
使得反演程序可以利用具有更小监测时间间隔的完整数据集进行快速不间断反演。一次反
演后，无需等待一组完整数据集的采集，加入新采集的、还未使用的数据进行下一次反演。
能够更快地将观测数据纳入反演方程。充分发挥反演计算快速的优势。

[0031] 2.本发明定义了数据响应因子，可以定量的计算新采集的数据与相同电极排列下的旧数据之间的差异。该方法可以更早的对地下是否发生较大异常做出判断。

[0032] 3.本发明提出了基于观测数据变化特征控制时间光滑约束的电阻率反演方法，根据数据响应因子，即新采集数据与相同电极排列的旧数据的差异大小，确定时间光滑约束
因子的大小。当该数据响应因子值大时，认为地下电阻率结构发生较大变化，此刻减小光滑度，有利于更好的发现瞬时变化。当该数据响应因子值较小时，认为地下电阻率结构发生微小变化，增大光滑度，减少反演伪像的产生。使得整个监测过程，对异常体的捕捉更加及时，定位更加准确。

附图说明

[0033] 构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

[0034] 图1是本发明提出的新型数据提取策略流程图；

[0035] 图2是基于观测数据变化特征控制时间光滑约束的四维电阻率反演流程图；

[0036] 图3是本发明进行数值模拟时使用的地电模型设计图；

[0037] 图4是本发明根据图3的模型设计图采取不同光滑约束得到的四维电阻率反演成像结果；
具体实施方式：

[0038] 下面结合附图与实施例对本发明作进一步说明。

[0039] 应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

[0040] 需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

[0041] 在本发明中，术语如“上”、“下”、“左”、“右”、“前”、“后”、“竖直”、“水平”、“侧”、“底”等指示的方位或位置关系为基于附图所示的方位或位置关系，只是为了便于叙述本发明各部件或元件结构关系而确定的关系词，并非特指本发明中任一部件或元件，不能理解为对本发明的限制。

[0042] 本发明中，术语如“固接”、“相连”、“连接”等应做广义理解，表示可以是固定连接，也可以是一体地连接或可拆卸连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的相关科研或技术人员，可以根据具体情况确定上述术语在本发明中的具体含义，不能理解为对本发明的限制。

[0043] 一种基于数据动态控制时间光滑约束的四维电阻率反演方法，包括以下步骤：

[0044] (1)定义采集的一系列数据集合D＝(D1,D2,…,DNc)T。其中，Di(i＝1,2,…,Nc)表示第i次采集过程得到的电阻数据集。

[0045] (2)定义 其中， 表示数据集Di中每一个元素，dnum表示一份完整数据集中所包含的元素个数。

[0046] (3)定义采集一组完整数据集用时为TC,一次反演解译所需时间为TS。令NCS＝[TC/TS]([*]是取整运算),将采集的完整数据集均分为NCS小份。

[0047] (4)对拆分后的小份数据进行重组，得到重组数据集DR＝(DR1,DR2,…,DRNk)T。

[0048] (5)定义数据响应因子 定量的表示了相邻两次观测数据变化的大小，作为控制时间光滑因子的依据。

[0049] (6)构造动态调整时间光滑度的四维电阻率反演的目标函数。

[0050] (7)根据数据响应因子，选取合适的时间光滑因子ωk。

[0051] (8)构造动态调整时间光滑度的四维电阻率反演方程，求解模型参数增量向量ΔM，迭代反演后，求得最终反演结果。

[0052] 所述步骤(1)中，对某区域进行监测，连续进行数据采集，使用某种特定的电极排列，如温纳装置、施伦贝尔装置、bipole-bipole(跨孔装置)或某几种装置的组合，对AB两点供电，测量MN的电位差，除以此时的供电电流，即可获得电阻数据。

[0053] 所述步骤(4)中，对拆分后的小份数据进行重组，将任意相邻的NCS小份数据重新组合为新的完整数据集，称为重组数据集。取Nk份重组数据集用于反演，用DRk(k＝1,2,...,Nk)表示，则新的数据集合为DR＝(DR1,DR2,…,DRNk)T，Nk表示用于同时反演的重组数据集的个数。

[0054] 所述步骤(5)中，定义 为当前时刻采集的数据，即采集到第in组完整数据集中第jn个元素。则此时的数据响应因子 拆分后每
一份数据集的个数为Nr＝dnum/NCS， 为第in组、in-1组采集完整
数据集中的部分元素。数据响应因子为相邻两次相同电极排列的数据的2范数，定量的表示相邻两次观测数据变化的大小。这里Random为符合正态分布的随机数，数学期望为0，标准差为σ。

[0055] 所述步骤(6)中，目标函数如下所示：

[0056] Φ＝||eTe||2+λψ+μΓ

[0057] 等式右端第一项是数据误差项，Wd为数据权重矩阵，e是实际观测数据与理论观测数据的差值向量。第二项是模型光滑约束项，第三项即时间光滑约束项。其中ψ、Γ分别为空间和时间上相应的目标函数，λ和μ用于调节空间和时间模型光滑约束对目标函数影响的权重。

[0058] 时间光滑约束项具体如下：

[0059]

[0060] 其中Mk是一组模型参数，与DRk一一对应，M＝[M1,...,MNk]T。ΔM是M在迭代中的模型参数增量向量，ωk称为时间光滑因子，用来控制时间光滑度的大小。tk(k＝1,2,...,Nk)表示不同的时刻，tk+1-tk为采集拆分后的1小份数据所用的时间间隔。CT称为时间光滑约束矩阵，具体形式如下：

[0061]

[0062] 所述步骤(7)中，通常的选取方案如下表。需要说明的是，实际时间光滑度的选取还需根据实际情况分析，本专利只是给出了一个可行的方案。

[0063] 表1时间光滑度的选取方案

[0064]

[0065]

[0066] 所述步骤(8)中，动态调整时间光滑度的四维电阻率反演方程如下：

[0067] (JTJ+λCTC+μCTTCT)ΔM＝{JT(DR-G(M))-μCTTCTM}

[0068] 其中J为偏导数矩阵，C表示光滑度矩阵，G表示正演过程。

[0069] 求解反演方程得到ΔM，通过迭代得到反演结果。

[0070] 作为一种典型实施例，(1)对某区域进行监测，连续进行数据采集，使用跨孔装置。本专利不对电法勘探的采集方式进行限定，本算例中使用了跨孔电阻率技术。

[0071] 定义采集的一系列数据集合D＝(D1,D2,…,DNc)T。其中，Di(i＝1,2,…,Nc)表示第i次采集过程得到的电阻数据集。

[0072] (2)定义 其中， 表示数据集Di中每一个元素，dnum表示一份完整数据集中所包含的元素个数。

[0073] (3)定义采集一组完整数据集用时为TC,一次反演解译所需时间为TS.令NCS＝[TC/TS]([*]是取整运算),将采集的完整数据集均分为NCS小份。

[0074] (4)对拆分后的小份数据进行重组，将任意相邻的NCS小份合为一个新的完整数据集，称为重组数据集，取现有的全部重组数据集用于反演，用DRk(k＝1,2,...)表示，则新的数据集合为DR＝(DR1,DR2,…,DRNk)T。Nk表示用于反演的重组数据集的个数。拆分重组过程如附图1所示，图中取NCS＝4。

[0075] (5)定义 为当前时刻采集的数据，即采集到第in组完整数据集中第jn个元素。则此时的数据响应因子 拆分后每一份数据集的
个数为Nr＝dnum/NCS， 为第in组、in-1组采集完整数据集中的部
分元素。数据响应因子为相邻两次相同电极排列的数据的2范数，定量的表示相邻两次观测数据变化的大小。这里Random为符合正态分布的随机数，数学期望为0，标准差为σ＝0.1(本算例中用值，实际选取要考虑实际情况，本专利不对此进行限定)。

[0076] (6)动态调整时间光滑度的四维电阻率反演的目标函数是：Φ＝||eTe||2+λψ+μΓ，等式右端第一项是数据误差项，Wd为数据权重矩阵，e是实际观测数据与理论观测数据的差值向量。第二项是模型光滑约束项，第三项即时间光滑约束项。其中ψ、Γ分别为空间和时间上相应的目标函数，λ和μ用于调节空间和时间模型光滑约束对目标函数影响的权重。

[0077] (7)时间光滑约束项具体如下：

[0078]

[0079] 其中Mk是一组模型参数，M＝[M1,...,MNk]T，与DR一一对应。ΔM是M在迭代中的模型参数增量向量，ωk称为时间光滑因子，用来控制时间光滑度的大小。tk(k＝1,2,...,Nk)表示不同的时刻，tk+1-tk为采集拆分后的1小份数据所用的时间间隔。CT称为时间光滑约束矩阵，具体形式如下：

[0080]

[0081] (8)根据数据响应因子，选取合适的时间光滑因子ωi。通常的选取方案是：当时，令ωi＝0；当 时,令ωi∈(0.005,0.02)；当 时，令ωi＝0.1。
需要说明的是，实际时间光滑度的选取还需根据实际情况分析，本专利只是给出了一个可
行的方案。

[0082] (9)动态调整时间光滑度的四维电阻率反演方程如下：

[0083] (JTJ+λCTC+μCTTCT)ΔM＝{JT(DR-G(M))-μCTTCTM}

[0084] 其中J为偏导数矩阵，C表示光滑度矩阵，G表示正演过程。

[0085] 求解反演方程得到ΔM，通过迭代得到反演结果。

[0086] A.地电模型图如附图3所示，反演区域设置为6m*14m,背景电阻率为800Ω·m。模型中添加了一个瞬时扩散的低阻体，其电阻率为50Ω·m，分别在t4和t5时刻发生突变，共3个阶段，阶段内地下情况均为静态的。各阶段的低阻体在模型中的位置和形态如图3所示。

[0087] B.正演使用矩形单元双线性插值的有限元方法，并采用自然边界条件。网格大小为0.5m*0.5m,电极间距设为0.5米,两根测线均布置28个电极。所用数据为bipole-bipole
(AM-BN)和pole–tripole(A-BMN和A-NMB)两种形式的电极排列。

[0088] C.如附图3所示，从t5时刻起,我们得到了2组重组数据集，代入反演可以得到2个反演模型结果(我们认为后一个反演结果对应于t4时刻的模型)；t6时刻后，反演数据加入了S2阶段对应的数据DS5；t7时刻后，反演数据逐步加入了S3阶段对应的数据DS6,DS7,DS8，到t9时刻,我们可以进行t5～t95个时刻的反演。所用的数据以及 如表2所示：

[0089] 表2不同时刻的数据变化情况

[0090]

[0091]

[0092] D.经过反演处理之后的成像效果如附图4所示，每一行对应不同的时间点，每一列分别对应不同的时间光滑约束因子，我们选取6种时间光滑权重ω∈(0,0.001,0.01,0.1,
0.5,2)，分别计算得到了5*6个反演结果图。模型正演结果与观测数据的拟合值
为模型正演结果,di为观测数据,

[0093] E.附图4第一行对应初始阶段S1，同时反演的两组数据相同，故调节时间光滑度仅产生微小变化。第二行S2的1小部分数据DS5进入反演程序，光滑度ω＝0.001和ω＝0.01反演结果中得到了S2的雏形，其它光滑度下均未得到，且ω＝0.1图像中异常体支离破碎。同样的情形在第3,4行同样出现，此时S3的对应的正演数据DS6加入反演，3种模型对应的正演数据混合起来同时反演，ω＝0.01仍然较好的表示出S3模型的形态和位置，但也有少量假
异常。数据的复杂程度高，出现假异常也不可避免，但数据加权促使反演结果倾向于模型
S3，调整时间光滑度使得反演结果更清晰。第5行中的反演结果均可以很好的描述异常体的位置和形态，但ω＝0.1仍有明显的假异常。不适当的时间光滑度权重会使得反演方程更加病态化。纵向来看，第3行的反演结果就逼近S3模型，第五行得到稳定结果。从数据变化幅度看，第二行 的值为2.51，后三行约为4，实质上是模型S2、S3与S1正演结果的数据变化幅
值。RMS值侧面反映了混乱的数据带来的影响，一方面是反演结果中有部分假异常的存在，另一方面，一个模型的正演结果不能等同于3个模型正演结果的混合，所以10％～20％的值就不奇怪了。

[0094] 该数值模拟表面，基于数据动态控制时间光滑约束的四维电阻率反演方法，可以大大减少反演过程中的假异常，提高了捕获潜在的小规模瞬态事件的能力，可以更早的、更加准确的发现异常体的形态和定位。

[0095] 以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

[0096] 上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。