CN113419280B - 基于改进椭圆拟合的叠前裂缝密度估算方法 - Google Patents

Abstract

本发明涉及一种基于改进椭圆拟合的叠前裂缝密度估算方法，包括以下步骤：取叠前道集一个时间点的数据，分别取方位角和振幅为椭圆拟合的坐标，利用拉依达法则剔除异常值，建立正演坐标点与实际坐标点之间的最小平方差和的目标函数，同时利用最小二乘线性拟合算法计算椭圆拟合系数，作为改进高斯牛顿优化迭代算法的初始值，利用迭代算法得到最优椭圆拟合系数，据此计算裂缝密度，完成整个叠前道集的裂缝密度估算。本发明估算方法具有计算效率高、拟合精度高、对噪声的鲁棒性强等优点，有效解决了现有基于线性估计最小二乘方法不能达到最优解，使用存在一定局限性问题。

Description

基于改进椭圆拟合的叠前裂缝密度估算方法

技术领域：

本发明涉及一种叠前裂缝密度估算方法，尤其是涉及一种基于改进椭圆拟合的叠前裂缝密度估算方法。

背景技术：

随着油气勘探技术的提高，裂缝性油气藏已经成为了一个重要的勘探新领域，储层裂缝的识别与裂缝的准确描述与预测是裂缝性储层勘探开发的关键。

目前，按照资料类型可以将裂缝预测技术分为两大类，一类是叠后预测技术，另外一类是叠前预测技术。目前叠后裂缝预测技术有相干分析技术、图像处理技术、曲率分析技术、倾角检测技术、盈利分析技术等。它们通常以叠后地震资料为基础输入，研究的是整个网状裂缝系统中那些较大尺度的裂缝，从预测效果上看预测是网状断裂系统的骨架和轮廓，在平面上呈线性分布，数量较少；而对较小的尺度、数量更庞大的小型或者微裂缝系统，叠后裂缝预测技术效果不好。

叠后裂缝预测的精度受限于叠后地震资料的信息量较小，缺乏偏移距和方位角的信息，而叠前地震资料具有数据量大、信息量大的特点，可以用以检测包含叠前观测方位-偏移距二维时空域中的微断裂的地震响应信息来检测断裂体系。目前，叠前裂缝检测主要是基于纵波反射波资料。在均匀各向同性介质中，振幅等地震波动力学及运动学属性均无方向和方位变换，在振幅方位角平面内的拟合图形是圆形。当存在裂缝性气藏时，由于裂缝油气藏以及裂缝的存在，使得介质具有各向异性，其平面属性拟合图形呈椭圆图形，因此可以通过方位角振幅属性进行椭圆拟合来进行裂缝预测。

目前常见的椭圆拟合算法是基于最小全局误差的最小二乘方法，该算法基于拟合预测坐标与实际坐标的最小方差建立目标函数，该算法原理简单，计算效率快的优点，最小二乘法是线性估计,已经默认了是线性的关系,并不能达到最优解，因此使用具有一定的局限性。

发明内容：

为了解决上述技术问题，本发明提供了一种基于改进椭圆拟合的叠前裂缝密度估算方法，所述方法包括以下步骤：

(1)取工区叠前道集一个时间点的数据，分别取方位角和振幅为椭圆拟合的坐标，

建立以方位角为X轴，以振幅为Y轴的笛卡尔二维坐标系，将笛卡尔坐标系变换至极坐标系中，变换式如下：

其中，上面两个式子表示一个点的极坐标变换式，x₁,x₂分别是原笛卡尔坐标系的值，z₁,z₂是椭圆的中心坐标，r是极坐标半径，

是坐标夹角；

(2)利用拉依达法则剔除异常值：

地震振幅数据符合正态分布，符合拉依达法则剔除异常值的应用范围，采用拉依达法则可以剔除大的振幅异常点，3σ准则又称为拉依达准则，其中σ是标准差，拉依达法则先假设一组观测数据只含有随机误差，对其进行计算处理得到标准偏差，按一定概率确定一个区间，超过3倍的标准差的误差点，就属于大的误差观测点，含有该误差的数据应予以剔除；

(3)建立正演坐标点与实际坐标点之间的最小平方差和的目标函数：

其中，m为数据组数，a为椭圆的短轴，b为椭圆长轴，α为坐标变换的旋转角度，每个观测点的坐标夹角是

g是椭圆正演坐标与实际坐标的误差；

(4)利用最小二乘线性拟合算法计算椭圆拟合系数，作为改进高斯牛顿优化迭代算法的初始值，利用迭代算法得到最优椭圆拟合系数：

对于(3-1)式，一共有m+5个待求系数，分别是a,b,z₁,z₂,α,

i＝1,2,…,m，其中a是椭圆的短轴，b是长轴，z₁,z₂是椭圆的中心坐标，α是坐标变换的旋转角度，

是第i个观测点的坐标夹角，分别求取待求系数的偏导数：

建立偏导矩阵(4-7)

设待求系数的集合是

具体步骤如下：

a)以最小二乘拟合的的解为迭代优化的初始值M₀，带入到(3-1)式得到观测误差ΔD＝[g₁,g₂,…,g_m]；

b)求初始参数正演与实际观测值之间的误差梯度ΔG＝J^TΔD；

c)计算海塞矩阵H＝J^TJ，此步算法可采用分块并行进行计算，加快运算速度；

d)计算拟合系数的扰动量ΔM＝H^-1ΔG，本发明算法中，计算海塞矩阵的逆矩阵时，利用海塞矩阵的对称性质，采用了三角矩阵分解递归迭代的快速解法，具体为：

①对拟合系数的扰动量算式进行变形得到：

HΔM＝ΔG (4-7)

②利用Cholesky分解法得H＝LL^T，则L是下三角矩阵；

③式(4-7)可变为

LL^TΔM＝LS (4-8-1)

LS＝ΔG (4-8-2)

式(4-8-2)采用递归迭代快速计算出S；然后式(4-8-1)可以简化为：

L^TΔM＝S (4-8-3)

同样可以采用快速递归迭代算法求解出椭圆拟合系数扰动量ΔM；

e)更新椭圆拟合的系数M＝M₀+ΔM；

f)M₀＝M，重复步骤1)～步骤5)，直至误差收敛域误差门槛值内；

(5)椭圆拟合系数计算裂缝密度：

设裂缝密度属性为F_Den＝a/b-1，其中F_Den是裂缝密度；

(6)完成步骤(2)～步骤(5)遍历叠前道集所有采样时间点，获得裂缝密度估算结果。

本发明的基于改进椭圆拟合的叠前裂缝密度估算方法，取得了以下技术效果：

(1)本发明是一种非线性全局优化迭代算法，具有拟合精度高、对噪声的鲁棒性强的优点；

(2)本发明改进了高斯牛顿算法的不足，采用了局部最优的最小二乘解作为迭代算法的初始值，采用海塞矩阵的三角矩阵分解递归求解代替原方法的海塞矩阵求逆求解过程，一阶偏导矩阵运算过程中可以采分块并行，故本发明算法计算效率快，拟合精度高；

(3)本发明算法计算过程中采用了偏导矩阵作为椭圆拟合系数扰动量的系数，间接实现了非线性映射关系，且偏导矩阵计算过程等价于拟合系数矩阵主特征值提取过程，故算法对噪声具有一定的鲁棒性。

本发明与上述背景技术相比较，还存在如下有益效果：

常规基于线性拟合的最小二乘方法裂缝密度估算方法相比，无法获得裂缝密度的最优估算，而将本发明“基于改进椭圆拟合的叠前裂缝密度估算”应用于裂缝性储层裂缝密度估算中，以传统的线性最小二乘拟合结果为初始值，并采用改进的高斯牛顿算法进行最优拟合，可以获得裂缝密度的最优估算，具有拟合精度高、抗噪能力强等优点。

附图说明：

图1为本发明的技术流程图；

图2a-2c为本发明利用拉依达法则对观测异常点进行了剔除前后的椭圆拟合对比图。其中图2a是对异常点剔除之前的椭圆拟合，图2b是对异常点剔除之后的椭圆拟合，图2c的表格是用于对异常点剔除的数据；

图3a-3c为本发明方法的效果图，其中，图3a是模拟角道集数据进行了椭圆拟合试验，其中X轴是方位角，Y轴是地震振幅，图中圆圈是原始叠前道集的某一时刻数据在方位角振幅二维平面上的分布；黑色虚线是采用最小二乘拟合得到的椭圆正演图形；黑色实线是本发明改进的高斯牛顿迭代法拟合得到的椭圆正演图形；图3b的数据同图2的数据相同；表3c列出的是图3a中的数据。

具体实施方式：

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图，以大庆油田松辽盆地三维实际地震工区为例，对本发明作进一步的详细描述。

本发明的基于改进椭圆拟合的叠前裂缝密度估算方法，包括以下步骤：

(1)取工区叠前道集一个时间点的数据，分别取方位角和振幅为椭圆拟合的坐标：

建立以方位角为X轴，以振幅为Y轴的笛卡尔二维坐标系，由于椭圆拟合在笛卡尔坐标系中会出现坐标重合的现象，所以将笛卡尔坐标系变换至极坐标系中，变换式如下：

上面两个式子表示一个点的极坐标变换式，x₁,x₂分别是原笛卡尔坐标系的值，z₁,z₂是椭圆的中心坐标，r是极坐标半径，

是坐标夹角。

(2)利用拉依达法则剔除异常值：

地震振幅数据符合正态分布，符合拉依达法则剔除异常值的应用范围，采用拉依达法则可以剔除大的振幅异常点。3σ准则又称为拉依达准则，其中σ是标准差，拉依达法则先假设一组观测数据只含有随机误差，对其进行计算处理得到标准偏差，按一定概率确定一个区间，超过3倍的标准差的误差点，就属于大的误差观测点，含有该误差的数据应予以剔除。

图2是利用拉依达法则对观测异常点进行了剔除前后的椭圆拟合对比图，其中2(a)是对异常点剔除之前的椭圆拟合，2(b)是对异常点剔除之后的椭圆拟合，2(c)的表格是用于对异常点剔除的数据；对比分析，异常点的剔除后再进行椭圆拟合更加符合数据特征。

(3)建立正演坐标点与实际坐标点之间的最小平方差和的目标函数：

设一共有m组数据，椭圆的短轴是a，长轴是b，坐标变换的旋转角度是α，每个观测点的坐标夹角是

g是椭圆正演坐标与实际坐标的误差，则建立椭圆拟合正演与观测的平方差的目标函数(3-1)。

(4)利用最小二乘线性拟合算法计算椭圆拟合系数，作为改进高斯牛顿优化迭代算法的初始值，利用迭代算法得到最优椭圆拟合系数：

对于(3-1)式，一共有m+5个待求系数，分别是a,b,z₁,z₂,α,

i＝1,2,…,m，其中a是椭圆的短轴，b是长轴，z₁,z₂是椭圆的中心坐标，α是坐标变换的旋转角度，

是第i个观测点的坐标夹角，分别求取待求系数的偏导数：

建立偏导矩阵(4-7)

本发明算法采用了改进的高斯牛顿优化算法进行椭圆拟合，设待求系数的集合是

具体步骤如下：

a)以最小二乘拟合的的解为迭代优化的初始值M₀，带入到(3-1)式得到观测误差ΔD＝[g₁,g₂,…,g_m]；

b)求初始参数正演与实际观测值之间的误差梯度ΔG＝J^TΔD；

c)计算海塞矩阵H＝J^TJ，此步算法可采用分块并行进行计算，加快运算速度；

d)计算拟合系数的扰动量ΔM＝H^-1ΔG，本发明算法中，计算海塞矩阵的逆矩阵时，利用海塞矩阵的对称性质，采用了三角矩阵分解递归迭代的快速解法，具体为：

①对拟合系数的扰动量算式进行变形得到：

HΔM＝ΔG (4-7)

②利用Cholesky分解法得H＝LL^T，则L是下三角矩阵；

③式(4-7)可变为

LL^TΔM＝LS (4-8-1)

LS＝ΔG (4-8-2)

式(4-8-2)采用递归迭代快速计算出S；然后式(4-8-1)可以简化为：

L^TΔM＝S (4-8-3)

同样可以采用快速递归迭代算法求解出椭圆拟合系数扰动量ΔM。

e)更新椭圆拟合的系数M＝M₀+ΔM；

f)M₀＝M，重复步骤1)～步骤5)，直至误差收敛域误差门槛值内。

(5)椭圆拟合系数计算裂缝密度：

设裂缝密度属性为F_Den＝a/b-1，其中F_Den是裂缝密度。

(6)完成步骤(2)～步骤(5)遍历叠前道集所有采样时间点，得到裂缝密度估算结果。

图3本发明方法的效果图，图3(a)是模拟角道集数据进行了椭圆拟合试验，其中X轴是方位角，Y轴是地震振幅，图中圆圈是原始叠前道集的某一时刻数据在方位角振幅二维平面上的分布；黑色虚线是采用最小二乘拟合得到的椭圆正演图形；黑色实线是本发明改进的高斯牛顿迭代法拟合得到的椭圆正演图形；3(b)的数据同图2的数据相同；3(c)是图3(a)的数据。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于改进椭圆拟合的叠前裂缝密度估算方法，其特征在于，包括以下步骤：

(1)取工区叠前道集一个时间点的数据，分别取方位角和振幅为椭圆拟合的坐标，

建立以方位角为X轴，以振幅为Y轴的笛卡尔二维坐标系，将笛卡尔坐标系变换至极坐标系中，变换式如下：

其中，上面两个式子表示一个点的极坐标变换式，x₁,x₂分别是原笛卡尔坐标系的值，z₁,z₂是椭圆的中心坐标，r是极坐标半径，

是坐标夹角；

(2)利用拉依达法则剔除异常值：

地震振幅数据符合正态分布，符合拉依达法则剔除异常值的应用范围，采用拉依达法则剔除大的振幅异常点，3σ准则又称为拉依达准则，其中σ是标准差，拉依达法则先假设一组观测数据只含有随机误差，对其进行计算处理得到标准偏差，按一定概率确定一个区间，超过3倍的标准差的误差点，就属于大的误差观测点，含有该误差的数据应予以剔除；

(3)建立正演坐标点与实际坐标点之间的最小平方差和的目标函数：

其中，m为数据组数，a为椭圆的短轴，b为椭圆长轴，α为坐标变换的旋转角度，每个观测点的坐标夹角是

g是椭圆正演坐标与实际坐标的误差；

(4)利用最小二乘线性拟合算法计算椭圆拟合系数，作为改进高斯牛顿优化迭代算法的初始值，利用迭代算法得到最优椭圆拟合系数：

对于(3-1)式，一共有m+5个待求系数，分别是

其中a是椭圆的短轴，b是长轴，z₁,z₂是椭圆的中心坐标，α是坐标变换的旋转角度，

是第i个观测点的坐标夹角，分别求取待求系数的偏导数：

建立偏导矩阵(4-7)

设待求系数的集合是

具体步骤如下：

a)以最小二乘拟合的解为迭代优化的初始值M₀，带入到(3-1)式得到观测误差ΔD＝[g₁,g₂,…,g_m]；

b)求初始参数正演与实际观测值之间的误差梯度ΔG＝J^TΔD；

c)计算海塞矩阵H＝J^TJ，此步算法采用分块并行进行计算，加快运算速度；

d)计算拟合系数的扰动量ΔM＝H^-1ΔG；

e)更新椭圆拟合的系数M＝M₀+ΔM；

f)M₀＝M，重复步骤a)～步骤e)，直至误差收敛域误差门槛值内；

(5)椭圆拟合系数计算裂缝密度：

设裂缝密度属性为F_Den＝a/b-1，其中F_Den是裂缝密度；

(6)完成步骤(2)～步骤(5)遍历叠前道集所有采样时间点，获得裂缝密度估算结果。

2.如权利要求1所述的叠前裂缝密度估算方法，其特征在于，其中步骤(4)中，计算拟合扰动量的计算中，计算海塞矩阵的逆矩阵时，利用海塞矩阵的对称性质，采用了三角矩阵分解递归迭代的快速解法，具体为：

①对拟合系数的扰动量算式进行变形得到：

HΔM＝ΔG (4-8)

②利用Cholesky分解法得H＝LL^T，则L是下三角矩阵；

③式(4-8)可变为

LL^TΔM＝LS (4-8-1)

LS＝ΔG (4-8-2)

式(4-8-2)采用递归迭代快速计算出S；然后式(4-8-1)简化为：

L^TΔM＝S (4-8-3)

同样采用快速递归迭代算法求解出椭圆拟合系数扰动量ΔM。

Images