CN116520419A

CN116520419A - 一种热流体裂缝通道识别方法

Info

Publication number: CN116520419A
Application number: CN202310491495.1A
Authority: CN
Inventors: 许超; 杨瑞召; 陈井瑞; 韩枫涛; 冯洋洋
Original assignee: China University of Mining and Technology Beijing CUMTB
Current assignee: China University of Mining and Technology Beijing CUMTB
Priority date: 2023-05-05
Filing date: 2023-05-05
Publication date: 2023-08-01
Anticipated expiration: 2043-05-05
Also published as: CN116520419B

Abstract

本发明公开了一种热流体裂缝通道识别方法，该方法包括：获取地表采集的天然地震数据；对天然地震数据进行滤波处理，通过滤波方法去除近地表中存在的相干噪音；对地下进行网格化，将每个网格点视为虚拟源，根据波动方程或射线追踪方法计算各个网格点到地表台站的旅行时；将连续记录的天然地震资料分成若干个时段，利用时频谱筛选出有效信号时间段，采用滑动时窗非线性叠加的方法计算各个网格点叠加后的相似系数；给定门槛值对三维网格点所有的相似系数进行阈值化处理，其中相似系数较大的代表着裂缝越发育。可以最大程度的提高被动地震中弱信号识别的信噪比，从而增强中深层地热流体裂缝通道的刻画能力。

Description

一种热流体裂缝通道识别方法

技术领域

本发明属于热流体通道识别领域，具体涉及一种热流体裂缝通道识别方法。

背景技术

热流体通道对于裂隙热储来说，是指地表水沿构造裂隙带下渗的深循环迳流通道，以及成为热流体后又沿断裂带上升至地表的通道；对于孔隙热储来说，是指穿透隔水阻热盖层的断裂构造，即被封闭的古沉积水上涌的通道。不管是哪种热流体通道，其都与深大断裂以及地质体存在的天然裂缝相关。裂缝中存在两种不同类型的力学作用，一方面，在粘性岩石中，裂缝在其尖端产生破裂和剪切时会产生短时间的声发射，一般这些信号的工业应用被称为微地震。另一方面，在地球岩石圈变形过程中会发生岩石-裂缝-流体耦合作用，这些会产生特殊的裂缝物理特征信号。这种信号表现为充满流体的裂缝的谐波共振。这种信号会随着局部应力的外部变化和内部变化引起，例如经过的地震波、构造变形、水力压裂等。

在当前背景下，地热勘探的区域通常是无法采用三维地震勘探的，一方面是地质条件复杂，另一方面是勘探需要的成本太高，而圈定热流体通道的规模以及空间形态分布特征对于后期井位的部署非常重要。

因此，发展一种基于被动地震手段识别中深层地热流体裂缝通道的方法具有重要意义。

发明内容

本发明所为了解决背景技术中存在的技术问题，基于被动地震资料，构建了一种利用信噪比极低的微破裂信号对地热田中的热流体通道进行识别和刻画的方法；包括预处理滤波、时频谱有效信号段筛选、建立速度模型计算走时信息、利用走时信息拉平信号并叠加和阈值化提取裂缝通道。

为了解决技术问题，本发明的技术方案是：

一种热流体裂缝通道识别方法，所述方法包括：

获取地表采集的天然地震数据，进行滤波处理，去除近地表中存在的相干噪音；

对观测目标区地下空间进行网格化，将每个网格点视为虚拟源，计算各个网格点到地表台站的走时并记录；

将连续记录的滤波后天然地震资料分成多个时段，利用时频分析自动筛选出含有有效信号的时段，针对每个时段采用滑动时窗的方法计算各个网格点叠加后的相似系数；

预设门槛值对三维网格点所有的相似系数进行阈值化处理，即得到地下裂缝三维空间形态分布特征，其中相似系数较大的代表着裂缝越发育。

进一步，所述滤波处理，具体包括：

采用倒谱滤波方法并在频率域上进行，先将时域信号通过傅里叶变换到频率域中，再利用一次傅里叶变换将频率信号转换到倒谱域当中，此时裂缝地震信号集中在倒谱域低频中，而噪音分布在全倒谱域当中，通过对倒谱域中的信号进行低通滤波，将有效信号分离出来，倒谱变换的具体如下：

设地震道是一个小波(w)与反射系数(r)在时间上的褶积，则地震信号(s)可表述为：

s(t)＝w(t)*r(t)

在频率域中：

S(f)＝W(f)×R(f)

倒谱是振幅谱对数的傅里叶变换：

ln[S]＝ln[W]+ln[R]

F[ln[S]]＝F[ln[W]]+F[ln[R]]

倒谱变换的特性通过使用低通滤波将小波从反射系数中分离出来，其中低频率包含频谱包络的平滑形状，而其余的频率包含反射系数的响应；地震记录中的单色信号变成了频率上的单一峰值，这些尖峰不会被倒谱滤波器去除，而是在倒谱域带通滤波后分布在一个较窄的频带内；因此，需要在倒谱滤波后，频率域信号恢复到时域信号之前，使用中值滤波器对滤波后的振幅谱进行平滑处理，以去除地震记录中无法被倒谱滤波去除的单色噪音。

进一步，选用最感兴趣目标点的地表检波器对去噪后的信号进行时频分析，通过时频谱自动筛选含有有效信号的时间段参与后续处理。

进一步，利用射线追踪计算走时，计算3D射线追踪的方法选用快速行进法FMM、最短路径法SPM和快速扫描法FSM中的其中之一。

进一步，所述计算各个网格点到地表台站的走时，具体包括：

采用快速扫描法计算裂缝地震法中单向深度偏移需要的走时信息，三维空间中的程函方程如下：

公式中的偏导数采用迎风差有限差分格式可表示为：

其中u_i,j,k表示为节点(i,j,k)到时值，Δx、Δy、Δz为离散网格边长；令空间步长h＝Δx＝Δy＝Δz，则有：

其中：

对于模型边界则采用单侧网格差分格式进行计算，如对于模型边界面x＝1上的网格节点值其差分格式为：

具体算法更新策略如下：

对某一网格点走时更新时，将其对应的从小到大排列为W、V、U，迭代后节点到时更新值为μ^*，令μ^*＝W+r(x,y,z)；考虑等网格划分模型h＝Δx＝Δy＝Δz；

若μ^*<V，那么μ^*就是该网格节点处的解；否则：

然后继续下一步判断；

若μ^*<U，那么μ^*就是该网格节点处的解，否则：

根据上述步骤得到所有更新节点处到时的唯一解；

其中，在地表存在起伏时，规则网格计算走时信息会产生一定的误差，近地表的影响通过高程静校正和剩余静校正减弱；于高程静校正，提取检波器与基准面之间的高程差，并利用近地表速度计算高程差的传播时间，在进行深度偏移之前，利用它来校正地震记录；

分析高程静校正、剩余静校正和速度模型校正的最佳方法是利用检波器上可以观察到的射孔炮的波形，结合射孔炮的已知位置信息对地震记录进行反演约束。

进一步，采用滑动时窗的方法计算各个网格点叠加后的相似系数，具体包括：

利用计算出的走时信息对地下每个虚拟源点即网格剖分节点进行单向旅行时偏移堆叠成像；裂缝层析成像将计算滑动时窗内地下每个网格点的Semblance系数，针对各个网格点的相似性度量在考虑了有效信号幅值的同时也考虑了相位信息，每个网格点的Semblance系数计算公式优化如下：

其中D代表Semblance相干系数，i代表地震道号，j代表滑动时窗内的时间点号，K代表所有需要被堆叠的地震道数，s_ij代表地震道i在时间j处的地震振幅；其中S_i(t)＝s_i(t)+iH(s_i(t))，S_i(t)是原始信号s_i(t)构建的复信号；

相似系数(D)是在一个滑动时窗内每个体单元网格点声发射相关性的度量；裂缝地震法将对地震道上的连续记录进行滑动时窗测量；每一个时间单元的能量扫描叠加可以获得一个地下空间能量扫描体，连续记录的滑动时窗测量将获得成千上万的地下空间能量扫描体；通过剔除那些由较强能量占主导的能量体，堆叠其它质量较高的能量体可以获得地下裂缝能量体的空间分布特征；其中能量较高的部分代表裂缝发育概率较大，裂缝发育较密集的区域。

进一步，进行阈值化处理得到的能量扫描叠加体，考虑台站硬件的物理内存限制导致连续的观测资料必定是分段的，并且时间的分段也可以是人为定义；在将不同的时间段进行叠加前，各个时间段叠加求和得到的能量扫描叠加体即以各个网格点的相似系数度量的需进行归一化处理，最后将所有有效时段的能量扫描体进行求和，根据实际情况选择一定的门槛值，保留大于门槛值的部分即是最终的热流体裂缝通道。

与现有技术相比，本发明的优点在于：

一种应用于地热领域热流体裂缝通道识别的方法，该方法包括：获取地表采集的天然地震数据；对天然地震数据进行滤波处理，通过滤波方法去除近地表中存在的相干噪音；利用时频谱自动筛选具有有效信号的时间段；对地下进行网格化，将每个网格点视为虚拟源，根据波动方程或射线追踪方法计算各个网格点到地表台站的旅行时；将连续记录的天然地震资料分成若干个时段，采用滑动时窗的方法计算各个网格点叠加后的相似系数；给定门槛值对三维网格点所有的相似系数进行阈值化处理，其中相似系数较大的代表着裂缝越发育。可以最大程度的提高被动地震中弱信号识别的信噪比，从而增强中深层地热流体裂缝通道的刻画能力，缺少其中任何一个处理流程都有可能导致裂缝层析成像的失败。

附图说明

图1、本发明一种热流体裂缝通道识别方法的主流程图；

图2、倒谱滤波预处理前后对比示意图

图3、基于时频谱自动化识别的有效信号段。

图4、裂缝能量扫描体示意图。

图5、Semblance算法优化后相似系数计算结果对比图。

图6、能量扫描叠加体的分层切片示意图。

图7、阈值化处理的裂缝体三维雕刻示意图。

具体实施方式

下面结合实施例描述本发明具体实施方式：

需要说明的是，本说明书所示意的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容能涵盖的范围内。

同时，本说明书中所引用的如“上”、“下”、“左”、“右”、“中间”及“一”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

实施例1：

本发明基于被动地震资料，如图1所示，构建了一种利用信噪比极低的微破裂信号对地热田中的热流体通道进行识别和刻画的方法。技术方案的处理流程包括预处理滤波、有效信号段筛选、建立速度模型计算走时信息、利用走时信息拉平信号并叠加、阈值化提取裂缝通道。

第一步，必须要对近地表采集到的天然地震资料进行预处理滤波。本技术方案的重要预处理步骤是需要对信号进行高质量的非共振信号分析和抑制。本发明中使用的是倒谱滤波方法(如图2所示)，利用它在不影响从深部发射并传播到地面检波器接受到的信号波形相位特征的基础上，抑制反射系数噪声和表面波噪声。倒谱滤波可以有效区分深层震源产生的子波波形和通过岩石传播产生的重复及畸变波。这里使用的倒谱滤波方法是在频率域上进行的。先将时域信号通过傅里叶变换到频率域当中，再利用一次傅里叶变换将频率信号转换到倒谱域当中，此时裂缝地震信号将主要集中在倒谱域低频当中，而噪音则广泛分布在全倒谱域当中，通过对倒谱域中的信号进行低通滤波，可以将有效信号分离出来。倒谱变换的描述如下：

假设地震道是一个小波(w)与反射系数(r)在时间上的褶积，则地震信号(s)可表述为：

s(t)＝w(r)*r(y)

在频率域中：

S(f)＝W(f)×R(f)

倒谱是振幅谱对数的傅里叶变换：

ln[S]＝ln[W]+ln[R]

F[ln[S]]＝F[ln[W]]+F[ln[R]]

倒谱变换的特性可以通过使用低通滤波将小波从反射系数中分离出来，其中低频率包含频谱包络的平滑形状，而其余的频率包含反射系数的响应。地震记录中的单色信号变成了频率上的单一峰值，这些尖峰不会被倒谱滤波器去除，而是在倒谱域带通滤波后分布在一个较窄的频带内。因此，通常需要在倒谱滤波后，频率域信号恢复到时域信号之前，使用中值滤波器对滤波后的振幅谱进行平滑处理，以去除地震记录中无法被倒谱滤波去除的单色噪音。

第二步，需要利用射线追踪计算走时。常用的计算3D射线追踪的方法有快速行进法(FMM)，最短路径法(SPM)，快速扫描法(FSM)。其中，快速扫描法适用范围更广，效率更高。这里使用快速扫描法计算裂缝地震法中单向深度偏移需要的走时信息。三维空间中的程函方程如下(Huang J W,2012)：

式(1)中的偏导数采用迎风差有限差分格式可表示为：

其中u_i,j,k表示为节点(i,j,k)到时值，Δx、Δy、Δz为离散网格边长。令空间步长h＝Δx＝Δy＝Δz，则有：

其中：

对于模型边界则采用单侧网格差分格式进行计算，例如对于模型边界面x＝1

上的网格节点值其差分格式为：

具体算法更新策略如下：

1、对某一网格点走时更新时，将其对应的从小到大排列为W、V、U，迭代后节点到时更新值为μ^*，令μ^*＝W+r(x,y,z)。考虑等网格划分模型h＝Δx＝Δy＝Δz；

2、若μ^*<V，那么μ^*就是该网格节点处的解。否则：

然后继续下一步判断。

3、若μ^*<U，那么μ^*就是该网格节点处的解。否则：

根据上述步骤可以得到所有更新节点处到时的唯一解。

可选的，在地表存在起伏时，规则网格计算走时信息会产生一定的误差，这种误差可以通过一定的技术手段减弱地形对地形成像结果的影响。近地表的影响可以通过高程静校正和剩余静校正减弱。于高程静校正，可以提取检波器与基准面之间的高程差，并利用近地表速度计算高程差的传播时间。在进行深度偏移之前，利用它来校正地震记录。

分析高程静校正、剩余静校正和速度模型校正的最佳方法是利用检波器上可以观察到的射孔炮的波形，结合射孔炮的已知位置信息对地震记录进行反演约束。这种分析方法通常在水力压裂工业中应用，在缺乏射孔资料的地区，通常只能进行高程校正处理。

第三步，选择最感兴趣目标点的地表检波器去噪后的记录进行时频分析，通过时频谱自动筛选含有有效信号的时间段参与后续处理，如图3所示为裂缝内Krauklis波引起的谐波共振(由裂缝尖端的错位和内部流体运移引发)信号。这些信号是刻画热流体裂缝通道有效信号的重要组成部分。

第四步，利用校正后的走时信息对地下每个虚拟源点(网格剖分节点)进行单向旅行时偏移堆叠成像(图4)。裂缝层析成像将计算滑动时窗内地下每个网格点的Semblance系数，本发明技术方案中针对各个网格点的相似性度量在考虑了有效信号幅值的同时也考虑了相位信息，每个网格点的Semblance系数计算公式优化如下：

其中D代表Semblance相干系数，i代表地震道号，j代表滑动时窗内的时间点号，K代表所有需要被堆叠的地震道数，s_ij代表地震道i在时间j处的地震振幅。其中S_i(t)＝s_i(t)+iH(s_i(t))，S_i(t)是原始信号s_i(t)构建的复信号。

相似系数(D)是在一个滑动时窗内每个体单元网格点声发射相关性的度量。裂缝地震法将对地震道上的连续记录进行滑动时窗测量。每一个时间单元的能量扫描叠加可以获得一个地下空间能量扫描体，连续记录的滑动时窗测量将获得成千上万的地下空间能量扫描体。通过剔除那些由较强能量占主导的能量体，堆叠其它质量较高的能量体可以获得地下裂缝能量体的空间分布特征。其中能量较高的部分代表裂缝发育概率较大，裂缝发育较密集的区域。如图5中a为旅行时校正后的5道地震记录，它们之间存在着随机剩余动校正和静校正，b为分段滑动时窗计算的各段时窗相似系数，其中虚线为传统裂缝地震成像中Semblance算法计算的相似系数在时间上的变化，而实线为本申请方法计算的相似系数在时间上的变化，可以显著的看到，改进后方法计算的相似系数在时间上的相对信噪比强于传统方法。

第五步，阈值化处理得到的能量扫描叠加体。考虑到台站硬件的物理内存限制导致连续的观测资料必定是分段的，并且时间的分段也可以是人为定义的。在将不同的时间段进行叠加前，各个时间段叠加求和得到的能量扫描叠加体(以各个网格点的相似系数度量的)需要进行归一化处理，最后将所有有效时段的能量扫描体进行求和，根据实际情况选择一定的门槛值，保留大于门槛值的部分即是最终的热流体裂缝通道。如图6所示为能量扫描叠加体的分层切片示意图，图7为阈值化处理的裂缝体三维雕刻示意图。

本发明的技术方案通过以上五步处理流程，可以最大程度的提高被动地震中弱信号识别的信噪比，从而增强中深层地热流体裂缝通道的刻画能力，缺少其中任何一个处理流程都有可能导致裂缝层析成像的失败。

上面对本发明优选实施方式作了详细说明，但是本发明不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

不脱离本发明的构思和范围可以做出许多其他改变和改型。应当理解，本发明不限于特定的实施方式，本发明的范围由所附权利要求限定。

Claims

1.一种热流体裂缝通道识别方法，其特征在于，所述方法包括：

2.根据权利要求1所述的一种热流体裂缝通道识别方法，其特征在于，所述滤波处理，具体包括：

s(t)＝w(t)*r(t)

在频率域中：

S(f)＝W(f)×R(f)

倒谱是振幅谱对数的傅里叶变换：

ln[S]＝ln[W]+ln[R]

F[ln[S]]＝F[ln[W]]+F[ln[R]]

3.根据权利要求2所述的一种热流体裂缝通道识别方法，其特征在于，所述计算各个网格点到地表台站的走时，具体包括：

公式中的偏导数采用迎风差有限差分格式可表示为：

其中u_i，j，k表示为节点(i，j，k)到时值，Δx、Δy、Δz为离散网格边长；令空间步长h＝Δx＝Δy＝Δz，则有：

其中：

具体算法更新策略如下：

对某一网格点走时更新时，将其对应的从小到大排列为W、V、U，迭代后节点到时更新值为μ^*，令μ^*＝W+r(x，y，z)；考虑等网格划分模型h＝Δx＝Δy＝Δz；

若μ^*<V，那么μ^*就是该网格节点处的解；否则：

然后继续下一步判断；

若μ^*<U，那么μ^*就是该网格节点处的解，否则：

根据上述步骤得到所有更新节点处到时的唯一解；

4.根据权利要求3所述的一种热流体裂缝通道识别方法，其特征在于，选用最感兴趣目标点的地表检波器对去噪后的信号进行时频分析，通过时频谱自动筛选含有有效信号的时间段参与后续处理。

5.根据权利要求1所述的一种热流体裂缝通道识别方法，其特征在于，采用滑动时窗的方法计算各个网格点叠加后的相似系数，具体包括：

6.根据权利要求1所述的一种热流体裂缝通道识别方法，其特征在于，进行上述阈值化处理得到的能量扫描叠加体，考虑台站硬件的物理内存限制导致连续的观测资料必定是分段的，并且时间的分段也可以是人为定义；在将不同的时间段进行叠加前，各个时间段叠加求和得到的能量扫描叠加体即以各个网格点的相似系数度量的需进行归一化处理，最后将所有有效时段的能量扫描体进行求和，根据实际情况选择一定的门槛值，保留大于门槛值的部分即是最终的热流体裂缝通道。