一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法与流程

本发明属于水文统计学技术领域，具体涉及一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法。

背景技术：

地下水状态跟人类社会的可持续发展息息相关。一个值得信赖的地下水水流和溶质运移模拟在地下水水位预报，溶质运移预测，地下水资源管理以及地下水环境危险性评估方面是具有非常重要的意义。地下水水流和溶质运移模拟的质量好坏在很大程度上取决于地下水水文参数刻画的好坏(例如，水文渗透系数，孔隙度，等等)。然而，相对于含水层的尺度而言，在实际野外得到的这些水文参数的试验测量数据是非常稀少地。而稀少的水文地质参数测量数据是很难被用来直接刻画研究区域水文地质参数场。所以，如何利用这些稀少的野外水文地质参数来尽可能准确的刻画整个含水层的水文地质参数场是一个非常大的挑战。特别是对于大尺度的研究区而言，在最近几十年来，同化实时状态数据的随机逆模拟已经被证明是一个非常有用的能够刻画水文地质参数场的工具。

并且，综合以上研究还可发现，当今主流的反演非线性高斯参数场的逆模拟方法主要还是集中在集合卡尔曼滤波、集合平滑法、逆序贯模拟法或基于马尔科夫蒙特卡洛的贝叶斯方法当中。其中，尽管基于mcmc的方法计算相当耗时并且采样效率低下，但它是经典的可以处理复杂参数分布及非线性问题的逆模拟方法。反观其他方法，他们都依赖于某种线性化或一阶、二阶矩近似，虽然这些方法通常比马尔科夫蒙特卡洛计算快得多，但随着非线性问题的增加，反演的结果会变得越来越不精确并且会造成不稳定的近似。因此，对于需要处理高精度的问题，以及要严格检验近似估算方法的先进性的时候，渐近精确方法(如马尔科夫蒙特卡洛)就是唯一的解决方案。但是马尔科夫蒙特卡洛方法计算非常耗时，难以解决大尺度高维度问题。

技术实现要素：

发明目的：为了克服背景技术的不足，本发明公开了一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，该方法通过耦合并行退火技术，通过交换冷链和热链的采样信息从而使得目标链能够探索到跟大的模型空间使得马尔科夫目标链能够探索到更大的模型空间，从而提高探索目标后验分布来提高具有高斯分布的水文地质参数场的反演效率。

技术方案：本发明的一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，包括以下步骤：

s1、建立初始高斯参数场、设定平行退火技术中需要的温度梯度以及预条件克兰克尼科尔森建议分布计算中的跳跃梯度；

s2、通过预条件克兰克尼科尔森建议分布得到当前的建议样本，用于计算建议样本的似然度；

s3、根据建议分布和当前分布的似然度计算所有马尔科夫链的当前的接受概率，再根据metroplishastings法则，判断新的样本是建议样本，还是当前样本；

s4、根据设定的交换频率，通过并行退火技术计算任意一对不同温度的马尔科夫链的交换接受概率，再根据metroplishastings法则，判断冷链、热链是否交换，使得目标链能够尽快的探索整个模型空间；

s5、重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样。

其中，s1具体为：

设立初始高斯场服从均值为0，方差为c的高斯分布，即

设立一个温度梯度，t1<t2<…<ti<…tn，t1＝1，设立一个跳跃因子梯度，β1<β2<…<βi<…βn，βn<1，

式中n且为总的马尔科夫链的个数，i表示马尔科夫链的编号，i∈(1,n)，并且交换建议频率为d。

进一步的，s2具体为：

在第k^th个样本迭代生成一个预条件克兰克尼科尔森建议分布且对于所有马尔科夫链i＝1,...,n，

进一步的，s3具体为：

对于每个马尔科夫链i，当满足接受概率时，则下一个样本否者下一个样本维持不变，即

式中为建议分布的对数似然度，为当前分布的对数似然度。

进一步的，s4具体为：

对于任意一对不同温度的马尔科夫链，第i^th和第j^th条马尔科夫链，如果且当满足交换接受概率时，这两条链进行交换否者维持不变，

有益效果：与现有技术相比，本发明的优点为：首先，本发明能够刻画水文地质参数高斯场，能够处理水文地质参数反演中的高维度线性及非线性高斯场问题；其次，本方法是一个渐进精确的采样方法，在处理地质统计学逆模拟问题中，比预条件克兰克尼科尔森马尔科夫蒙特卡洛方法效率更高。

附图说明

图1为实施例中观测井和抽水井分布示意图；

图2为本发明中贝叶斯逆模拟方法示意图；

图3为实施例中对数水力渗透系数参考场的频率直方图；

图4为实施例中pcn-mcmc和pcn-pt采样的对数水力渗透系数频率直方图；

图5为实施例中pcn-mcmc和pcn-pt采样的均方根误差；

图6为实施例中pcn-mcmc和pcn-pt潜在尺度因子r值的变化曲线。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

为了验证本发明是可以加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，在一个饱和，稳定流的承压含水层中验证，该含水层空间二维为5000m*5000m，厚度为50m，且离散为100*100*1个单元格，每个单元格为50m*50m*50m，东、西边界为隔水边界且给定水位分别为20m和0m，南、北边界为隔水边界，并如图1所示，给定25口观测井和4口抽水井，圆圈代表观察井，星号代表抽水井，且抽水井#1、#2、#3、#4的抽水120m³/d，70m³/d，90m³/d，90m³/d。通过序贯高斯模拟[deutschandjournel,1998]，并依据均值为2.5ln[m/d]，标准差为2ln[m/d]，最大相关长度和最小相关长度分别为2000m和1500m，角度为135度的指数变异函数模型生成高维度的服从高斯分布的对数水文渗透系数参考场。

通过地下水流动模拟器modflow[mcdonaldandharbaugh,1988]计算如下地下水稳定流流动方程式：

式中：表示散度算子；k为水力渗透系数[m/d]；为拉卜拉算子；h为水位值[m]；q为含水层单位体积的体积流量[1/d]。

pcn-pt：耦合并行退火技术的预条件克兰克尼科尔森马尔科夫蒙特卡洛方法

pcn-mcmc：传统的预条件克兰克尼科尔森马尔科夫蒙特卡洛方法

如图2所示的本发明的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法(pcn-pt)，包括步骤如下：

s1、建立初始高斯参数场、设定平行退火技术中需要的温度梯度以及预条件克兰克尼科尔森建议分布计算中的跳跃梯度；

s2、通过预条件克兰克尼科尔森建议分布得到当前的建议样本，用于计算建议样本的似然度；

s3、根据建议分布和当前分布的似然度计算所有马尔科夫链的当前的接受概率，再根据metroplishastings法则，判断新的样本是建议样本，还是当前样本；

s4、根据设定的交换频率，通过并行退火技术计算任意一对不同温度的马尔科夫链的交换接受概率，再根据metroplishastings法则，判断冷链、热链是否交换，使得目标链能够尽快的探索整个模型空间；

s5、重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样。

具体实施时：

(1)设立初始高斯场服从均值为0，方差为c的高斯分布，即设立一个温度梯度且t1<t2<…<ti<…tn，且t1＝1；设立一个跳跃因子梯度β1<β2<…<βi<…βn，且βn<1，式中n且为总的马尔科夫链的个数，i表示马尔科夫链的编号，i∈(1,n)，并且交换建议频率为d。

(2)在第k^th个样本迭代生成一个预条件克兰克尼科尔森建议分布且对于所有马尔科夫链i＝1,...,n，

(3)对于每个马尔科夫链i，当满足接受概率时，则下一个样本否者下一个样本维持不变，即

式中为建议分布的对数似然度；为当前分布的对数似然度。

(4)对于任意一对不同温度的马尔科夫链，例如，第i^th和第j^th条马尔科夫链，如果且当满足交换接受概率时，这两条链进行交换否者维持不变

(5)重复以上过程进行下一个高斯场样本采样。

对于本发明的方法pcn-pt，设定20条并行的马尔科夫链，对于温度梯度为1.76^i-1，跳跃因子成指数衰减且为0.78^20-i+。

pcn-mcmc的跳跃因子为0.005。采样的样本总数为800000个。

如图3所示，显示的是对数水力渗透系数参考场的频率直方图，如图4所示，显示的是通过pcn-mcmc和pcn-pt采样的对数水力渗透系数场的频率直方图(已扣除马尔科夫链的burn-in阶段)，可以看出这两种方法采样之后都能获得参考场的高斯分布特征且统计信息基本一致。如图5所示，显示的是pcn-mcmc和pcn-pt样本的均方根误差，可以看出通过pcn-pt方法计算的样本均分根误差能够更快的减小并且收敛到均方根误差均值。如图6所示，显示的是两种方法在随着马尔科夫链增长的的潜在尺度因子(potentialscalefactor)r变化曲线(当r值小于1.2是表明采样样本达到稳定分布)，可以看出，虽然在这两种方法中，r值都随着马尔科夫链的增长而减小，但是pcn-pt方法采样计算的r值在采样结束前大部分都能小于1.2，而pcn-mcmc计算的r值在采样结束之前，大部分都还大于1.2，由此可知，发明的新方法pcn-pt能够更快的达到一个稳定分布。

通过以上分析发明的新方法pcn-pt能够更快速高效的处理高维度的非线性高斯场问题。