CN105373673A - 一种自然电场监测数据动态反演方法及系统 - Google Patents

Abstract

本发明提供了一种自然电场监测数据动态反演方法及系统，该方法包括如下步骤：建立污染监测地的状态模型；建立污染监测地的观测模型；获取第k-1时刻的先验估计污染物扩散演化模型数据；根据状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；获取第k时刻的实测自然电场数据，并根据状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。本发明将污染物扩散演化规律与自然电场实测数据结合，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。

Description

一种自然电场监测数据动态反演方法及系统

技术领域

本发明涉及环境监测领域，尤其涉及一种自然电场监测数据动态反演方法及系统。

背景技术

随着我国经济的持续高速发展，在工业化、城镇化过程中的环境污染问题也日益严重。特别是越来越多的生活和工业垃圾填埋场、工矿废水废渣库等，严重威胁周边的土壤与地下水安全。土壤与地下水等浅地表污染具有隐蔽性和滞后性，还具有累积效应和长期残留影响性，治理非常困难，即时监测即时控制污染源的扩散是保护土壤和地下水不受污染或少受污染的积极方法。与钻孔测试、定期采样分析等传统环境监测手段相比，地球物理方法具有效率高、成本低、快速无损、无二次污染等优点。其中自然电场法野外观测更加快捷便利，且对多种与污染扩散伴生的渗流运动、氧化还原以及微生物活动信号相当敏感，特别适合地下污染物监测。

对于监测采集到的自然电场时序数据的反演解译，目前常规的做法是对各时刻的观测数据单独进行反演处理。这种对时序数据进行静态“单帧”处理的方式可能会带来明显的反演偏差。这是由于污染区中存在多种成因的自然电场与观测数据的关系是非线性的，并且缺乏必要的先验信息对反演进行有效约束。常规反演方法只对某个时刻的数据孤立地进行反演解释，没有把相邻时刻的信息关联起来整体进行分析，缺乏先验知识导致反演结果可靠性低或反演算法不稳定。

现有技术中，地下介质中的自然电场分布满足泊松方程：

其中，σ是污染监测地的地下介质电导率，是污染监测地的自然电场电势分布，j是污染监测地的自然电场电流密度分布。

在污染物监测中，需要通过地表观测到的电势来反演计算地下介质的电导率分布。常规的反演计算就是通过一定的策略来估计符合要求的模型，反演问题即为寻找适当的模型m以满足：

其中，F(m)为模型m的正演计算数据，ε为允许的拟合差。由于正演模型与实际物理模型的差异以及观测数据的噪音影响等因素，这类问题往往没有确定解或存在多解。地球物理反演问题的不适定性在引入Tikhonov正则化后，可以得到适当的改善，求解稳定。

其中，λ为正则化因子，由约束系数矩阵C以及模型m和参考模型m⁰之间的差构成：

最终可以获得反演问题的最小二乘迭代计算方程：

A^TΔd＝(A^TA+λC^TC)Δm，

其中，A为灵敏度矩阵，Δd为观测数据d_ob与正演计算数据f(m)之间的残差，Δm为每次迭代的模型增量。由于每次迭代均需要进行一次正演计算，一次反演过程需要反复多次进行正演计算。如图5a所示，对于监测时序数据的反演，目前常规的做法是对各时刻的观测数据单独进行反演处理，以获得各时刻的静态“单帧”模型。

发明内容

针对现有技术中的缺陷，本发明提供一种自然电场监测数据动态反演方法及系统，用于实现将污染物扩散演化规律与自然电场实测数据结合，为监测地下污染扩散情况提供动态反演结果。

一方面，本发明提供一种自然电场监测数据动态反演方法，所述方法包括：

根据孔隙介质流体运动规律，建立污染监测地的状态模型，所述状态模型表示第k时刻和第k-1时刻先验估计污染物扩散演化模型数据的关系；

根据自然电场正演模拟方法，建立污染监测地的观测模型，所述观测模型表示第k时刻自然电场实测数据与第k时刻先验估计污染物扩散演化模型数据的关系；

初始化自然电场正演模型，获取初始化先验估计污染物扩散演化模型数据M₀，并根据所述状态模型，获取第k-1时刻的先验估计污染物扩散演化模型数据；

根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；

获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。

优选地，

在建立状态模型之前，所述方法还包括：

根据孔隙介质流体运动规律，通过求解理查德方程，建立污染监测地的污染物扩散演化模型；

在建立观测模型之前，所述方法还包括：

根据自然电场正演模拟方法，通过求解泊松方程，建立污染监测地的自然电场正演模型。

优选地，根据孔隙介质流体运动规律，建立污染监测地动态系统的状态模型，所述状态模型公式如下：

M_k＝H(M_k-1)+w_k，

其中，M_k为第k时刻的先验估计污染物扩散演化模型数据，M_k-1为第k-1时刻的先验估计污染物扩散演化模型数据，H为污染监测地的污染物扩散演化模型，w_k为演化误差或噪声。

优选地，根据自然电场正演模拟方法，建立污染监测地的观测模型，所述观测模型公式如下：

Z_k＝F(M_k)+v_k，

其中，Z_k为第k时刻的实测自然电场数据，F为污染监测地的自然电场正演模型，v_k为观测误差或噪声。

优选地，所述根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据，包括：

获取第k-1时刻的先验估计污染物扩散演化模型数据和先验估计协方差P_k-1|k-1；

根据卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1，所述卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=H_k{\hat{M}}_{k-1|k-1},$

$P_{k|k-1}=H_k{P}_{k-1|k-1}{H}_k^T+Q_k,$

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，Q_k为演化过程协方差矩阵。

优选地，所述获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据，包括：

获取第k时刻的实测自然电场数据Z_k；

根据卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k，所述卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k{\tilde{y}}_k,$

$P_{k|k}=(I-K_k{F}_k)P_{k|k-1},$

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，为实测自然电场数据的数据残差，K_k为最优卡尔曼增益，

${\tilde{y}}_k=Z_k-F_k{\hat{M}}_{k|k-1},$

$K_k=P_{k|k-1}{F}_k^T{S}_k^{-1},$

其中，S_k为实测自然电场数据的数据残差协方差，

$S_k=F_k{P}_{k|k-1}{F}_k^T+R_k,$

其中，F_k为第k时刻的自然电场正演模型F，R_k为观测过程协方差矩阵。

优选地，所述根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据，还包括：

设置预设加权点，对污染物扩散演化模型和自然电场正演模型进行无迹变换；

获取每个加权点第k-1时刻的先验估计污染物扩散演化模型数据 ${\hat{M}}_{i,k-1|k-1;}$

根据基于无迹变换的卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1，所述基于无迹变换的卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i{\hat{M}}_{i,k-1|k-1},$

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}}^T,$

其中，i为加权点的编号，2n为加权点的个数，为第i个加权点第k时刻的先验估计污染物扩散演化模型数据，W_i为每个加权点的权重系数。

优选地，所述获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据，包括：

获取第k时刻的实测自然电场数据Z_k；

根据基于无迹变换的卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k，所述基于无迹变换的卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k(Z_k-{\hat{M}}_{k|k-1}),$

$P_{k|k}=P_{k|k-1}-K_k{P}_{z_k{z}_k}{K_k}^T,$

其中，K_k为最优卡尔曼增益，为过程参数，

$K_k=P_{M_k{Z}_k}{P}_{Z_k{Z}_k}^{-1},$

$P_{Z_k{Z}_k}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{Z_{i,k}-{\hat{z}}_{k|k-1}\}{\{Z_{i,k-1|k-1}-{\hat{z}}_{k|k-1}\}}^T,$

$P_{M_k{Z}_k}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{Z_{i,k}-{\hat{z}}_{k|k-1}\}}^T,$

${\hat{z}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i{Z}_{i,k},$

其中，为过程参数，Z_i,k为第i个加权点第k时刻的实测自然电场数据。

另一方面，本发明提供一种自然电场监测数据动态反演系统，所述系统包括：

模型单元，用于根据孔隙介质流体运动规律，建立污染监测地的状态模型，并根据自然电场正演模拟方法，建立污染监测地的观测模型；

反演单元，用于根据初始化先验估计污染物扩散演化模型数据M₀，通过所述模型单元建立的状态模型和观测模型，获取第k时刻的先验估计污染物扩散演化模型数据，并根据第k时刻的实测自然电场数据，对所述第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据；

所述状态模型表示第k时刻和第k-1时刻先验估计污染物扩散演化模型数据的关系，所述观测模型表示第k时刻自然电场实测数据与第k时刻先验估计污染物扩散演化模型数据的关系。

优选地，所述反演单元包括：

预测模块，用于根据第k-1时刻的先验估计污染物扩散演化模型数据，通过所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；

更新模块，用于根据第k时刻的实测自然电场数据，通过所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。

由上述技术方案可知，本发明提供了一种自然电场监测数据动态反演方法及系统，通过建立污染监测地的状态模型和观测模型，利用卡尔曼滤波递归估计模型或基于无迹变换的卡尔曼滤波递归估计模型，实现数据监测的预测过程和更新过程。其中，预测过程是通过第k-1时刻的先验估计污染物扩散演化模型数据，根据状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；更新过程是通过第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。本发明将污染物扩散演化规律与自然电场实测数据结合，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。

附图说明

为了更清楚地说明本公开实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些图获得其他的附图。

图1为本发明一实施例提供的一种自然电场监测数据动态反演方法的流程示意图；

图2为本发明另一实施例提供的一种自然电场监测数据动态反演方法的流程示意图；

图3为本发明另一实施例提供的一种自然电场监测数据动态反演方法的流程示意图；

图4为本发明一实施例提供的一种自然电场监测数据动态反演系统的结构示意图。

图5a为现有技术中传统“孤立”反演模式图；

图5b为本发明的模型演化与观测数据结合反演模式图；

图6为本发明基于卡尔曼滤波的自然电场时序数据动态反演流程图；

图7为本发明基于无迹变换的污染监测自然电场时序数据动态反演流程；

图8为本发明动态反演结果示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1示出了本发明一实施例提供的一种自然电场监测数据动态反演方法的流程示意图，如图1所示，一种自然电场监测数据动态反演方法，所述方法包括：

S11、根据孔隙介质流体运动规律，建立污染监测地的状态模型。

可以理解的是，在建立状态模型之前，应当建立污染监测地的污染物扩散演化模型。本实施例中，根据孔隙介质流体运动规律，通过求解理查德方程，建立污染监测地的污染物扩散演化模型。非均匀介质中的理查德方程一般难以求解或者没有解析解，现有技术中，通常利用数值计算方法如有限元法等由计算机来计算实现，所述理查德方程公式如下：

$\mathrm{\φ}\frac{\∂S}{\∂t}-\▿\[\frac{K\left(S\right)}{{\mathrm{\ρ}}_{c}g}\▿P_c\left(S\right)+K\left(S\right)\hat{Z}\]=0,$

可以理解的是，污染监测地的污染物扩散过程是一个动态系统过程，第k时刻的先验估计污染物扩散演化模型数据可以通过第k-1时刻的先验估计污染物扩散演化模型数据按照状态模型演化获取，所述状态模型公式如下：

M_k＝H(M_k-1)+w_k，w_k～N(0～Q_k)

其中，M_k为第k时刻的先验估计污染物扩散演化模型数据，M_k-1为第k-1时刻的先验估计污染物扩散演化模型数据，H为污染监测地的污染物扩散演化模型，w_k为演化误差或噪声，w_k服从协方差为Q_k的零均值正态分布。

S12、根据自然电场正演模拟方法，建立污染监测地的观测模型。

可以理解的是，在建立观测模型之前，应当建立污染监测地的自然电场正演模型。本实施例中，根据自然电场正演模拟方法，通过求解泊松方程，建立污染监测地的自然电场正演模型。非均匀介质中的泊松方程一般很难直接求解或没有解析解，通常利用数值计算方法如有限元法等由计算机来计算实现，所述泊松方程公式如下：

其中，σ是污染监测地的地下介质电导率，是污染监测地的自然电场电势分布，j是污染监测地的自然电场电流密度分布。

可以理解的是，污染监测地的地下介质电性结构与污染监测地的自然电场分布具有密切关系，通过第k时刻的实测自然电场数据反演获取第k时刻的先验估计污染物扩散演化模型数据，所述观测模型公式如下：

Z_k＝F(M_k)+v_k，v_k～N(0～R_k)

其中，Z_k为第k时刻的实测自然电场数据，F为污染监测地的自然电场正演模型，v_k为观测误差或噪声，v_k服从协方差为R_k的零均值正态分布。

S13、初始化自然电场正演模型，获取初始化先验估计污染物扩散演化模型数据M₀。

可以理解的是，初始化即为选取某一时刻作为起始点，在合理范围内随机给定模型参数，获取该时刻的先验估计污染物扩散演化模型数据M₀。

S14、根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据。

可以理解的是，在获取M₀的基础上，依据状态模型，可以获取第k-1时刻的先验估计污染物扩散演化模型数据。该过程称为预测，如图5b所示，是利用前一时刻的状态估计来对当前时刻的状态进行预测估计，称之为先验状态估计。

S15、获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。

可以理解的是，在获取实测自然电场数据Z_k的基础上，依据状态模型和观测模型，可以对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。该过程称为更新，如图5b所示，是利用当前时刻的观测数据对先验状态估计进行校正优化，优化后的估计称之为后验状态估计。

本实施例将污染物扩散演化规律与自然电场实测数据结合，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。

图2示出了本发明另一实施例提供的一种自然电场监测数据动态反演方法的流程示意图，如图2所示，一种自然电场监测数据动态反演方法，所述方法包括：

S21、根据孔隙介质流体运动规律，建立污染监测地的状态模型。

S22、根据自然电场正演模拟方法，建立污染监测地的观测模型。

S23、初始化自然电场正演模型，获取初始化先验估计污染物扩散演化模型数据M₀。

可以理解的是，步骤S21-S23与图1所示的实施例步骤相同，在此不做赘述。

S24、根据卡尔曼滤波递归估计模型，获取第k-1时刻的先验估计污染物扩散演化模型数据和先验估计协方差P_k-1|k-1。

可以理解的是，如图6所示，卡尔曼滤波估计只需要系统前一时刻的状态估计和当前状态的观测数据来计算动态过程当前的状态估计。卡尔曼滤波估计中每次循环包含“预测”和“更新”两个步骤。

在本实施例中，基于卡尔曼滤波的时序数据反演，避免常规反演方法中每次观测数据反演均需要多次迭代进行正演计算，大大减少监测数据处理的计算代价。

S25、根据卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1。

可以理解的是，如图6所示，预测是利用前一时刻的状态估计来对当前时刻的状态进行预测估计，称之为先验状态估计。一般情况下，预测包括先验状态估计预测和先验估计协方差预测。

在本实施例中，所述卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=H_k{\hat{M}}_{k-1|k-1},$

$P_{k|k-1}=H_k{P}_{k-1|k-1}{H}_k^T+Q_k,$

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，Q_k为演化过程协方差矩阵。

S26、获取第k时刻的实测自然电场数据Z_k，根据卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k。

可以理解的是，如图6所示，更新是利用当前时刻的观测数据对先验状态估计进行校正优化，优化后的估计称之为后验状态估计。一般情况下，更新包括状态估计更新和估计协方差的更新。

具体来说，卡尔曼滤波估计的状态用和P_k|k两个变量来表示，前者为根据更新到k时刻的观测数据对k时刻系统状态的后验状态估计，后者为后验估计误差协方差矩阵，用来表示后验状态估计的精度。

在本实施例中，所述卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k{\tilde{y}}_k,$

P_k|k＝(I-K_kF_k)P_k|k-1，

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，为实测自然电场数据的数据残差，K_k为最优卡尔曼增益，

${\tilde{y}}_k=Z_k-F_k{\hat{M}}_{k|k-1},$

$K_k=P_{k|k-1}{F}_k^T{S}_k^{-1},$

其中，S_k为实测自然电场数据的数据残差协方差，

$S_k=F_k{P}_{k|k-1}{F}_k^T+R_k,$

其中，F_k为第k时刻的自然电场正演模型F，R_k为观测过程协方差矩阵。

在本实施例中，将污染监测地的自然电场实测数据不断输入上述卡尔曼滤波递归估计模型中，以不断校正优化模型估计，在模式解与实际观测之间找到一个最优解，同时继续为动态过程提供初始条件，不断循环迭代地输出优化模型估计。最终获得一系列监测目标区域的地下介质模型估计，实现对时序自然电场数据的动态反演。

本实施例将污染物扩散演化规律与自然电场实测数据结合，利用卡尔曼滤波递归估计进行预测和更新，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。基于卡尔曼滤波的时序数据反演，避免常规反演方法中每次观测数据反演均需要多次迭代进行正演计算，大大减少监测数据处理的计算代价。

图3示出了本发明另一实施例提供的一种自然电场监测数据动态反演方法的流程示意图，如图3所示，一种自然电场监测数据动态反演方法，所述方法包括：

S31、根据孔隙介质流体运动规律，建立污染监测地的状态模型。

S32、根据自然电场正演模拟方法，建立污染监测地的观测模型。

S33、初始化自然电场正演模型，获取初始化先验估计污染物扩散演化模型数据M₀。

可以理解的是，步骤S31-S33与图1所示的实施例步骤相同，在此不做赘述。

S34、设置预设加权点，对污染物扩散演化模型和自然电场正演模型进行无迹变换。

可以理解的是，由于污染物扩散演化和自然电场正演都具有明显的多维非线性特征，也就是说，污染物扩散演化模型H和自然电场正演模型F都是多维非线性函数。传统的非线性滤波方法需要将非线性系统线性化，如扩展卡尔曼滤波采用一阶线性化截断，忽略其余高阶项，带来了精度不高、稳定性差等缺点。

在本实施例中，如图7所示，无迹变换通过设置加权点(Sigma点)来近似表示多维目标采样点，对非线性函数的概率密度分布进行近似，用一系列确定样本来逼近状态的后验概率密度，无须对非线性函数直接线性化。无迹变换模拟非线性函数的概率分布而非模拟任意的非线性函数本身，比较容易实现。函数的高斯分布由一系列确定的加权点(Sigma点)来替代，可以避免计算雅克比矩阵，大大减少计算工作量。因此采用无迹变换来实现动态系统过程状态估计及其协方差矩阵的高效计算与非线性传递。

S35、根据基于无迹变换的卡尔曼滤波递归估计模型，获取每个加权点第k-1时刻的先验估计污染物扩散演化模型数据

S36、根据基于无迹变换的卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1。

在本实施例中，所述基于无迹变换的卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i{\hat{M}}_{i,k-1|k-1},$

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}}^T,$

其中，i为加权点的编号，2n为加权点的个数，为第i个加权点第k时刻的先验估计污染物扩散演化模型数据，W_i为每个加权点的权重系数。

S37、获取第k时刻的实测自然电场数据Z_k，根据基于无迹变换的卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k。

在本实施例中，所述基于无迹变换的卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k(Z_k-{\hat{M}}_{k|k-1}),$

$P_{k|k}=P_{k|k-1}-K_k{P}_{z_k{z}_k}{K_k}^T,$

其中，K_k为最优卡尔曼增益，为过程参数，

$K_k=P_{M_k{Z}_k}{P}_{Z_k{Z}_k}^{-1},$

$P_{Z_k{Z}_k}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{Z_{i,k}-{\hat{z}}_{k|k-1}\}{\{Z_{i,k-1|k-1}-{\hat{z}}_{k|k-1}\}}^T,$

$P_{M_k{Z}_k}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{Z_{i,k}-{\hat{z}}_{k|k-1}\}}^T,$

${\hat{z}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i{Z}_{i,k},$

其中，为过程参数，Z_i,k为第i个加权点第k时刻的实测自然电场数据。

S38、重复步骤S35-S37，获得各监测时刻的后验估计污染物扩散演化模型数据，实现对污染监测地的污染物扩散情况进行监测。

可以理解的是，每一个实测自然电场数据对应一个监测时刻。在本实施例中，将污染监测地的自然电场实测数据不断输入上述基于无迹变换的卡尔曼滤波递归估计模型中，以不断校正优化模型估计，在模式解与实际观测之间找到一个最优解，同时继续为动态过程提供初始条件，不断循环迭代地输出优化模型估计。最终获得一系列监测目标区域的地下介质模型估计，实现对时序自然电场数据的动态反演，如图8所示的动态反演结果示意图，左侧为模拟真实模型，右侧为反演结果模型。

本实施例将污染物扩散演化规律与自然电场实测数据结合，利用基于无迹变换的卡尔曼滤波递归估计进行预测和更新，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。采用无迹变换对多维非线性函数进行线性化处理，实现动态过程状态估计及其协方差矩阵的高效计算与非线性传递，可以避免计算雅克比矩阵，进一步减少计算工作量，使实时在线处理监测数据成为可能。

图4示出了本发明一实施例提供的一种自然电场监测数据动态反演系统的结构示意图，如图4所示，一种自然电场监测数据动态反演系统40，所述系统包括：

模型单元41，用于根据孔隙介质流体运动规律，建立污染监测地的状态模型，并根据自然电场正演模拟方法，建立污染监测地的观测模型；

反演单元42，用于根据初始化先验估计污染物扩散演化模型数据M0，通过所述模型单元41建立的状态模型和观测模型，获取第k时刻的先验估计污染物扩散演化模型数据，并根据第k时刻的实测自然电场数据，对所述第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据；

所述状态模型表示第k时刻和第k-1时刻先验估计污染物扩散演化模型数据的关系，所述观测模型表示第k时刻自然电场实测数据与第k时刻先验估计污染物扩散演化模型数据的关系。

优选地，所述反演单元42包括：

预测模块421，用于根据第k-1时刻的先验估计污染物扩散演化模型数据，通过所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；

更新模块422，用于根据第k时刻的实测自然电场数据，通过所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。

本实施例将污染物扩散演化规律与自然电场实测数据结合，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。

综上所述，本发明提供了一种自然电场监测数据动态反演方法及系统，通过建立污染监测地的状态模型和观测模型，利用卡尔曼滤波递归估计模型或基于无迹变换的卡尔曼滤波递归估计模型，实现数据监测的预测过程和更新过程。其中，预测过程是通过第k-1时刻的先验估计污染物扩散演化模型数据，根据状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；更新过程是通过第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。本发明将污染物扩散演化规律与自然电场实测数据结合，充分利用时序自然电场实测数据中的先后关联信息，增加自然电场实测数据反演中的先验知识与算法稳定性，为监测地下污染扩散情况提供动态反演结果。

本领域普通技术人员可以理解：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (10)

1.一种自然电场监测数据动态反演方法，其特征在于，所述方法包括：

根据孔隙介质流体运动规律，建立污染监测地的状态模型，所述状态模型表示第k时刻和第k-1时刻先验估计污染物扩散演化模型数据的关系；

根据自然电场正演模拟方法，建立污染监测地的观测模型，所述观测模型表示第k时刻自然电场实测数据与第k时刻先验估计污染物扩散演化模型数据的关系；

初始化自然电场正演模型，获取初始化先验估计污染物扩散演化模型数据M₀，并根据所述状态模型，获取第k-1时刻的先验估计污染物扩散演化模型数据；

根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；

获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。

2.如权利要求1所述的方法，其特征在于，

在建立状态模型之前，所述方法还包括：

根据孔隙介质流体运动规律，通过求解理查德方程，建立污染监测地的污染物扩散演化模型；

在建立观测模型之前，所述方法还包括：

根据自然电场正演模拟方法，通过求解泊松方程，建立污染监测地的自然电场正演模型。

3.如权利要求2所述的方法，其特征在于，根据孔隙介质流体运动规律，建立污染监测地动态系统的状态模型，所述状态模型公式如下：

M_k＝H(M_k-1)+w_k，

其中，M_k为第k时刻的先验估计污染物扩散演化模型数据，M_k-1为第k-1时刻的先验估计污染物扩散演化模型数据，H为污染监测地的污染物扩散演化模型，w_k为演化误差或噪声。

4.如权利要求3所述的方法，其特征在于，根据自然电场正演模拟方法，建立污染监测地的观测模型，所述观测模型公式如下：

Z_k＝F(M_k)+v_k，

其中，Z_k为第k时刻的实测自然电场数据，F为污染监测地的自然电场正演模型，v_k为观测误差或噪声。

5.如权利要求4所述的方法，其特征在于，所述根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据，包括：

获取第k-1时刻的先验估计污染物扩散演化模型数据和先验估计协方差P_k-1|k-1；

根据卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1，所述卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=H_k{\hat{M}}_{k-1|k-1},$

$P_{k|k-1}=H_k{P}_{k-1|k-1}{H}_k^T+Q_k,$

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，Q_k为演化过程协方差矩阵。

6.如权利要求5所述的方法，其特征在于，所述获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据，包括：

获取第k时刻的实测自然电场数据Z_k；

根据卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k，所述卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k{\tilde{y}}_k,$

P_k|k＝(I-K_kF_k)P_k|k-1，

其中，H_k为第k时刻的污染监测地的污染物扩散演化模型H，为实测自然电场数据的数据残差，K_k为最优卡尔曼增益，

${\tilde{y}}_k=Z_k-F_k{\hat{M}}_{k|k-1},$

$K_k=P_{k|k-1}{F}_k^T{S}_k^{-1},$

其中，S_k为实测自然电场数据的数据残差协方差，

$S_k=F_k{P}_{k|k-1}{F}_k^T+R_k,$

其中，F_k为第k时刻的自然电场正演模型F，R_k为观测过程协方差矩阵。

7.如权利要求4所述的方法，其特征在于，所述根据所述第k-1时刻的先验估计污染物扩散演化模型数据以及所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据，还包括：

设置预设加权点，对污染物扩散演化模型和自然电场正演模型进行无迹变换；

获取每个加权点第k-1时刻的先验估计污染物扩散演化模型数据

根据基于无迹变换的卡尔曼滤波递归估计模型的预测公式，获取第k时刻的先验估计污染物扩散演化模型数据以及先验估计协方差P_k|k-1，所述基于无迹变换的卡尔曼滤波递归估计模型的预测公式如下所述：

${\hat{M}}_{k|k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i{\hat{M}}_{i,k-1|k-1},$

$P_{k|k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}}^T,$

其中，i为加权点的编号，2n为加权点的个数，为第i个加权点第k时刻的先验估计污染物扩散演化模型数据，W_i为每个加权点的权重系数。

8.如权利要求7所述的方法，其特征在于，所述获取第k时刻的实测自然电场数据，并根据所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据，包括：

获取第k时刻的实测自然电场数据Z_k；

根据基于无迹变换的卡尔曼滤波递归估计模型的更新公式，获取第k时刻的后验估计污染物扩散演化模型数据以及后验估计协方差P_k|k，所述基于无迹变换的卡尔曼滤波递归估计模型的更新公式如下所述：

${\hat{M}}_{k|k}={\hat{M}}_{k|k-1}+K_k(Z_k-{\hat{M}}_{k|k-1}),$

$P_{k|k}=P_{k|k-1}-K_k{P}_{z_k{z}_k}{K_k}^T,$

其中，K_k为最优卡尔曼增益，为过程参数，

$K_k=P_{M_k{Z}_k}{P}_{Z_k{Z}_k}^{-1},$

$P_{Z_k{Z}_k}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i\{Z_{i,k}-{\hat{z}}_{k|k-1}\}{\{Z_{i,k-1|k-1}-{\hat{z}}_{k|k-1}\}}^T,$

$P_{M_k{Z}_k}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i\{{\hat{M}}_{i,k-1|k-1}-{\hat{M}}_{k|k-1}\}{\{Z_{i,k}-{\hat{z}}_{k|k-1}\}}^T,$

${\hat{z}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{W}_i{Z}_{i,k},$

其中，为过程参数，Z_i,k为第i个加权点第k时刻的实测自然电场数据。

9.一种自然电场监测数据动态反演系统，其特征在于，所述系统包括：

模型单元，用于根据孔隙介质流体运动规律，建立污染监测地的状态模型，并根据自然电场正演模拟方法，建立污染监测地的观测模型；

反演单元，用于根据初始化先验估计污染物扩散演化模型数据M0，通过所述模型单元建立的状态模型和观测模型，获取第k时刻的先验估计污染物扩散演化模型数据，并根据第k时刻的实测自然电场数据，对所述第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据；

所述状态模型表示第k时刻和第k-1时刻先验估计污染物扩散演化模型数据的关系，所述观测模型表示第k时刻自然电场实测数据与第k时刻先验估计污染物扩散演化模型数据的关系。

10.如权利要求9所述的系统，其特征在于，所述反演单元包括：

预测模块，用于根据第k-1时刻的先验估计污染物扩散演化模型数据，通过所述状态模型，获取第k时刻的先验估计污染物扩散演化模型数据；

更新模块，用于根据第k时刻的实测自然电场数据，通过所述状态模型和观测模型，对第k时刻的先验估计污染物扩散演化模型数据进行修正，获取第k时刻的后验估计污染物扩散演化模型数据。