基于Schwarz Christoffel变换的曲流河地质建模方法

姬安召， 王玉风， 张光生

(陇东学院能源工程学院， 庆阳 745100)

在油藏地质建模时，如何根据现有的井位数据建立出符合曲流河沉积特征的地质模型，这对油藏综合地质评价以及油藏数值模拟都具有重要意义。然而诸多学者从沉积学的角度对曲流河河道的变迁、点坝砂体内部结构解剖进行了大量研究，这些方法将从半定量与半定性的方面重构的曲流河沉积及演化的历史，使得沉积相的分布情况与砂体的构型与实际沉积演化情况更加吻合[1-2]。在建立曲流河地质建模时，这些点坝砂体控制约束地质模型的属性数据，按照地质统计学方法预测储层属性数据。为此，诸多学者对常规地质统计学方法进行改进，对精细构筑曲流河地质模型进行研究。在建立曲流河砂体过程中，将多点地质统计学方法与砂体平面展布图像相结合，曲流河砂体与夹层在训练图像与旋转数据体的双重因素约束情况下，重构了曲流河夹层合理发育规模及产状[3-4]。采取层次模拟、试验筛选、分级预测交互方式建立点坝砂体形态，然后利用多点地质统计学方法建立了点坝内部结构模型[5-6]。上述这些研究成果对构建曲流河形态及三维地质模型这有重要的意义。

在Schwarz Christoffel变换方面，诸多学者做了大量的研究。文献[7]研究了带状不规则封闭区域到规则矩形区域的映射关系，其核心解决方案需要借助中间的直线带状过渡区域来实现，这个直线带状过渡区域由复数域中的对数函数实现映射的，而规则矩形区域的映射函数是由复数域中的第一类雅可比椭圆函数实现的。在Schwarz Christoffel映射数值计算过程中，因为过渡的区域与实际区域，过渡区域与规则矩形区域差异形状差异大，因次需要选择合理的点位映射对应关系，若点位选择不恰当，则会引起点位的集聚现象，甚至导致计算结果不收敛[8]。从带状不规则封闭区域到直线带状过渡区域映射计算时，需要计算积分方程，该积分方程的积分起点与积分终点均为奇异点，因此解决奇异积分的计算尤为关键，文献[9]给出的解决方案，即采用高斯雅克比积分方法可以解决Schwarz Christoffel变换中的奇异积分计算问题。但高斯雅可比积分关键是要确定被积函数积分节点数量与积分节点的权值，文献[10-11]中对不规则封闭带状区域到直线带状过渡区域的Schwarz Christoffel映射关系的被积函数做了近似处理，得到了被积函数的形式，并给出了被积函数的节点及权值的计算方法。在不规则封闭带状区域到直线带状过渡区域的Schwarz Christoffel映射数值计算过程中，无论是边界点位的映射计算问题，还是映射区域内部点位的映射计算问题，都要受到循黎曼原理约束，一般计算思路是将通过积分方程的计算得到映射点位，然后建立映射点位前后的对应的函数关系，将该函数关系作为优化的目标来进行求解。然而常规的优化方法都是基于实变量函数而言，即使有约束条件，该约束条件可以转化为线性或非线性的等式或不等式约束条件，而二维平面上的Schwarz Christoffel映射问题的约束条件不但要遵循黎曼对应原理，还与计算过程中选择已知点位有关，因次想要用常规的线性或非线性的等式或不等式约很困难，在解决多边形区域到上半平面Schwarz Christoffel映射问题时，文献[12]提出了遵循黎曼对应原理的一维的实参数(上半平面x轴)与二维平面(多边形区域)点位的对应变化关系，文献[13-14]在此基础之上，建立复参数(二维平面)与实参数(一维空间点位的对应距离)的变换关系，解决了边界映射参数的约束问题。文献[15]提出了实际井位与矩形中井位的映射计算方法，这为Schwarz Christoffel变换的建模方法提供的理论计算方法。虽然目前商用的地质建模软件Petrel中提供了主方向可变的变差函数分析方法，但需要事先确定变差函数的主方向，并且以固定的数据编码格式输入到软件中。

借鉴前人的研究方法及思路，结合曲流河具有改道改向的特点，将曲流河轮廓看作不规则封闭带状区域，通过Schwarz Christoffel变换方法，将曲流河河道的方向映射为规则区域矩形的长边，垂直于河道方向映射为规则矩阵区域的短边。应用常规的地质统计学方法进行矩形区域油藏物性参数的预测时，就可以减少实际曲流河中因河流改道改向而造成物性分布方向多样性的影响。根据Schwarz Christoffel变换原理可知，该映射方法是可逆的，并且遵循保形映射，因此可以将矩形区域油藏模拟的基本参数，按照点位的对应关系可以还原到实际的油藏区域。

1 基本映射数学模型

在复平面w(不规则封闭带状区域)上有N(N>4)边形，它的顶点与内角分别为wb(k)和παk(k=1,2,…,N)。将直线带状过渡区域z边界上的点映射到不规则封闭带状区域w的Schwarz Christoffel映射公式[7]为

(1)

式(1)中：A为伸缩系数；C为映射中心，m；παk为复平面w多边形内角，rad；wb(k)为不规则封闭带状区域边界点，m；zb(k)为带状过渡区域边界点，m；f(ξ)为分段函数，具体表达式为

(2)

式(2)中：i为虚数单位；M为直线带状过渡区域下边界点的总个数，个；N为复平面w不规则封闭带状区域顶点的总个数，个；θ+为规则直线带状过渡区域左边的无限远点的角度，rad，则θ+=π；θ-为规则直线带状过渡区域右边的无限远点的角度，rad，则θ-=-π。对于式(1)的奇异积分方程，本文中采用高斯雅克比积分方法进行计算。

根据第一类椭圆函数关系，可以得出规则的矩形区域边界到直线带状过渡区域边界的Schwarz Christoffel映射关系，可表示[7]为

(3)

式(3)中：ub(k)为规则矩形区域边界点位，m；zb(k)为直线带状过渡区域的边界点位，m；l为第一类椭圆函数的映射模量；sn()为第一类椭圆函数。文献[7-8]中提出了不规则封闭带状区域到矩形区域的Schwarz Christoffel映射的数值计算问题。

2 曲流河内部点位与矩形区域内部点位的映射计算模型

在这个过程中，不规则带状封闭区域到直线带状过渡区域以及规则矩形到直线带状过渡区域的边界的映射点、伸缩系数都必须是已知的。在上述条件下，才能根据不规则封闭带状区域内部的点位计算对应的规则矩形区域内部的保形映射点位。具体计算关系可表示为

(4)

式(4)中：win(n)为距离wb(n)最近的不规则封闭带状区域内部的第n个点(win(n)为已知参数)，m；zin(n)为直线带状过渡区域内部的点位，是win(n)对应的映射点位(zb(k)为已知参数，zin(n)为待求解参数)，m；下标n为多边形内部已知点的编号。zin(n)的求解公式为

(5)

式(5)中：uin(n)为规则矩形区域内部的点位(uin(n)为待求解参数)，m；l为映射模量。式(4)的数值计算与式(1)的计算方法类似，但文献[7]中实参数与复参数的对应关建立复杂，根据文献[15]提供的方法，采用粒子群(partical swarm optimization, PSO)优化算法[16-17]求解二维空间平面点位映射问题，避免复杂的实参数与复参数的复杂关系的建立。式(5)采用牛顿迭代法对雅克比椭圆函数的进行计算。

3 映射模型的缩放比例问题

3.1 规则矩形区域顶点选方法的约定

根据不规则封闭带状区域的形状，将其边界映射到规则矩形边界时，首先在不规则封闭带状区域边界上选择两点，这两点被映射为规则矩形长边的两个顶点，然后在不规则封闭带状区域边界选择其他两点，这两点被映射为矩形短边的两个顶点。选择这四个点在不规则封闭带状区域中的形状尽可能与规则矩形区域的形状匹配。基于以上考虑，规定选择的这四点的编号依次为规则矩形的第一顶点、第二顶点、第三顶点和第四顶点。

3.2 规则封闭带状区域到带状过渡区域的映射比例的确定

根据油藏边界(不规则封闭带状区域边界)到直线带状过渡区域边界的对应关系，油藏边界上选择的第一顶点到第二顶点之间点被映射到直线带状过渡区域的下边界，油藏边界上选择的第三顶点到第四顶点之间点被映射到直线带状过渡区域的上边界。根据上述的映射对应关系，记油藏边界选择的第一顶点到第二顶点之间所有点相邻两点之间的距离之和为L1，记油藏边界选择的第三顶点到第四顶点之间所有点相邻两点之间的距离之和为L2，将L1与L2的平均值作为直线带状过渡区域x方向的伸缩系数。同理，记油藏边界选择的第二顶点到第三顶点之间所有点相邻两点之间的距离之和为L3，记油藏边界选择的第三顶点到第四顶点之间所有点相邻两点之间的距离之和为L4，将L3与L4的平均值作为直线带状过渡区域y方向的伸缩系数。

3.3 直线带状过渡区域到矩形区域的映射比例的确定

根据映射对应关系式(5)可知，直线带状过渡区域上下边界的长度与规则矩形区域的高度对应，直线带状过渡区域的高度与矩形区域的宽对应。根据映射模量l，对第一类完全椭圆函数积分进行计算，可得规则矩形区域的四顶点，计算第一类完全椭圆函数可得矩形的高度Kp(l′)和宽度2K(l)，具体的计算公式为

(6)

(7)

式中：K(l)为矩形宽度的一半；Kp(l′)为矩形的高度；其中l和l′的关系为l=1-l′。文献[18]中给出了第一类完全椭圆函数数值积分的具体计算方法。

4 实例分析

X油田位于鄂尔多斯盆地，截至2017年7月，研究区钻井78口，其中直井50口，水平井28口，目前采油井70口，累积采油24.48×104m3。X砂岩油藏为属于河流相沉积，属于小型的地层加背斜构造复合圈闭，油藏整体无统一的油水界面，含油面积79.6 km2。含油范围既受岩性尖灭线限制，也受构造控制。

4.1 曲流河边界到矩形油藏边界映射的计算

4.1.1 油藏边界映射点位的选择

X砂岩油藏局部区域属于曲流河沉积，河道开始由西南向东北发展，中间经过改道后又转向东方向。根据河道的大致走向，确定出研究区的边界，井位的分布如图1所示，油藏边界离散了38个点。按照3.1节给出的点位选择方法，确定了油藏边界四个顶点将被映射为规则矩形区域的区域顶点，37号点为第一顶点，16号点为第二顶点，17号点为第三顶点，36号点为第四顶点。

4.1.2 曲流河边界到矩形区域边界点位计算结果

根据文献[7]的求解思路，计算出图1的映射模量为5.577 054。根据文献[14]高斯雅可比积分方法以及积分路径的确定方法，采用高斯雅克比积分方法计算式(2)的积分，按照Levenberg Marquardt算法计算式(2)的积分方程组，得到直线带状过渡区域边界映射点位和规则矩形区域边界点位的分布情况，如图2所示。计算得到的伸缩系数A=3.93×106+1.29×106i。图2(a)中的直线带状过渡区域的下边界37～16号点被映射到规则矩形区域的高度的一条边上，直线带状过渡区域的高度方向点位被映射到规则矩形区域的宽度方向，但图2(b)中是没有考虑实际油藏与映射后油藏的比例缩放。

图1 油藏边界映射为矩形区域的顶点选择结果Fig.1 Selection result of the reservoir boundary mapped to the rectangular region vertex

图2 油藏边界-带状过渡区域边界-矩形边界映射图Fig.2 Reservoir boundary banded transition region boundary rectangular boundary map

4.1.3 精度评定结果

采用Levenberg Marquardt算法时，积分方程迭代计算的绝对误差与迭代次数的关系如图3曲线1所示。在迭代前11次中，迭代次数与绝对误差呈线性关系下降。迭代次数从11～29到次时，绝对误差下降缓慢，下降到10-2数量级，但迭代次数达到30次以后，绝对误差下降很快，说明在接近映射点位真值的局部范围内，Levenberg Marquardt算法在计算不规则封闭带状区域到直线带状过渡区域边界点位映射时的效果较好。在X油藏实例的计算中，采用Levenberg Marquardt算法计算积分方程式(2)的绝对误差达到6.959×10-14m。

牛顿迭代法计算复参数雅克比椭圆函数的绝对误差与迭代次数的关系如图3曲线2所示。从曲线2可以看出，牛顿迭代法计算是的效率很高，在迭代5次时，边界上38个点的绝对误差达到10-14数量级，主要原因是雅克比椭圆函数的导函数是存在的，在计算的区域中不存在奇点。采用牛顿迭代法计算雅克比椭圆函数的绝对误差达到1.776×10-15m。最后通过综合精度评价，得出整理映射计算的绝对误差为6.451×10-10m。

图3 油藏边界-带状过渡区域-矩形边界映射计算绝对误差与迭代次数关系Fig.3 Relationship between absolute error and iteration times of reservoir boundary banded transition region rectangular boundary mapping calculation

4.2 油藏中的井位到矩形区域中的井位映射计算

通过上述计算得到直线带状过渡区域中38口井的井位，如图4所示。通过式(4)和式(5)的计算，得到X油藏38口井在矩形区域中的位置，如图5所示。这里特别需要说明的是，从油藏边界到带状过渡区域边界在映射过程中存在这井位的缩放问题，同时从直线带状过渡区域到矩形区域映射也存在井位的缩放问题。但这两者对井位缩放比例是不一样的，前者是井位的整体缩放，而后者在x轴方向和y轴方向缩放的比例是不一致的。图5中给出的矩形区域的映射井位是没有考虑x轴方向和y轴方向的伸缩系数，但在后续矩形区域油藏地质建模时需要考虑比例的缩放问题。

通过表1可以看出，在不考虑缩放比例情况下，计算出X砂岩油藏的井位在直线带状过渡区域中的点位绝对误差也在10-11m以下。这样的精度是能够满足工程计算的需要，同时也说明本文中所给出的PSO算法在求解积分方程式(4)是可行的。通过表1中的第3、第7列可以看出，在设定积分方程点位精度为10-10m情况下，迭代次数没有超过设下，通过多种随机性和确定性模拟的方法对比，最后选用了沉积相控制下的序贯高斯模拟方法，模拟了X砂岩油藏在映射后规则矩形区域中的孔隙度的分布，如图6所示，孔隙度沿着规则矩形区域的高度方向(即映射后的河道方向)呈现条带状的分布。定500次，最大的PSO迭代次数61次。在表1中第4、第8列反映了牛顿迭代法计算雅克比椭圆函数的绝对误差，其绝对误差在10-15数量级以下。

图4 从油藏内部井位到带状过渡区域映射井位分布Fig.4 Mapping well location distribution from internal well location of reservoir to banded transition area

图5 带状过渡区域到矩形油藏映射井位示意图Fig.5 Schematic diagram of mapping well location from banded transition area to rectangular reservoir

表1 PSO算法求解积分方程与牛顿迭代法求解雅克比椭圆函数绝对误差数据表Table 1 Absolute error data of solving integral equation by PSO algorithm and Jacobian elliptic function by Newton iterative method

4.3 矩形区中地质模型建立及物性参数的分析

根据不规则封闭带状区域到带状过渡区域的映射比例的确定方法，计算得到油藏边界37号点到16号点之间的距离为19 863.6 m，36号点到17号点之间的距离为24 596.6 m，二者之间的平均值为22 230.1 m，即x方向的伸缩系数为22 230.1。计算得到油藏边界37号点到36号点之间的距离为2 498 m，16号点到17号点之间的距离为2 239 m，二者之间的平均值为2 364 m，即y方向的伸缩系数为2 364。根据式(6)和式(7)计算得到直线带状过渡区域到矩形区域的映射比例参数K=1.57，Kp=18.9。因此在建立规则矩形区域的构造模型、地层模型时，只要涉及平面位置的参数，需要对平面位置参数的x与y值乘以对应的伸缩系数即可。例如对映射后规则矩形区域的井位及油藏边界点数据，对x坐标乘以22 230.1，再除以2K；对y坐标乘以2 364，再除以Kp，通过上述处理，就得到考虑伸缩系数的矩形区域的数据。在地层及沉积相的控制。

4.4 矩形区域地质模型的还原

4.4.1 原地质模型

根据映射前的原始数据，在地层与沉积相的约束性下，使用序贯高斯模拟技术，模拟了X砂岩油藏原始的孔隙度参数场的分布情况，如图7所示。在模拟的过程中，拟合了南西北东向的变差函数，因次模拟的孔隙度分布主要主变程方向，而河道是沿着砂体是随着河道的改向而变化的，这与实际的砂体分布有差异。根据图7(b)～图7(d)可以看出，研究区2、3、5各个小层都呈现同样的趋势。

图6 矩形区域中X砂岩油藏岩心实验孔隙度分布模拟结果Fig.6 Simulation results of porosity distribution of X sandstone reservoir in rectangular area

图7 X砂岩油藏映射变换前的孔隙度分布Fig.7 Porosity distribution of X sandstone reservoir before mapping transformation

4.4.2 矩形油藏建立的注意事项

(1)采用矩形区域映射数据进行建模时，可以将上下左右四个边界可以向中心移动一个网格步长的距离。这样做的目的是保证在矩形中形成的网格完全位于映射边界内部，以免在后续网格数据还原过程中有网格节点处于映射边界的外部。

(2)将PETREL软件建立的网格数据导出为ECLIPSE格式的十进制数据时，应包含角点网格数据的关键字COORD。PETREL软件默认的网格编号是从左上角开始，即系统默认的坐标原点在左上角，因此导出的x方向的数据没有问题，而导出的y方向的数据均为负值。实际矩形边界的映射数据y值均为正值，实映射数据以北方向作为y轴的正方向，而PETREL软件中的y轴的正方向为南方向，因此y值出现负值情况。因此必须取消y值前面的符号才能进行映射计算。

(3)模型数据缩放比例的还原。图8为考虑缩放比例以后建立的地质模型，因此在矩形油藏网格及井位的数据还原式，对x、y坐标数据分别乘以矩形区域x、y方向的缩放系数2K与Kp，然后再除以带状区域的x、y方向的缩放系数2 462.986和22 230.606。

(4)特殊点的处理，若矩形区域的映射数据包含0值，则将其重置为一个很小的数，如0.000 001，以免在积分是出现奇异值问题。

将矩形区域的地质模型通过Schwarz Christoffel逆变换，将属性参数还原到原油藏的实际边界区域中，图8反映了还原后的孔隙度分布情况。通过图8可以看出，进行映射变换处理，孔隙度变化方向基本沿着河道方向，符合基本的地质规律。矩形区域中的规则正交网格也被映射到实际油藏的模型中，也保持了一定的正交性。

图8 从矩阵区域还原到实际油藏边界的X砂岩油藏孔隙度分布Fig.8 Porosity distribution of X sandstone reservoir restored from matrix area to actual reservoir boundary

5 结论

(1)采用粒子群优化算法求解了从不规则油藏内部点位到直线带状过渡区域内部点位的数值映射计算，粒子群算法迭代的最大次数为61次，映射井位的绝对误差在10-11m以下,说明粒子群算法在计算不规则油藏内部点位到直线带状过渡区域内部点位是可行的。

(2)采用牛顿法对雅克比椭圆函数进行求解得到矩形区域的映射点位，映射的点位绝对误差在10-15数量级以下，说明牛顿法在计算雅克比椭圆函数是可行的。

(3)从油藏边界到直线带状过渡区域边界存在井位的缩放，但缩放的大小与伸缩系数A的大小有关，是对井位的整体缩放。从直线带状过渡区域到规则矩形区域映射对井位的缩放，在x轴方向和y轴方向缩放的比例是不一致的，主要取决于雅克比椭圆函数映射模量的大小，而映射模量的大小取决于矩形边界的四个顶点。

(4)根据不规则封闭带状区域到规则矩形区域Schwarz Christoffel变换基本原理，建立曲流河模型边界映射到规则矩形区域边界的数学模型。以X砂岩油藏为例，选择了该区域中38口井，建立的映射前后的地质模型，并对物性参数进行了分析。通过分析对比，采用Schwarz Christoffel映射后建立的地质模型及分析的物性参数更加符合地质规律。