随机模拟原理

随机模拟是以随机函数理论为基础的。随机函数由一个区域化变量的分布函数和协方差函数来表征。一个随机函数Z(x)有无数个可能的实现，模拟的基本思想是从一个随机函数抽取多个可能的实现。在模拟中可进行条件限制，即用观测点的采样支队模拟过程进行条件限制，使得采样点的模拟值与实测值相同，这就是我们常说的条件模拟。

随机模拟与克里格等插值方法有较大差别，主要表现在以下三个方面。

(1)插值法考虑局部估计值的精确程度，力图对估计点的未知值做出最优(估计方差最小)和无偏(估计值的均值与观测点的均值相同)的估计，不考虑估计值的空间相关性，而模拟首先考虑的是结果的整体性质和模拟值的统计空间相关性，其次才是局部估计值的精度。

(2)插值法给出观测值的平滑估计，而削弱了观测数据的离散性，忽略了数据之间的细微变化;而条件随机模拟通过在插值模型中系统的加上了“随机噪音”，使得产生的结果比插值模型真实得多。

(3)插值法只插身一个实现，在随机模拟中则产生许多可选的模型，各模型之间的差别正是空间不确定性的反映。

(一)随机模拟的主要方法

随机模拟的基本思想非常简单，但依据不同的研究对象、不同的样品信息来源和不同的地质复杂程度，随机模拟的方法多种多样。概括起来，主要有以下五种基本类型。

1．序贯模拟

序贯模拟是一种思路，高斯型、指示型模拟都可以采用，即是将u邻域内的所有已知数据(包括原始数据和先前已模拟实现的数据)作为后续模拟的条件数据，在此基础上进行下一步的模拟。考虑N个随机变量Zi的联合分布，这里N可以很大，Zi可以代表一个区域内离散在N个网格点的同一属性，也可以是同一点的N个不同属性。模拟时的抽样是在变量的ccdf中进行的。

使用序贯模拟方法，需要确定越来越复杂的累积条件概率分布函数(ccdf)。在实际计算过程中，由于数据的屏蔽作用，只有最近的数据作为条件，然而即使比较远的原始数据也应该加以重视。

2．高斯模拟

高斯随机函数:F(u1，…，uk;Z1，…，Zk) =Prob{Z(uk≤Zk)}的任何k元cdf完全由来自协方差函数:Cov(u，u') =E{Z(u)Z(u')} －E{Z(u)Z(u')}E{Z(u')}的知识确定。这样，高斯模型就显得容易实现，因而得到了广泛应用。但要注意在使用过程中要先对原始数据进行正态转换。

高斯模拟是将地质变量作为符合高斯分布的随机变量，空间上作为一个高斯随机场，以高斯随机函数来描述，即随机函数是多元正态的，当且仅当:

(1)所有Y(¯u) 的子集，如也是多元正态的。

(2)随机函数Y(¯u)的成分的所有线性组合也是正态分布的。

(3)给定任何其他子集的实现，随机函数Y(¯u)的任何子集的条件分布也是多元正态分布。

当储层参数具有连续的特性时，我们可以将其分解为一个连续的或近似连续的随机曲面外加一个具有空间相关性的高斯随机噪声的模型。

设Z(* )是一个高斯过程，且为一个回归趋势面m(x，β)外加一个零均值的高斯噪声［ε(x)］，即:

这里β(β1，…，βp)T是一个向量;fj是回归函数，常取f1(x)=1;βp是回归参数。

协方差函数为:

Cov{z(x)，z(y)} =Cov{ε(x)，ε(y)} =σ2k(x，y)　(6-39)

可以用来度量噪声的空间连续性。如果要对x处进行估计，从高斯假设条件可知:

这里K是一个n×n的以K(xi，xj)为元素的矩阵，称为协方差矩阵，此外还有:

k=(K(x，x1)，…，K(x，xn))T

Zdata=(Z(x1)，…，Z(xn))T　(6-41)

而F是一个n×p的矩阵，第i行为:

f(xi)T=(f1(xi)，f2(xi)，…，fp(xi))　(6-42)

从而在观测数据给定的条件下，随机变量Z(x)的分布可以写成如下的高斯分布:

Z(x)|data∝N{m(x，β) +k TK－1(Zdata－mc(x，β))，σ2(K(x，x) －k TK－1k)} (6-43)这里，mc(x，β)是一个向量，其分量为m(xi，β)。由于mc(x，β) =Fβ，所以一般来讲，该向量要随着回归函数fi(x)的改变而变化。

利用最小二乘法对β进行估计，使得(Zdata－Fβ)TK－1(Zdata－Fβ)达到最小。当协方差矩阵K已经确定时，β的最小二乘估计βLS可写成如下形式:

βLS=HFTK－1Zdata　(6-44)

最后，在x处的内插结果可写成:

Z(x) =m(x，βLS) +k TK－1(Zdata－mc(x，βLS))

=f T(x)βLS+k TK－1(Zdata－FβLS)　(6-45)

顺序高斯模拟算法包括下列步骤。

(1)将所有条件数据(硬数据和已模拟实现的数据)进行正态变换，从非正态分布变换为正态分布，作为先验条件概率分布。

(2)对变换后的条件数据进行变差函数计算与建模。

(3)对模拟网格，定义一个模拟的随机路径。

(4)对每个网格点，根据邻近的条件数据(包括井的硬数据和已模拟的数据)和地震波阻抗数据进行协克里格估计，得到后验条件概率分布。

(5)从后验条件概率中进行随机抽样，得到该网格点上的一个模拟值。

(6)将模拟值加到条件数据中去。

(7)重复(4) ～(7)，直至所有网格点模拟一遍。

(8)进行高斯条件模拟值的正态逆变换。

3．指示模拟

指示随机函数模型适用于k元点统计量所控制的离散变量的模拟，其最重要的应用在于它可以给出条件概率最直接的评价，这正是序贯模拟时所必需的。它对离散变量和连续变量都可以进行模拟。

我们在直角坐标系下，将地质模型离散为n+m个网格单元，并以整型数i=1，…，n+m进行标记，对每个网格i，我们以砂/泥两种岩相为例，定义一个二值指示变量:

假设xi为满足二阶平稳条件的随机场，并具有稳定的均值:

E{xi} =p(xi=1} =πsand　(6-47)

该均值与模型的砂岩比例相关，其空间上的协方差为:

E{［xi－πsand］［xj－πsand］} =C(hij)　(6-48)

这里，hij是网格i与j之间的内部距离矢量。一般来讲，协方差与i和j之间的距离、矢量hij的方向有关，这就使我们得以模拟岩相在空间上的各向异性。

在实际应用当中，对于砂/泥模型矢量X={x1，…，xn+m}，设井的数据作为条件数据x1，…，xn，则待估数据为xn+1，…，xn+m。

对于砂/泥两种岩相，分解式(6-50)可用图6-19来表示。结合式(6-49)、式(6-50)，顺序指示模拟过程包括以下6个基本步骤。

(1)随机地提取一个未模拟的网格i。

(2)用指示克里格(IK)估计砂/泥局部先验概率分布;

对于砂岩的克里格方程组为:

(3)结合(2)和(3)获得砂/泥的局部后验概率分布。

(4)对局部后验概率分布进行随机抽样，得到网格i的一个模拟值xi。

(5)将(4)所得到的模拟值xi作为模拟xi+1时的条件数据:

(6)重复(1) ～(5)，直至所有网格都被模拟一次。

上述步骤可用图6-19来说明，图中X1、X2、X3为井数据，X10为第10个模拟网格处的地震属性，X4～X9为已经模拟完成的实现值。这样，在模拟到第10个网格时，以X1～X9为条件数据，通过指示克里格估计得到由条件数据产生第10个网格处的局部概率分布(式6-49);最后，通过对该分布的随机抽样来得到该网格处的一个随机实现。图中的箭头表示模拟的随机路径，由于网格一经模拟之后，就将被加入到条件数据集，所以，避免了抽样不均匀的问题，且每个网格仅被模拟一次，在整个模拟过程中，条件数据集是在不断累加的。

4．布尔模拟

布尔模拟是用于离散变量的随机模拟算法，示性点过程模拟是布尔模拟的扩展。在地学中主要用于推断，如相模拟、矿体分布等。

5．退火模拟

退火模拟是一种条件随机模拟方法，可以综合再现两点统计量和复杂多元空间统计量，连续型和离散型变量都适用。开始时产生一初始图像，同时将条件数据置于最近结点，在其余结点随机地从直方图中抽取样本，然后序贯地交换一对非条件结点数据，计算目标函数(实验与模型变差函数的平均平方差)，目标函数减小时接受交换，直到进一步的交换不能使目标函数降低或达到最小值为止。

(二)随机模拟的方法分类

随机模拟方法有许多种，各种随机模拟方法在其基本原理、复杂程度和应用条件诸方面均有所不同。不同学者从不同角度对随机模拟方法进行了分类。分类主要以数据的分布类型、变量的类型、模拟的结果是否忠实原始数据、模拟实现的过程，以及模拟所设计的变量数来加以区分，如表6-1所示。

图6-19 模拟实现示意图

按照变量的特征，将随机模拟分为离散性和连续性随机模拟。离散性随机模拟方法主要用于模拟离散型变量，如岩性、沉积相等;连续性随机模拟方法主要用于模拟连续型变量，如孔隙度、渗透率等。

从数据分布类型角度看，各种高斯模拟算法(分形随机函数法、矩阵方法、频率域方法、序贯高斯模拟、基于各种克里格估计的高斯模拟等)适用于服从正态(高斯)分布的数据;利用高斯得分转换可将非高斯分布的数据转换为高斯分布的数据。其他模拟算法(如各种指示模拟方法、模拟退火方法、转向带法、布尔模拟方法等)则不要求数据为高斯分布。

从变量类型角度看，各种基于指示克里格的模拟(序贯高斯模拟、基于Markov-Bayes方法的模拟、概率场模拟)、截断高斯模拟、布尔模拟适用于类型变量的模拟;序贯高斯模拟、序贯高斯协同模拟、分形随机函数法、矩阵方法、频率域方法适用于连续型变量的模拟;模拟退火模拟(包括Metropolis算法、热浴法、协同模拟退火法、并行模拟退火法)既适用于类型变量，也适用于连续型变量。

从模拟实现的过程看，可以分为基于目标的模拟和基于像元的模拟。前者包括标点过程法、Markov随机域法、截断高斯模拟、两点直方图法、随机成因模拟法和指示模拟法，用来表现各种地质特征的空间分布;后者包括模拟退火法、序贯指示模拟、分形随机域法、Markov随机域法、LU分解法、转向带法，用来模拟各种连续性参数(岩石物性参数)，以及离散性参数(地质特征参数)的变化。

从参与模拟的变量数目看，有单变量模拟和多变量模拟。前者还包括结合软硬数据的单变量模拟(基于同位协同克里格的高斯模拟、基于全协同克里格的高斯模拟、序贯高斯协同模拟和基于Markov-Bayes模型的模拟)和结合二级变量的单变量模拟(基于泛克里格的模拟、基于外部漂移泛克里格的模拟和基于局部变化均值的简单克里格的模拟)，后者包括K个相关变量的联合条件模拟、基于Markov型协区域化模型的多变量联合模拟、相关变量组成的矢量模拟、化为独立变量的方法和从条件分布中抽取的方法。

表6-1 模拟方法特征