1. 引言
储层随机建模技术在现代油藏描述中得到越来越广泛的应用 [1] [2]，该项技术能够更好地刻画储层的非均质性和定量描述储层中的不确定性 [3] [4] [5]。储层随机建模方法按照模拟对象的不同，可以分为基于像元(Pixel)、目标(Object)和样式(Pattern)三类 [6] [7] [8]。基于像元的建模方法在模拟的过程中是一个网格一个网格逐个的模拟，模拟的顺序是随机的，因此建立的模型通常会含有一部分“噪点”，也就是分布很离散、缺乏地质意义的一些相对孤立的网格点或少量网格点组合。针对这一问题，加拿大学者Clayton提出了MAPS平滑方法 [9]，该方法通过一个加权平均窗口对初始模拟结果进行平滑，能够过滤掉认为不符合地质认识的孤立的网格点。国内学者尹艳树在此基础上提出了一种基于信息度的平滑方法 [10]，在权值确定中考虑了条件数据的影响，进一步完善了该平滑方法。在商业化软件中也增加了平滑处理的模块。在多数情况下，平滑方法能够取得不错的效果。当模型中含有薄的呈线性(二维)或片状(三维)分布的地质体时，例如薄的泥岩夹层，这些地质体对渗流又有重要的影响 [11] [12] [13] [14]，需要在模型中进行表征。MAPS方法会把这种地质体完全平滑掉，而不仅仅是把噪点去掉。针对这一问题，对平滑方法MAPS进行了改进，在存在薄的地质体时，针对性地设计平滑的窗口和权值，取得了较好的效果。同时对于边界上的噪点也进行了处理，原始的MAPS方法没有考虑这一问题。
2. 问题的提出
序贯指示模拟(SISIM)、多点地质统计学SNESIM等基于像元的建模算法，在模拟离散变量(如岩性、岩相等)过程中，采用的是根据随机路径确定的网格点逐点进行模拟。因为是单个网格单元的模拟，容易导致模拟结果出现相对离散的点，这些离散的点不符合地质认识，是所谓的“噪点”。Clayton教授提出了一种模型平滑方法MAPS可以在一定程度上解决该问题。在MAPS方法中，主要采用了一种平滑窗口，如二维情况可以采用3 × 3、5 × 5窗口对模型进行平滑处理。图1为平滑前后的对比，一般情况下MAPS可以取得较好的效果。但是在某些特殊情况下，现有的平滑方法不适用。如图2所示，该图为阿尔伯塔大学地质统计学授课中的平滑原理示意图，采用3 × 3的模板，权值分布采用的是各向同性的方法，即中心点权值最大，距离中心点距离相同的网格的权值一样。当平滑窗口处于如图中红色虚线框所示位置时，中心点所属相类型的概率采用加权平均的方法。例如平滑窗口共有9个网格，权值之和为15，而属于代码1的有两个网格，其权值分别为1和3，因此中心点属于代码1的概率为4/15 = 0.27。同理，属于代码2的概率为(2 + 1 + 2 + 2)/15 = 0.47，属于代码3的概率为(1 + 2 + 1)/15 = 0.27。因此，平滑后红色虚线框内的中心网格的属性将用概率最大的代码2所替代，这样就去除了离散的点。在研究中发现，当平滑窗口位于左下角时，如图中蓝色虚线框所示，则中心点属于代码2的概率为8/15 = 0.53，而属于代码3的概率为7/15 = 0.47，因此中心点会被代码3所替代。继续平移平滑窗口(向右平移一个网格)，下一个平滑窗口中心点的绿色(代码3)也将被白色(代码2)所替代。如此重复下去，整个绿色条带都将被白色所替代。通过图2可以直观的看出，绿色条带并不是孤立的，连续性很好，只是因为呈线性分布，所以这种平滑方法针对线性分布的地质体是不适用的。
>原始SISIM模拟结果 
1s 平滑后结果
Figure 1. Comparison before and after smoothing
图1. 平滑前后对比
Figure 2. Schematic diagram of smoothing
图2. 平滑处理的原理示意图
3. 方法的改进
改进处理方法主要针对MAPS方法中平滑窗口进行重新设计，针对特殊的地质情况设计合理的平滑窗口和权值。原始MAPS算法中平滑窗口为各向同性，而泥岩夹层的分布主要呈线性分布(二维)或薄片状分布(三维)，不同方向的连续性差别较大。采用各向同性的平滑窗口不仅把孤立的点平滑掉，而且也会把线性分布的泥岩平滑掉。以薄层泥岩的平滑处理为例说明，将各向同性平滑窗口改进为长方形的各向异性平滑窗口，如图3所示，窗口不再是规则的正方形网格，而是根据地质体的空间分布特征来设置，如本例中水平方向设置为5个网格，而在垂直方向上设置为3个网格，并且权值不再是各向同性分布，例如本例中中间网格的权值明显高于上下层网格的权值，这种窗口适合于图2中呈线性分布地质体的平滑。在实际应用中需要根据地质体的特点来设计合理的平滑窗口。
使用改进后的平滑窗口对图2进行平滑，当平滑窗口位于图中黄色虚线时。平滑窗口共有15个网格，权值之和为59。属于代码1概率为1/59 = 0.0169，属于代码2的概率为9/59 = 0.1525，属于代码3的概率为49/59 = 0.8305。平滑后图中绿色相带就保留下来了，而图中红色虚线内的离散点按照此平滑窗口也会正确的被白色所替代。提高平滑窗口中心点及其左右两边网格的权重可以满足了线性分布岩相类型的平滑处理。
下面采用通用模板进一步说明改进平滑窗口的工作原理以及需要注意的问题：假设岩相模型中单层泥岩分布代码如下，0代表背景或者其他岩相类型，1代表泥岩夹层，红色表示泥岩分布的两端，1的数量为L，表示泥岩长度。
Figure 4. Distribution of single mudstone layer
图4. 单层泥岩分布
使用图5所示3*m的平滑窗口去平滑图4所示泥岩夹层。当窗口中心点K移动到图4 中左边红色1处时，属于泥岩的权值为K + K1 + … + Kn，属于其他相的权值为K1 + … + Kn + sum(‘1’)，sum(‘1’)表示平滑窗口中所有权值为1的和。为保证泥岩夹层两端(图4中红色1)不会被其他相类型(图4代码0)替代，属于其他相的权值应该小于属于泥岩的权值，所以平滑窗口应该满足K > sum(‘1’)，这样平滑过程中泥岩夹层两端不会被平滑掉，泥岩夹层中间也不会被平滑掉，这样整个泥岩夹层就完整保留下来。图5所示平滑窗口长度为m，图4所示泥岩夹层长度L，为保证窗口中心点K移动到图4中左边红色1处时，窗口边缘点Kn在图4右边红色1的左边，需要满足n ≥ L，其中m = 2*n + 1，这样上述权值才能正确计算，泥岩夹层才能保留下来。所以若想在平滑过程中不保留长度为L以下的泥岩夹层，则让n < L，并且Kn > K-sum(‘1’)，这样在移动窗口的过程中所有泥岩夹层的权值都小于背景相权值而被平滑掉。
原始MAPS方法在平滑过程中当模板中心点K位于模型的左上角时，平滑窗口的左上部分在泥岩夹层模型外部，只有中心点K右面及下面平滑窗口权值参与计算。这样平滑窗口中的权值没有被完全使用，会损失部分信息。在改进方法中对原始泥岩夹层模型进行向外扩充，比如使用3 × 5平滑窗口对网格数量为200 × 200的泥岩夹层模型进行平滑，事先在模型上面和下面各增加一层背景相，在模型左面及右面各增加两层背景相，整个模型的网格数量为202 × 204。这样平滑窗口在202 × 204模型网格内部平滑就不会出现平滑窗口溢出到模型外部的情况，使得平滑窗口中每个权值都得到计算，避免了模型边界处的平滑误差。
4. 改进前后的对比
在老油田SISIM方法经常被用于泥岩夹层模型的建立，一般来说由于网格大小的限制以及泥岩夹层厚度分布的特点，泥岩夹层一般呈条带状分布。如图6所示例子，网格大小为10 × 1米，泥岩夹层的厚度为1~2米，平均长度为200米，泥岩比例为10%。采用SISIM模拟，变差函数设置为水平方向变程为200米，垂直方向2米，理论模型为球状模型。图6左边为模拟的结果，可以看出夹层呈水平条带分布，局部存在一些离散的点或是横向延伸很短的线。右边为泥岩夹层的长度分布直方图。在统计泥岩延伸长度的过程中，由于相邻单层网格的泥岩会出现在垂向上叠加而形成较长的泥岩夹层的情况，单独统计某一层网格上的泥岩长度与实际泥岩长度不符，所以本文设置统计泥岩长度规则如下：所有能够与某一泥岩网格在上下左右或对角上相邻的泥岩网格都属于同一个泥岩夹层，网格在统计过程中不会被重复计算。统计夹层长度算法先遍历整个模型，提取出所有独立的相互不连通的泥岩夹层，再计算每个泥岩夹层的最大长度。计算的原始模型中泥岩长度平均值为140.4米，比输入的期望值200米(水平方向变程)小了近1/3，主要原因是因为存在一些离散的噪点。采用MAPS方法平滑，平滑样板3 × 3，平滑结果如图7所示，泥岩长度平均值为92.6米，小于平滑前的泥岩长度平均值。对比图6和图7可以看出，采用MAPS平滑泥岩夹层模型有如下特点，对于单层的泥岩夹层，MAPS方法会将其完全平滑为背景相，对于离散噪点MAPS方法去除不完全，处理后不仅仅导致泥岩长度变短，而且很多泥岩夹层都被平滑掉了，也大大改变了原始的泥岩百分比。
针对MAPS方法平滑存在的不足，采用上节提出的改进方法，按照改进平滑窗口设置规则，如果需要平滑掉长度为30米及30米以下的泥岩夹层，设置3 × 9大小的模板[1 1 1 1 1 1 1 1 1; 4 5 6 7 20 7 6 5 4; 1 1 1 1 1 1 1 1 1] 对泥岩夹层模型(图6)进行平滑，迭代两次。平滑后的结果如图8所示。可以看出平滑后泥岩夹层最短长度为30米，所有30米及以下的泥岩夹层噪点都被平滑掉了。新平滑方法保留了原始泥岩夹层的分布特点，针对噪点及认为不合理的泥岩短线进行了剔除。统计直方图显示平均泥岩长度为164米，相对于MAPS处理后模型的140米，更接近期望值200米。
Figure 6. Interlayer distribution model simulated by SISIM
图6. SISIM模拟夹层分布模型
Figure 7. Smooth effect of MAPS method
图7. MAPS方法平滑效果
Figure 8. Smoothing results with mudstone length greater than 30 m
图8. 泥岩长度大于30米的平滑结果
Figure 9. Smoothing results with mudstone length greater than 40 m
图9. 泥岩长度大于40米的平滑结果
该算法设计灵活，可以根据实际资料的情况设置不同的下限值，例如设置3 × 11大小的模板[1 1 1 1 1 1 1 1 1 1 1; 4 5 6 7 8 24 8 7 6 5 4; 1 1 1 1 1 1 1 1 1 1 1]对泥岩夹层模型(图6)进行平滑，可以平滑掉长度40米及40米以下的泥岩夹层，平滑结果如图9所示。对比图8、图9和图6可以看出，新的平滑方法可以消除给定长度的噪点，平滑处理后的泥岩夹层的平均长度更接近模拟前预设的期望值，也就是说新的平滑方法效果明显。通过改变平滑窗口的长度能够保留长度在指定值以上的泥岩夹层，改进算法更具有针对性。
5. 结论
基于像元的建模方法如SISIM、SNESIM等经常被用于建立离散变量的模型，如微相、岩相、岩性模型等。由于算法本身的特点，会导致模拟结果中会出现一些不符合地质认识的离散的“噪点”，为了使模型更合理，通常需要做平滑后处理。Clayton教授提出的MAPS平滑方法在多数情况下均能取得较好的效果，但是当存在薄的地质体时，例如薄的泥岩夹层，该平滑方法会把线性或片状分布的夹层完全平滑掉。针对这种情况，本文提出了改进方法，根据地质体的特点设计平滑窗口，并根据需要平滑掉的最短泥岩长度设计相应的权值，改进后的平滑方法不仅能够去除“噪点”，而且能够克服原始MAPS方法把薄的地质体平滑掉的不足。通过对比分析发现，改进后方法得到的地质模型更接近地质真实情况。本文给出的算例是二维模型，提出的算法可以很容易扩展到三维空间。
致谢
本文是国家科技重大专项任务(辫状河储层构型建模算法模块编写及知识库应用系统研制，编号2016ZX05011-001-004)的部分研究成果。
NOTES
*通信作者。