CN100543497C - 核测井数据多尺度处理方法 - Google Patents

Abstract

本发明涉及一种核测井数据多尺度处理方法，包括以下步骤：1)读入测井数据；2)对测井数据进行滤波；3)对滤波后的数据进行归一化处理；4)对经过滤波和归一化处理后测井数据进行多尺度分解、融合和重构；5)显示曲线和图像。本发明解决了现有技术只能在一维深度域中对采集到的数据进行处理和解释的局限性，本发明方法可在多维空间反映地层油、气、水分布的多尺度特性，提高了描述剩余油气饱和度分布的精确度。

Description

核测井数据多尺度处理方法

技术领域

本发明涉及一种核测井数据处理和分析方法。

背景技术

多尺度估计理论及其应用研究已成为信息技术的重要分支，是科技界的一个热点问题。测井所研究的对象、过程都具有多尺度特性：油藏、储层、岩心、岩屑、薄片、孔隙、裂缝等实体，其空间尺度分布从数千米到几微米；地层的电、声和核物理参数也构成了多维空间；测井系统本身是一个典型的多类传感器复杂系统。

核测井是地球物理测井的重要分支，是测井复杂系统中的一个子系统。每一种核测井仪器都要安装一个或多个同类或非同类传感器，每个传感器(探测器)具有各自的探测范围，可在不同尺度上对地层的核物理特性进行观测。但直到目前尚未见到此类研究成果或专利发布。

核测井采集的原始数据主要有三种结构：第一种，每个深度点测一个数，在一维空间中构成一条曲线，如自然伽马总计数率；第二种，每个深度点测一幅能谱，在深度-能量域中构成一幅二维数字图像，如各种伽马能谱测井；第三种，每个深度点测一幅时间谱，在深度-时间域中构成一幅二维数字图像，如热中子寿命测井。这三种类型的数据都包含多种尺度的信息，而现有的方法只在一维深度域中对所测参数进行处理和解释，忽略了地层参数在不同尺度上的响应，不能反映剩余油分布的非均匀性。

发明内容

本发明提出了一种核测井数据多尺度处理方法，解决了现有技术只能在一维深度域中对采集到的数据进行处理和解释的局限性，本发明方法可在多维空间反映地层油、气、水分布的多尺度特性，提高了描述剩余油气饱和度分布的精确度。

本发明的技术解决方案是：

一种核测井数据多尺度处理方法，包括以下步骤：

1]读入测井数据：

所述的测井数据包括以下数据中的至少一组：

①自然伽马测井：自然伽马总计数率、自然伽马能谱和铀、钍、钾含量；

②散射伽马测井：散射伽马总计数率、散射伽马能谱、密度和岩性指数；

③同位素中子源测井：中子计数率、俘获伽马计数率、俘获伽马能谱和孔隙度；

④脉冲中子测井：热中子总计数率和热中子时间谱；脉冲中子测井俘获伽马总计数率、俘获伽马时间谱和能谱；脉冲中子测井非弹伽马总计数率和能谱；中子活化伽马计数率、伽马时间谱和能谱；

2]对测井数据进行滤波：

对读入的测井数据进行深度域的纵向滤波和沿与深度轴垂直的时间轴或能量轴的横向滤波；所述纵向和横向滤波的方法包括Kalman滤波、多点拟合滤波、多点平滑滤波或多尺度滤波；

3]对滤波后的数据进行归一化处理：

所述归一化处理是指计算每个采样点的归一化测井值：

$f\′=\frac{f-f_{\min }}{f_{\max }-f_{\min }}\×100$

式中，f_min为测井曲线的最小值；f_max为测井曲线的最大值；f为测井曲线在每个采样点的实测数据；

4]对经过滤波和归一化处理后测井数据进行多尺度分解、融合和重构：

所述对测井数据进行融合和重构包括在数据、物理参数和地质参数三个层次上分别对单条测井曲线和多条测井曲线进行多尺度分解、融合和重构；

所述对单条测井曲线进行多尺度分解、融合和重构的步骤包括：

4.1.1]选择合适的多尺度分解法对测井曲线进行多尺度分解，所述横向时间轴上对时间谱多尺度分解的方法包括多指数拟合法和小波分解法；所述横向能量轴上对能量谱多尺度分解的方法包括解谱法和小波分解法；所述纵向深度轴上对总计数率、时窗计数率、能窗计数率、计数率比值、物理参数和地质参数测井曲线的多尺度分解的方法是指小波分解法；所述单条测井曲线包括：总计数率、时窗计数率、能窗计数率、物理参数或地质参数；

4.1.2]做出小波时频色谱图，通过颜色的变化显示出不同深度和不同尺度下小波系数的大小，进而使内部的能量聚集与分布得以清晰展示；

4.1.3]以小波模极大重构算法为基础，选用低频系数加权、高频系数取绝对值较大的融合规则，在各个层次上对测井数据进行小波重构得到融合小波金字塔；

4.1.4]对融合后的小波金字塔进行小波逆变换实现各个层次的数据重构；

4.1.5]采用熵、均值和方差作为评判标准，对融合数据进行定量评价：均值、方差数据反映峰值信息；信息熵数据反映空间细节信息；峰值信噪比、相关系数数据反映频谱信息；

所述对多条测井曲线的数据进行多尺度分解、融合和重构的步骤包括：

4.2.1]对多条原始测井曲线分别进行多尺度分解，即对原始数据分别进行低、高通滤波，使多列原始数据各自分解为含有不同频率成分的低频逼近信号和高频细节信号；所述横向时间轴上对时间谱多尺度分解的方法包括多指数拟合法和小波分解法；所述横向能量轴上对能量谱多尺度分解的方法包括解谱法和小波分解法；所述纵向深度轴上对总计数率、时窗计数率、能窗计数率、计数率比值、物理参数和地质参数测井曲线的多尺度分解的方法是指小波分解法；所述多条测井曲线包括不同源距探测器的总计数率、时窗计数率、能窗计数率、物理参数和地质参数；

4.2.2]根据需要，对低频逼近信号重复上面的过程，建立各数据的小波塔形分解；

4.2.3]对各个分解层进行融合处理，最终得到融合小波金字塔；

4.2.4]对融合后的小波金字塔通过小波逆变换法进行数据重构，所得的重构数据即为该多列测井曲线经融合重构的数据；

5]显示曲线和图像

所述的曲线和图像包括以下曲线和图像中的至少一种：

①原始数据曲线，包括滤波后的总计数率、时窗计数率、能窗计数率、计数率比值；

②经标准井检验选定的，多尺度分解得到的对含油饱和度或岩性分辨能力好的分量曲线；

③经标准井检验选定的，多尺度融合得到的对含油饱和度或岩性分辨能力好的融合、重构曲线；物理和地质参数随深度变化的曲线，包括热中子寿命、宏观俘获截面、碳氧比、含水饱和度；深度-时间、深度-中子寿命或深度-伽马能谱二维图像；

上述核测井数据多尺度处理方法还可包括进行地质解释的步骤：

所述进行地质解释的方法包括以下方法中的至少一种：

①用自然伽马曲线和自然电位融合曲线判定岩性；

②根据粗尺度信号的变化划分厚层和层组，根据细尺度信号的变化对厚层和层组进行细分；

③用近、远探测器计数率比值曲线指示孔隙度的变化；

④用深度-时间、深度-中子寿命二维图像，判定井眼流体和近井区地层中子寿命的径向分布，经统计分析和标准井刻度建立不同岩性油、气、水层的识别模式，作为资料解释的依据；

⑤根据热中子总计数率与俘获伽马总计数率的差别和热中子寿命或宏观截面色谱图，定性、定量识别和评价油、气、水层；

⑥用热中子宏观截面∑或碳氧比可计算油气水饱和度，并可获得储层在开发过程中剩余油饱和度分布的变化；

⑦对各种尺度上的低频和高频信号进行观察、分析和比较，可确定各个尺度的低频和高频信号与地层的空间和物理参数之间的相关性，依此可对不同级别的沉积旋回、层序界面和沉积相做定性和定量研究，并可用于厚层细分和描述剩余油气饱和度空间分布的细节做精细描述。

上述多尺度滤波的步骤包括：

2.1]选取具有一定的紧支撑性、对称性和平滑性的正交小波基；

2.2]选择一个小波并确定其分解层次N，然后对数据进行N层小波分解；

2.3]对第1层到第N层的每一层高频系数选择一个阈值进行软阈值量化处理；

2.4]根据小波分解的第N层低频系数和经过量化处理后的第1层到第N层的高频系数，进行测井数据的重构，重构后的数据已消除了统计不确定性和其它高频噪声；

所述对数据进行N层小波分解的步骤包括：

2.2.1]在某尺度i上，对给定的核测井信号序列 $x(i,k)\∈V_i\⋐l^2\left(Z\right)$      (k∈Z)，通过一个脉冲响应为h(k)的低通滤波器，获得粗尺度上(低频段)的平滑信号：

x_V(i-1，k)∈V_i-1

x_V(i-1，l)＝∑h(2l-k)x(i，k)

2.2.2]信号x(i，k)在低通滤波器中丢失的“细节信号”由x(i，k)通过一个脉冲响应是g(k)的高通滤波器得到细尺度(高频段)的细节信号：

x_D(i-1，k)∈D_i-1

x_D(i-1，l)＝∑g(2l-k)x(i，k)

下标D表示x(i，k)在细节信号空间D_i-1上的投影；

2.2.3]对平滑信号逐次进行N层的多尺度分解，得到2N个不同的频带，其中包含N个高频信号和N个低频信号，从低频信号中获得测井曲线中包含的地层概貌信息，从高频信号中获得地层的细节；

所述步骤4]中的小波模极大重构算法的步骤包括：

设测井仪器接采集到的信号为f(x)，Morlet连续小波变换定义为：

$W_f(a,b)=<f,{\mathrm{\&psi;}}_{\mathrm{ab}}>=\frac{1}{\sqrt{\left|a\right|}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}f\left(x\right)\mathrm{\&psi;}\left(\frac{x-b}{a}\right)\mathrm{dx}$

式中：Morlet小波基函数为

$\mathrm{\&psi;}\left(x\right)=e^{-x^2/2}\&CenterDot;e^{i{\mathrm{\&omega;}}_{0}x}$

a和b分别为尺度因子和位移因子；

对信号进行小波分析，可以得到不同尺度a的评估信号在不同空间段的大量系数，即相应的小波系数W_f(a，b)，并用色谱图的方式表示；在色谱图中，用颜色的变化来表示系数的大小；可选择得到某一尺度a时的小波系数曲线，这条曲线能够直观显示小波系数与被分析信号之间的相似程度，且可通过控制a的大小，使小波系数曲线和某些参数具有一定的相关性。

本发明方法的优点是：

(1)单条测井曲线经多尺度分解，从深度域转换到深度-尺度域，再经融合重构可消除噪声，提高对岩性和油气水饱和度的分辨能力。

(2)本发明放弃了用中子寿命或∑单一物理参数求取含油饱和度的传统方法，而将热中子寿命τ(或∑)看成是反映井眼-地层介质系统物理性质的一种尺度，将测得的时域中的数据转换到寿命域，实现了尺度的连续分解，时间道的分段处理就是对不同尺度范围内的信息重构。

(3)在深度-时间域和深度-寿命域中对核测井数据做多尺度分解和在选定的尺度域中实现数据融合重组和可视化处理，可在不同尺度或尺度组合中分析地层的精细结构，更准确地识别剩余油气和水淹层，并可压制统计涨落。

(4)用脉冲中子双谱饱和度测井仪，可同时采集热中子时间谱、俘获伽马时间谱、俘获伽马能谱和自然伽马能谱和温度、压力、套管接箍等辅助资料，获取的信息量大。

附图说明

图1：对两个同类探测器在同一尺度上采集的两列数据，做多尺度分解、融合、重构的原理示意图。

图2：从核测井原始数据，通过做多尺度分解、融合、重构，最后对含油气层做出多尺度地质解释的工作流程图。

图3：脉冲中子测井热中子时间谱多尺度滤波处理流程图。

图4：脉冲中子测井热中子时间谱多尺度滤波处理效果对比图：图4(a)原始热中子道计数时间谱；图4(b)纵向滤波后的时间谱；图4(c)纵-横双向滤波后的时间谱。

图5：多尺度分析识别油水层原理图：图5(a)原始热中子计数率；图5(b)纵向滤波结果；图5(c)横向滤波结果。

图6：τ曲线多尺度分析识别油水层的实例。该井为中国东部某油田的一口井，用Morlet小波基对中子寿命测井τ曲线进行多尺度分析，做出时频色谱图和不同尺度上的小波系数曲线。测井曲线经小波变换后成为深度—尺度域二维空间的函数，并且在不同的尺度和深度域显示出不同的幅值特征和周期性。对于中子寿命测井曲线，含油饱和度越高热中子随时间的衰减越慢，中子寿命测井曲线幅值越大，小波系数曲线幅值响应也较大。因此，色谱图中与较大小波系数幅度值相对应的灰色区域显示为油层，尺度范围为7—10。选定第8尺度为识别油水层的最佳尺度，依据小波系数曲线的响应幅度和周期振荡性特征，可认定4个含油储层。

图7：寿命-电阻率融合曲线多尺度分析定性识别油水层的实例。含油饱和度高的地层中子寿命长电阻率高，两列物理参数数列有很好的相关性，故可加以融合。利用Morlet小波基对融合曲线进行多尺度分析(尺度变化范围为1—10)，得到小波时频色谱图和不同尺度上的小波系数曲线。色谱图中沿浅-灰-深色区域方向显示小波系数从低值变化到高值，深色区域与油层具有一定的对应性。最佳尺度为第5尺度，小波系数幅值大并且周期振荡与油层相对应。小波系数曲线在28号地层幅值较大且振荡剧烈，定性为油层，其余为油水同层。23号层下部小波系数响应不明显，而色谱图上有较大响应，故应结合小波系数和色谱图共同来识别油水层。与单条测井曲线相比，融合曲线与解释结论比较吻合，识别效果较好。

具体实施方式

本发明是一种核测井数据多尺度处理方法，包括从核测井原始数据，通过多尺度分解、融合、重构，最后对含油气层做出多尺度地质解释的全过程，如附图2所示，包括以下步骤：

1.读入核测井数据

可读入的核测井数据包括：

1.1 自然伽马测井：自然伽马总计数率、自然伽马能谱和铀、钍、钾含量；

1.2 散射伽马测井：散射伽马总计数率、散射伽马能谱、密度和岩性指数；

1.3 同位素中子源测井：中子计数率、俘获伽马计数率、俘获伽马能谱和孔隙度；

1.4 脉冲中子测井：热中子总计数率和热中子时间谱；脉冲中子测井俘获伽马总计数率、俘获伽马时间谱和能谱；脉冲中子测井非弹伽马总计数率和能谱；中子活化伽马计数率、伽马时间谱和能谱。

根据施工目的可读入上述某一种核测井数据，也可是几种测井数据的组合。

2.对读入的测井数据进行横向和纵向多尺度滤波

多尺度滤波的步骤包括：选取具有一定的紧支撑性、对称性和平滑性的正交小波基，利用正交小波具有很强的去数据相关性，使信号的能量在小波域内集中于有用信号。确定分解层次N，然后对数据进行N层小波分解。因为噪声具有一定的频带宽度，在较高分解层中噪声的小波系数可以忽略，因此选择合适的分解尺度既是保证噪声完全去除的前提，同时又减少了计算代价。实际应用中小波分解层次取3至5层即可。

2.1 对数据进行N层小波分解的具体步骤包括：

在某尺度i上，对给定的核测井信号序列 $x(i,k)\&Element;V_i\&Subset;l^2\left(Z\right)$ (k∈z)，通过一个脉冲响应为h(k)的低通滤波器，获得粗尺度上的平滑信号

x_V(i-1，k)∈V_i-1                     (1)

$x_V(i-1,l)=\underset{k}{\mathrm{\&Sigma;}}h(2l-k)x(i,k)---\left(2\right)$

信号x(i，k)在低通滤波器中丢失的“细节信号”可以由x(i，k)通过一个脉冲响应是g(k)的高通滤波器得到

x_D(i-1，k)∈D_i-1                          (3)

$x_D(i-1,l)=\underset{k}{\mathrm{\&Sigma;}}g(2l-k)x(i,k)---\left(4\right)$

下标D表示x(i，k)在细节信号空间D_i-1上的投影。

循环进行上述步骤，进行N层的多尺度分解，得到2N个不同的频带，其中包含N个高频信号和N个低频信号。从低频信号中获得测井曲线中包含的地层概貌信息，从高频信号中获得地层的细节。

对第1层到第N层的每一层高频系数选择一个阈值进行软阈值量化处理。软阈值处理将较小的小波系数置零但对较大的小波系数向零做了收缩，可去掉测井曲线中与地层信息无关的统计起伏和其它高频噪声，去噪后的数据能生成更为光滑的曲线。

2.2 时间谱多尺度滤波的具体实施方式：

脉冲中子-中子PNN测井仪的中子源脉冲宽度很小而脉冲间隔过长，导致热中子道计数统计起伏较大。为了抑制时间道计数的统计不精确性和其它高频噪声，对采集的时间谱数据进行多尺度纵向和横向滤波。

纵向滤波是对不同深度点同一时间记录道上的热中子计数进行滤波，如表1所示。横向滤波是对同一深度点的60个时间计录道上的热中子计数进行滤波，如表2所示。

时间谱多尺度滤波过程分两步进行，首先对原始热中子道计数进行纵向滤波，然后对纵向滤波后的重构热中子计数进行横向多尺度滤波处理，工作流程如附图3所示。

深度/m            时间道计数(1—60道)

         Ch1  Ch2  Ch3  Ch4  Ch5  Ch6  Ch7..........Ch56 Ch57 Ch58 Ch59 Ch60

2284.5    39   34   63   69   66   66   50...........0    0    0    0    0

2284.4    35   31   51   69   67   58   58...........0    0    0    0    0

2284.3    36   36   51   55   58   53   45...........0    0    0    0    0

2284.2    47   48   69   77   76   75   80...........0    0    0    0    0

....

2221.2    37   44   52   51   19   10   9........... 0    0    0    0    0

2221.1    43   56   58   57   39   9    11...........0    1    1    1    0

2221.0    44   50   52   53   28   16   9........... 1    1    0    1    0

2220.9    46   58   70   46   3716      5........... 1    0    0    0    0

                         表1  纵向滤波示意图

深度/m                  时间道计数(1—60道)

        Ch1  Ch2  Ch3  Ch4  Ch5  Ch6  Ch7..........Ch56 Ch57 Ch58 Ch59 Ch60

2284.5   39   34   63   69   66   66   50...........0    0    0    0    0

2284.4   35   31   51   69   67   58   58...........0    0    0    0    0

2284.3   36   36   51   55   58   53   45...........0    0    0    0    0

2284.2   47   48   69   77   76   75   80...........0    0    0    0    0

....

2221.2   37   44   52   51   19   10   9........... 0    0    0    0    0

2221.1   43   56   58   57   39   9    11...........0    1    1    1    0

2221.0   44   50   52   53   28   16   9........... 1    1    0    1    0

2220.9   46   58   70   46   37   16   5........... 1    0    0    0    0

                   表2横向滤波示意图

选取大港油田某井的PNN热中子时间谱进行多尺度滤波处理。选用的短源距探测器记录的时间谱道计数，采样时间从快中子束发射30μs后到1800μs。

(1)纵向多尺度滤波

分别用db4，coif4，sym4，bior2.2四种小波基对纵向热中子计数进行处理后发现与sym4小波相对应的误差最小，故选sym4小波为最优小波基。此外，希望纵向滤波后的计数率曲线仍能够较好地反映原始地层的信息，保留地层中岩性发生突变的边缘特征，所以选用硬阈值处理方法。

(2)横向多尺度滤波

由于横向热中子随着时间的分布呈指数衰减规律，故选bior2.2小波为最优小波基。此外，希望横向滤波后的计数率曲线能比较光滑，有利于地层参数的求解，所以选用软阈值处理方法。

PNN测井热中子时间谱多尺度滤波处理前后的数据如附图4所示。图4(a)为三个深度点的原始热中子计数率，图4(b)为纵向滤波的结果，图4(c)为横向滤波的结果。从图中可看出，多尺度滤波后的时间谱数据变化比较光滑，对核测井信号的统计涨落具有很好的抑制作用，满足对时间谱后处理的需要。因此，多尺度滤波过程是对PNN测井数据处理的重要基础。

3.归一化处理

不同的测井数据有不同的量纲，因此要对滤波后的数据进行归一化处理，步骤包括：

$f\&prime;=\frac{f-f_{\min }}{f_{\max }-f_{\min }}\&times;100---\left(11\right)$

式中f_min为测井曲线的最小值，f_max为测井曲线的最大值，f为测井曲线的采样点实测数据。

4.多尺度分解

对归一化处理后的数据进行多尺度分解：横向时间轴多尺度分解的方法包括多指数拟合法和小波分解法；横向能量轴多尺度分解的方法包括解谱法和小波分解法；纵向深度轴多尺度分解的方法取小波分解法。

5.数据融合和重构

对测井数据进行包括在原始数据、物理参数和地质参数三个层次上分别进行时间谱或能谱融合和重构的步骤：

5.1 首先以小波模极大重构算法为基础，选用低频系数加权、高频系数取绝对值较大的融合规则，对同类或多类测井数据进行小波融合得到融合小波金字塔。然后对融合后的小波金字塔进行小波逆变换实现各个层次的数据重构。最后采用熵、均值和方差作为评判标准，对融合数据进行定量评价：均值、方差数据反映峰值信息；信息熵数据反映空间细节信息；峰值信噪比、相关系数数据反映频谱信息。

5.2 小波模极大重构算法的步骤包括：

当用测井曲线多尺度分析方法识别油水层时，设测井仪器接收到的信号为f(x)，Morlet连续小波变换定义为：

$W_f(a,b)=<f,{\mathrm{\&psi;}}_{\mathrm{ab}}>=\frac{1}{\sqrt{\left|a\right|}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}f\left(x\right)\mathrm{\&psi;}\left(\frac{x-b}{a}\right)\mathrm{dx}---\left(12\right)$

式中：Morlet小波基函数为

$\mathrm{\&psi;}\left(x\right)=e^{-x^2/2}\&CenterDot;e^{i{\mathrm{\&omega;}}_{0}x}---\left(13\right)$

其中，a和b分别为尺度因子和位移因子；改变a和b值，小波变换对不同频率在空间域上的取样步长也随之改变。

在低频时小波变换的空间分辨率较差，而频率分辨率较高；在高频时小波变换的空间分辨率较高，而频率分辨率较低。这符合低频信号变化缓慢而高频信号变化迅速的特点。

对信号进行小波分析，可以得到不同尺度a的评估信号在不同空间段的大量系数，即相应的小波系数W_f(a，b)，在色谱图中用颜色的变化来表示系数的大小。

某一尺度a时的小波系数曲线，能够直观显示小波系数与被分析信号之间的相似程度，且可以通过控制a的大小，使小波系数曲线和某些参数具有一定的相关性。

以小波模极大重构算法为基础，选用低频系数加权、高频系数取绝对值较大的融合规则，实现同类或多类测井信号数据融合，并采用熵、均值和方差作为评判标准，对融合数据进行定量评价。

按照融合规则经过小波重构就可获得融合数据。如脉冲中子双谱饱和度测井仪可在不同尺度上测得热中子和俘获伽马两个时间谱，分别记做数据D1和数据D2。对这两列数据分别进行分解，即对原始数据进行低、高通滤波，使每个原始数据都分解为含有不同频率成分的两个子数据，记为a1、d1和a2、d2。a1和a2是低频逼近信号，而d1和d2为高频细节信号。再根据需要，对低频子数据重复上面的过程，即可建立各数据的小波塔形分解。然后对各个分解层进行融合处理，最终得到融合小波金字塔。

数据融合效果可采用三类统计参数来进行分析与评价：第一类反映峰值信息，如均值、方差；第二类反映空间细节信息，如信息熵；第三类反映频谱信息，如峰值信噪比、相关系数。信噪比越大，说明融合后的数据精确度越高。融合曲线突出了测井曲线的公共信息，增加了地层信息量，提高了测井资料的可信度。从小波基的基本性质出发，选取几种常用小波基函数进行分析对比，以原始信号和重构信号的相关性作为小波基的选取依据，通过比较相关系数的大小，选定适合核测井多尺度分析方法的最优小波基，构建时频色谱图定出与油水层信息所对应的小波尺度范围。选取反映储层信息的最佳尺度，依据最佳尺度下小波系数曲线的幅值和小波系数振荡趋势区分油水层。

对信号进行小波分析，可得到不同尺度a上评估信号在不同空间段的大量系数，并且可以用色谱图显示出来。同一尺度上相同的颜色对应着物性相同或相近的储层，可用于识别油、气、水层。

由深度-寿命域τ分布图可求出热中子寿命加权平均值、热中子宏观俘获截面和含水饱和度。

6.显示曲线和图像

的曲线和图像包括：原始数据曲线，包括滤波后的总计数率和近、远探测器计数率比值曲线；经标准井检验选定的，多尺度分解得到的对含油饱和度或岩性分辨能力好的分量曲线；经标准井检验选定的，多尺度融合得到的对含油饱和度或岩性分辨能力好的融合曲线；物理和地质参数随深度的变化曲线，包括热中子寿命、宏观俘获截面、碳氧比、含水饱和度；深度-时间、深度-中子寿命或深度-伽马能谱二维图像；通过融合数据计算出地层的加权热中子寿命τ和宏观俘获截面∑值，并做出中子计数率、热中子寿命或宏观截面色谱图，直观反应储层岩性和油气水饱和度的分布状态。

7.进行地质解释

地质解释包括划分地质剖面和识别评价油气水层：

7.1 划分地质剖面了解含油气层段的地质背景：

用自然伽马和自然电位融合曲线判定岩性，对各种尺度上的低频和高频信号进行观察、分析和比较，可确定各个尺度的低频和高频信号与地层的空间和物理参数之间的相关性。低频信号能够较好地反映测井曲线中包含的地层概貌信息，而高频信号能反映地层的细节，依此可对地层不同级别的沉积旋回、层序界面和沉积相做定性和定量研究，用近、远探测器计数率比值曲线指示孔隙度的变化。在查明油气层所处的地质环境的基础上，根据粗尺度信号的变化划分厚层和层组，根据细尺度信号的变化对厚层和层组进行细分，给出岩石物性参数空间分布的细节特征。

7.2.识别和评价油气水层

识别油水层的方法：

7.2.1 用深度-时间、深度-中子寿命二维图像，判定井眼流体和近井区地层中子寿命的径向分布，经统计分析和标准井刻度建立不同岩性油、气、水层的识别模式，作为资料解释的基础。根据热中子总计数率与俘获伽马总计数率的差别和热中子寿命或宏观截面色谱图，定性、定量识别和评价油气水层。用热中子宏观截面∑或碳氧比可计算油气水饱和度，并可获得储层在开发过程中剩余油饱和度分布的变化。

7.2.2 中子寿命测井数据多尺度分析识别油气水层

7.2.2.1 中子寿命测井数据多尺度分析识别油气水层的原理：

中子寿命分布本身就是尺度的连续分解，时间道的分段处理就是对不同尺度范围内的信息重构。选用合适的小波基，如Morlet小波，提取单一参数曲线(如τ曲线)以及融合曲线的时频和幅值信息可用以分辨油水层。

Morlet小波是一种连续小波，它克服了离散小波在时间和尺度空间上对信号离散过于粗糙的缺点，可以反映信号中多种连续变化的周期信息。由于测井曲线中隐含的油水层信息的幅值和时频差异较小，故宜采用Morlet连续小波分析。

设测井仪器接收到的信号为f(x)，Morlet连续小波变换定义为：

$W_f(a,b)=<f,{\mathrm{\&psi;}}_{\mathrm{ab}}>=\frac{1}{\sqrt{\left|a\right|}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}f\left(x\right)\mathrm{\&psi;}\left(\frac{x-b}{a}\right)\mathrm{dx}---\left(14\right)$

式中：Morlet小波基函数为

$\mathrm{\&psi;}\left(x\right)=e^{-x^2/2}\&CenterDot;e^{i{\mathrm{\&omega;}}_{0}x},---\left(15\right)$

a和b分别为尺度因子和位移因子。

由于a和b值的变化，使小波变换对不同频率在空间域上的取样步长是调节性的，即在低频时小波变换的空间分辨率较差，而频率分辨率较高；在高频时小波变换的空间分辨率较高，而频率分辨率较低，这正符合低频信号变化缓慢而高频信号变化迅速的特点。

对信号进行小波分析后，可以得到不同尺度a的评估信号在不同空间段的大量系数，即相应的小波系数W_f(a，b)。并且可以用色谱图的方式把它表示出来，在色谱图中，系数的大小是以颜色的变化来表示的，如附图5b。

根据需要可以得到某一尺度a时的小波系数曲线(如图5c，d)。这些曲线能够直观地显示小波系数与被分析信号之间的相似程度，并可通过控制a的大小使小波系数曲线和某些参数具有一定的相关性。若它们之间有很好的相关性，就可以根据小波系数曲线快速、准确地判断相关参数的变化情况，可用以储层划分和油水层的识别。

小波变换不仅能反映信号的空间分辨能力，而且还能反映信号的幅值信息。附图5(a)表示在[0，2000]区间上分布的信号序列，信号1为cos(x)，分布区间[0，700)；信号2为cos(4x)，分布区间[700，1300)；信号3为4cos(x)，分布区间为[1300，2000]。其中信号1和2的幅值相同，频率不同；信号1和3的频率相同，幅值不同。图(b)是对(a)中序列选用Morlet小波多尺度分析后得到的时频色谱图，横坐标同图(a)，纵轴表示尺度值。颜色从黄到红再到蓝，表示小波系数从低值到高值变化。图中直观显示该序列在尺度域中存在两个周期，信号1和3对应的尺度相同，但由于幅值不同，颜色也不同。信号3幅值大，在色谱图上响应为蓝色；信号1的幅值较小，与红色相对应。图(c)和(d)分别表示尺度值为50和12的小波系数曲线。尺度值选取不同，对信号中频率的响应也不同。

中子寿命测井曲线幅值较大处与油层信息对应，其小波系数曲线幅值响应也较大；时频色谱图中，同一尺度上相同的颜色代表物性相同或相近的储层。因此，可借助于不同尺度下的小波系数曲线，反映测井信号的油水层信息，为利用多尺度分析方法定性识别油水层提供了理论依据。

7.2.2.2 中子寿命(τ)曲线多尺度分析定性识别油水层的实例：

选取大港—Y井3050-3120m深度段为研究对象。以该井热中子寿命(τ)曲线为例，定性识别油水层。

利用Morlet小波基对中子寿命测井曲线进行多尺度分析，做出时频色谱图和不同尺度上的小波系数曲线，如附图6所示。时频色谱图的尺度范围为1-10，从左到右依次增加，颜色由黄色到红色代表小波系数从低值到高值。从图中直观看出，测井曲线经小波变换后已成为深度—尺度域二维空间的函数，并且在不同的尺度和深度域显示出不同的幅值特征和周期性。对于中子寿命测井曲线，含油饱和度越高，热中子随时间的衰减越慢，中子寿命测井曲线幅值越大，小波系数曲线幅值响应也较大。因此，色谱图中与较大小波系数幅度值相对应的红色区域与油层相对应，对应的尺度范围为7-10。

选定第8尺度为识别油水层的最佳尺度，依据小波系数曲线的响应幅度和周期振荡性特征，可定性识别4个含油储层。

由图中可看出，不同的测井曲线作多尺度分解后，其识别油水层界面的效果也不相同。结合现场解释结论可知，中子寿命测井曲线识别第23和28储层的效果较好，而对识别25层效果不明显，应结合不同测井曲线定性识别油水层。

7.2.2.3 电阻率-中子寿命多尺度数据融合识别油水层：

采用基于小波模极大的数据融合方法，融合规则如下：

高频系数选取绝对值较大的小波系数，算法为

$d_n^{\left(j\right)}=\left\{\begin{array}{cc}{d}_{1n}^{\left(j\right)},& \mathrm{if}(\mathrm{abs}\left(d_{1n}^{\left(j\right)}\right)>\mathrm{abs}\left(d_{2n}^{\left(j\right)}\right))\\ d_{2n}^{\left(j\right)},& \mathrm{else}\end{array}\right.---\left(15\right)$

低频系数采用加权的方案，算法为： $a_n^{\left(j\right)}=\mathrm{\&alpha;}\&CenterDot;a_{1n}^{\left(j\right)}+\mathrm{\&beta;}\&CenterDot;a_{2n}^{\left(j\right)}$

其中，α和β是两组测井数据经主成分分析后得到的主成分值。

对来自不同探测器的测井数据进行融合，使之优势互补，以减小分辨油水层的多解性。在常规测井曲线中，电阻率曲线是划分油水层的重要依据，因此选取寿命曲线与电阻率曲线进行融合，并对融合后的曲线进行多尺度分析，以提高识别油水层的可信度。

选取大港—Y井3050-3120m深度段为研究对象，利用上述融合规则对电阻率Rt和中子寿命τ曲线做多尺度数据融合，得到融合曲线如附图6所示，而评价标准如表3。从表3中可以看出，融合曲线不仅峰值信息得到加强，而且包含了更多的空间细节信息。由于融合曲线与τ曲线具有很好的相关性，因此融合曲线的高值也与油层相对应。

表3 τ、RT融合数据评价标准

利用Morlet小波基对融合曲线进行多尺度分析(尺度变化范围为1—10)，得到小波时频色谱图和不同尺度上的小波系数曲线，如附图7所示。色谱图中黄-红-蓝代表小波系数值从低值变化到高值，蓝色区域与油层具有一定的对应性。最佳尺度为第5尺度，小波系数幅值大并且周期振荡与油层相对应。小波系数曲线在28层幅值较大且振荡剧烈，定性识别为油层，其余为油水同层。23层下部小波系数响应不明显，而色谱图上有较大响应，故应结合小波系数和色谱图共同来定性识别油水层。与单条测井曲线相比，融合曲线与解释结论比较吻合，识别效果较好。

本发明基本原理：

本发明是一种以多尺度估计理论为基础的核测井数据处理和分析方法，通过研究核测井信号的多尺度特性和不同尺度上的信息融合，从原始测井信号中提取出地层的本质特征，建立测井数据向地质目标的映射关系，提高测井资料的可信度，以期将测井资料更好地用于地层评价。

核测井中采集的随深度变化的测井曲线，如自然伽马、脉冲中子测井时窗、能窗或总计数率及各种单参数曲线，都是深度域的多尺度信息组合。本发明以多尺度估计理论为基础，建立了核测井多尺度系统，给出系统的模型参数，采用相关系数法选取最优小波基，实现核测井信号的系统描述和测井信号的多尺度分析，将测井数据从一维深度域拓展为二维深度—尺度域。依据多尺度分解和重构算法，采用不同的融合规则，实现同类和多类测井数据融合。

脉冲中子-中子PNN测井的最终目的是定量确定地层的剩余油饱和度，因此得到比较可靠的地层热中子寿命值是测井解释中的重要环节。本发明将热中子寿命看成是反映井眼-地层介质系统物理性质的一种尺度，针对PNN测井仪器短源距采集到的不同时间道热中子计数率，利用多尺度分析方法对时间谱数据进行了处理及应用研究。

在对原始时间谱数据进行纵、横向多尺度滤波的基础上，利用颜色代表热中子计数率，做出能够反映不同深度和不同时间中子计数率变化的色谱图。色谱图直观反映了热中子随时间的衰减快慢，可定性判断不同材料套管井段，识别储集层。此外，通过色谱图直观反映的热中子计数率在时间道上的横向延伸，可用来判断能准确反映地层真实信息的起止时间道。

选用Morlet小波基对中子寿命曲线以及融合曲线进行多尺度分析，依据小波时频色谱信息和小波系数幅度振荡特征可定性识别油水层。

基于热中子计数率随时间呈指数规律衰减，实现了热中子寿命(τ)谱的多指数反演算法。首先对蒙特卡罗模拟数据的热中子寿命谱进行了研究，分析了地层水矿化度，孔隙度和饱和度对地层τ本征值的影响。此外研究了井眼、时间道的选取对τ谱呈现双峰和单峰的影响。τ谱用于实际测井资料处理时可定性分辨岩性与油水层，由此可能引发寿命测井解释方法的改进和市场价值的提高。

Claims (1)

1、一种核测井数据多尺度处理方法，包括以下步骤：

1]读入测井数据：

所述的测井数据包括以下数据中的至少一组：

①自然伽马测井：自然伽马总计数率、自然伽马能谱和铀、钍、钾含量；

②散射伽马测井：散射伽马总计数率、散射伽马能谱、密度和岩性指数；

③同位素中子源测井：中子计数率、俘获伽马计数率、俘获伽马能谱和孔隙度；

④脉冲中子测井：热中子总计数率和热中子时间谱；脉冲中子测井俘获伽马总计数率、俘获伽马时间谱和能谱；脉冲中子测井非弹伽马总计数率和能谱；中子活化伽马计数率、伽马时间谱和能谱；

2]对测井数据进行滤波：

对读入的测井数据进行深度域的纵向滤波和沿与深度轴垂直的时间轴或能量轴的横向滤波；所述纵向和横向滤波的方法包括Kalman滤波、多点拟合滤波、多点平滑滤波或多尺度滤波；

3]对滤波后的数据进行归一化处理：

所述归一化处理是指计算每个采样点的归一化测井值：

$f\&prime;=\frac{f-f_{\min }}{f_{\max }-f_{\min }}\&times;100$

式中，f_min为测井曲线的最小值；f_max为测井曲线的最大值；f为测井曲线在每个采样点的实测数据；

4]对经过滤波和归一化处理后测井数据进行多尺度分解、融合和重构：

所述对测井数据进行融合和重构包括在数据、物理参数和地质参数三个层次上分别对单条测井曲线和多条测井曲线进行多尺度分解、融合和重构；

所述对单条测井曲线进行多尺度分解、融合和重构的步骤包括：

4.1.1]选择合适的多尺度分解法对测井曲线进行多尺度分解，所述横向时间轴上对时间谱多尺度分解的方法包括多指数拟合法和小波分解法；所述横向能量轴上对能量谱多尺度分解的方法包括解谱法和小波分解法；所述纵向深度轴上对总计数率、时窗计数率、能窗计数率、计数率比值、物理参数和地质参数测井曲线的多尺度分解的方法是指小波分解法；所述单条测井曲线包括：总计数率、时窗计数率、能窗计数率、物理参数或地质参数；

4.1.2]作出小波时频色谱图，通过颜色的变化显示出不同深度和不同尺度下小波系数的大小，进而使内部的能量聚集与分布得以清晰展示；

4.1.3]以小波模极大重构算法为基础，选用低频系数加权、高频系数取绝对值较大的融合规则，在各个层次上对测井数据进行小波重构得到融合小波金字塔；

4.1.4]对融合后的小波金字塔进行小波逆变换实现各个层次的数据重构；

4.1.5]采用熵、均值和方差作为评判标准，对融合数据进行定量评价：均值、方差数据反映峰值信息；信息熵数据反映空间细节信息；峰值信噪比、相关系数数据反映频谱信息；

所述对多条测井曲线的数据进行多尺度分解、融合和重构的步骤包括：

4.2.1]对多条原始测井曲线分别进行多尺度分解，即对原始数据分别进行低、高通滤波，使多列原始数据各自分解为含有不同频率成分的低频逼近信号和高频细节信号；所述横向时间轴上对时间谱多尺度分解的方法包括多指数拟合法和小波分解法；所述横向能量轴上对能量谱多尺度分解的方法包括解谱法和小波分解法；所述纵向深度轴上对总计数率、时窗计数率、能窗计数率、计数率比值、物理参数和地质参数测井曲线的多尺度分解的方法是指小波分解法；所述多条测井曲线包括不同源距探测器的总计数率、时窗计数率、能窗计数率、物理参数和地质参数；

4.2.2]根据需要，对低频逼近信号重复上面的过程，建立各数据的小波塔形分解；

4.2.3]对各个分解层进行融合处理，最终得到融合小波金字塔；

4.2.4]对融合后的小波金字塔通过小波逆变换法进行数据重构，所得的重构数据即为该多列测井曲线经融合重构的数据；

5]显示曲线和图像：

所述的曲线和图像包括以下曲线和图像中的至少一种：

原始数据曲线，包括滤波后的总计数率、时窗计数率、能窗计数率、计数率比值；

经标准井检验选定的，多尺度分解得到的对含油饱和度或岩性分辨能力好的分量曲线；

经标准井检验选定的，多尺度融合得到的对含油饱和度或岩性分辨能力好的融合、重构曲线；物理和地质参数随深度变化的曲线，包括热中子寿命、宏观俘获截面、碳氧比、含水饱和度；深度-时间、深度-中子寿命或深度-伽马能谱二维图像；

6]进行地质解释：

所述进行地质解释的方法包括以下方法中的至少一种：

用自然伽马曲线和自然电位融合曲线判定岩性；

根据粗尺度信号的变化划分厚层和层组，根据细尺度信号的变化对厚层和层组进行细分；

用近、远探测器计数率比值曲线指示孔隙度的变化；

用深度-时间、深度-中子寿命二维图像，判定井眼流体和近井区地层中子寿命的径向分布，经统计分析和标准井刻度建立不同岩性油、气、水层的识别模式，作为资料解释的依据；

根据热中子总计数率与俘获伽马总计数率的差别和热中子寿命或宏观截面色谱图，定性、定量识别和评价油、气、水层；

用热中子宏观截面∑或碳氧比可计算油气水饱和度，并可获得储层在开发过程中剩余油饱和度分布的变化；

对各种尺度上的低频和高频信号进行观察、分析和比较，可确定各个尺度的低频和高频信号与地层的空间和物理参数之间的相关性，依此可对不同级别的沉积旋回、层序界面和沉积相做定性和定量研究，并可用于厚层细分和描述剩余油气饱和度空间分布的细节做精细描述；

所述步骤2]中的多尺度滤波的步骤包括：

2.1]选取具有一定的紧支撑性、对称性和平滑性的正交小波基；

2.2]选择一个小波并确定其分解层次N，然后对数据进行N层小波分解；

2.3]对第1层到第N层的每一层高频系数选择一个阈值进行软阈值量化处理；

2.4]根据小波分解的第N层低频系数和经过量化处理后的第1层到第N层的高频系数，进行测井数据的重构，重构后的数据已消除了统计不确定性和其它高频噪声；

所述对数据进行N层小波分解的步骤包括：

2.2.1]在某尺度i上，对给定的核测井信号序列 $x(i,k)\&Element;V_i\&Subset;l^2\left(Z\right)$ (k∈Z)，通过一个脉冲响应为h(k)的低通滤波器，获得粗尺度上(低频段)的平滑信号：

x_V(i-1，k)∈V_i-1

$x_V(i-1,l)=\underset{k}{\mathrm{\&Sigma;}}h(2l-k)x(i,k)$

2.2.2]信号x(i，k)在低通滤波器中丢失的“细节信号”由x(i，k)通过一个脉冲响应是g(k)的高通滤波器得到细尺度(高频段)的细节信号：

x_D(i-1，k)∈D_i-1

$x_D(i-1,l)=\underset{k}{\mathrm{\&Sigma;}}g(2l-k)x(i,k)$

下标D表示x(i，k)在细节信号空间D_i-1上的投影；

2.2.3]对平滑信号逐次进行N层的多尺度分解，得到2N个不同的频带，其中包含N个高频信号和N个低频信号，从低频信号中获得测井曲线中包含的地层概貌信息，从高频信号中获得地层的细节；

所述步骤4]中的小波模极大重构算法的步骤包括：

设测井仪器接采集到的信号为f(x)，Morlet连续小波变换定义为：

$W_f(a,b)=<f,{\mathrm{\&psi;}}_{\mathrm{ab}}>=\frac{1}{\sqrt{\left|a\right|}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}f\left(x\right)\mathrm{\&psi;}\left(\frac{x-b}{a}\right)\mathrm{dx}$

式中：Morlet小波基函数为

$\mathrm{\&psi;}\left(x\right)=e^{-x^2/2}\&CenterDot;e^{i{\mathrm{\&omega;}}_{0}x}$

a和b分别为尺度因子和位移因子；

对信号进行小波分析，可以得到不同尺度a的评估信号在不同空间段的大量系数，即相应的小波系数W_f(a，b)，并用色谱图的方式表示；在色谱图中，用颜色的变化来表示系数的大小；可选择得到某一尺度a时的小波系数曲线，这条曲线能够直观显示小波系数与被分析信号之间的相似程度，且可通过控制a的大小，使小波系数曲线和某些参数具有一定的相关性。

Images