CN103513277B - 一种地震地层裂隙裂缝密度反演方法及系统 - Google Patents

Abstract

本发明提供了一种地震地层裂隙裂缝密度反演方法及系统，通过全波形反演充分利用了采集到的地震波信息，并自动考虑了地震波在地下的传播行为，能给出高精度地震反演结果。本发明通过HTI介质中弹性全波形反演理论，提出了全波形裂隙裂缝密度反演方法，实现了地震地下地层裂隙裂缝密度的定量预测。

Description

一种地震地层裂隙裂缝密度反演方法及系统

技术领域

本发明涉及地震勘探领域，特别涉及一种地震地层裂隙裂缝密度反演方法及系统。

背景技术

裂隙裂缝是油气的储集空间和渗滤通道，也是水、煤层气富集、储存、运移的场所。探明地下裂隙裂缝分布和发育程度对于油气预测、提高油气产能以及预防瓦斯爆炸、顶\底板突水和煤层气开发等都至关重要。

描述裂隙裂缝分布和发育程度的参数主要有裂隙裂缝的走向、裂隙裂缝的密度。目前地震法裂隙裂缝探测主要有横波分裂法和极化图法以及纵波方位AVOA法和方位旅行时法等。理论上，横波对地层中裂隙裂缝的存在比纵波更为敏感，横波极化图法和横波分裂法都可以很好地预测裂隙裂缝的走向；在地层厚度确定的情况下，横波分裂法中两个分裂横波旅行时的差反映了裂隙裂缝的发育程度，如时差大，裂隙裂缝发育程度高，地层中裂隙裂缝密度大。但目前采集的横波数据的质量整体都比较差，这就使测量的准确度大大降低，不能获得很好的横波裂缝预测结果。纵波地震勘探采集的数据质量整体来说都比较高，是目前地震探测的主要手段。纵波方位AVO法和方位旅行时(或方位层速度)法裂隙裂缝预测都采用椭圆拟合方法，椭圆的长轴指示裂隙裂缝的走向，椭圆的长轴和短轴之比反映裂隙裂缝发育程度即裂隙裂缝密度。但无论是分裂横波的旅行时差还是纵波椭圆长短轴的比都只能给出地层中裂隙裂缝密度的相对大小(即发育程度)，不能满足对地层渗透能力的定量评价和裂缝型油气藏油气储量计算要求。

发明内容

本发明的主要目的在于解决现有技术中存在的问题，提供地震地层裂隙裂缝密度反演方法及系统。

本发明的目的是通过下述技术方案予以实现的：

本发明一方面提供一种地震地层裂隙裂缝密度反演方法，包括以下步骤：

步骤一，采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

步骤二，在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并 计算地层的质量密度；

步骤三，在所述三分量炮集数据中径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度；

步骤四，根据所述深度域纵波速度、所述深度域横波速度以及所述地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

步骤五，利用所述拉梅常数、剪切模量和地层质量密度为初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

步骤六，通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙裂缝走向模型，其中，所述裂隙裂缝走向与正北方向的夹角记为θ；

步骤七，利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

步骤八，计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，当E<δ时，输出反演结果，并停止反演，否则，进入步骤九；δ取值范围为：1.0e^-4～1.0e^-6；

步骤九，将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；

步骤十，利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，计算反演公式中的系数张量C_pjkl(x),下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向；

步骤十一，将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向；T为地震记录长度，Ns为用于反演的炮的总数；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……， Ns；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

步骤十二，根据所述裂隙裂缝密度梯度，按以下公式迭代修改裂隙裂缝密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为所述第n次迭代修改后的裂隙裂缝密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙裂缝密度梯度；

步骤十三，将修改后的裂隙裂缝密度模型作为初始裂隙裂缝密度模型，即取裂隙裂缝密度e₀＝e_n+1,进入步骤七。

本发明另外提供一种地震地层裂隙裂缝密度反演系统，包括以下单元：

预处理单元，用于采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

深度域纵波速度及地层质量密度计算单元，用于在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并计算地层的质量密度；

深度域横波速度计算单元，用于在所述三分量炮集数据中径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度；

初始模型计算单元，用根据所述深度域纵波速度、所述深度域横波速度以及所述地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

模型计算单元，用于利用所述拉梅常数、剪切模量和地层质量密度为初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

裂隙走向模型计算单元，用于通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙裂缝走向模型，其中，所述裂隙裂缝走向与正北方向的夹角记为θ；

正演模拟单元，用于利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

反演结果输出单元，用于计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，当E<δ时，输出反演结果，并停止反演，否则，进入逆时传播单元；δ取值范围为：1.0e^-4～1.0e^-6；

逆时传播单元，用于将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；

系数张量计算单元，用于利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，计算反演公式中的系数张量C_pjkl(x)(其中下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向)；

裂隙裂缝密度梯度值计算单元，用于将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以步骤10获得的系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向；T为地震记录长度，Ns为用于反演的炮的总数；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……，Ns；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

迭代修改单元，用于根据所述裂隙裂缝密度梯度，按以下公式迭代修改裂隙裂缝密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为所述第n次迭代修改后的裂隙裂缝密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙裂缝密度梯度；

循环单元，用于将修改后的裂隙裂缝密度模型作为初始裂隙裂缝密度模型，即取裂隙裂缝密度e₀＝e_n+1，进入正演模拟单元。

本发明能够达到以下有益效果：全波形反演充分利用了采集到的地震波信息，并自动考虑了地震波在地下的传播行为，能给出高精度地震反演结果。本发明通过HTI介质中弹性全波形反演理论，提出了全波形裂隙裂缝密度反演方法，实现了地震地下地层裂隙裂缝密度的定量预测。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成 对本发明的限定。在附图中：

图1a为实施例一含断层的单层裂隙裂缝地质模型裂隙密度横向变化示意图；

图1b为实施例一含断层的单层裂隙裂缝地质模型地质剖面示意图；

图1c为实施例一含断层的单层裂隙裂缝地质模型层2裂隙示意图；

图2为实施例一近断层沿深度方向实际裂隙密度和反演裂隙裂缝密度图，其中实线表示实际裂隙密度，虚线表示初始模型，点划线表示反演裂隙密度；

图3为实施例一沿层实际裂隙密度和反演裂隙裂缝密度图，其中实线表示实际裂隙密度，虚线表示初始模型，点画线表示反演裂隙密度；

图4为实施例一地震地层裂隙裂缝密度反演方法流程图；

图5为单炮集实验数据采集图；

图6为实施例一中某炮数据的示意图；

图7为实施例二五层介质裂隙裂缝地质模型示意图；

图8为实施例二中某炮数据的示意图；

图9为实施例二沿深度方向实际裂隙裂缝密度和反演裂隙裂缝密度图；

图10为本发明一种地震地层裂隙裂缝密度反演系统的结构图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

针对前述现有技术中所存在的问题，本发明通过充分利用采集到的地震波信息，并同时考虑了地震波在地下的传播行为，通过HTI介质中弹性全波形反演理论，提出了全波形裂隙裂缝密度反演方法，实现了地震地下地层裂隙裂缝密度的定量预测。

以下结合具体实施例对本发明进行详细阐述：

实施例一：

本实施例为含有断层的三层地质模型，模型大小为1000m×800m；网格大小为：DX＝DZ＝10m，断层落差为10m，上下两层地层相同，中间层发育有裂隙，裂隙层厚度为30m，裂隙走向与x₁轴平行，模型参数如表1所示：

表1

相应的模型示意图如图1a、图1b以及图1c所示。裂缝密度岩深度方向变化如图2实线所示，断层附近裂隙发育，向两侧裂隙密度逐渐减小，如图3实线所示。

如图4所示，为本实施例一种地震地层裂隙裂缝密度反演方法流程图，包括以下步骤：

步骤401，采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

本实施例在上述模型中采集15炮三分量地震数据，单炮数据的采集图如图5所示；图6为该15炮地震数据中某炮的数据。其中，每炮101个检波点，检波点间距10m，最小炮检距0m，最大炮检距1000m，时间采样率1ms，采集长度500ms，对其进行预处理后，得到相应的三分量炮集数据。

预处理的过程包括两水平分量旋转，三分量静校正、地表一致性振幅补偿和叠前去除噪音等。

步骤402，在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并计算地层的质量密度；

在该15炮三分量炮集数据中的垂直分量数据中抽取共中心点道集数据，进行纵波速度分析；

计算时间域纵波速度：其中t_pi为第i层底纵波双程旅行时，V_pi为第i层的纵波速度，V_r,i为第i层纵波均方根速度；

需要说明的是，每一个层都有顶和底，因而第i层也会有顶和底。

进行时间深度转换：

其中t_pi为第i层底纵波双程旅行时，h_i为第i层的厚度，V_pi为第i层纵波速度；

本实施例中，地层质量密度用Gardner公式计算：

$\mathrm{\&rho;}=0.31V_p^{0.25};$

其中ρ为质量密度，V_p为纵波速度；

所述深度域纵波速度为所有地层纵波速度的集合：

V_p＝(V_p1,V_p2,…,V_pi,…,V_pN)_，其中，N为模型中地层的层数，V_pi为第i层纵波速度，

在本例中V_p＝(V_p1,V_p2,V_p3)。

步骤403，在所述三分量炮集数据中的径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度；

在15炮水平分量数据中抽取共转换点道集数据；

计算横波层速度：其中r_i为所述转换波速度分析获得的第i层纵横波速度比，V_pi为第i层纵波速度，V_si为第i层横波层速度；

进行时间-深度转换：

h_i＝(t_si-t_s(i-1))V_si；

其中，t_si为第i层底横波旅行时，h_i为第i层的厚度，V_si为第i层横波层速度；

需要说明的是，每一个层都有顶和底，因而第i层也会有顶和底。

所述第i层底横波旅行时根据以下公式计算：

$t_{si}=t_{ci}-\frac{1}{2}{t}_{pi};$

其中，t_ci为第i层转换波旅行时，t_pi为第i层纵波双程旅行时。

步骤404，根据所述深度域纵波速度、所述深度域横波速度以及所述地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

拉梅常数根据以下公式计算：

$\mathrm{\&lambda;}=\mathrm{\&rho;}(V_p^2-2V_s^2);$

其中λ为拉梅常数，ρ为质量密度，V_p为纵波速度，V_s为横波速度；

所述剪切模量根据以下公式计算：

μ＝ρV_s ²；

其中μ为剪切模量，ρ为质量密度，Vs为横波速度，为所有地层横波速度的集合，公式表述如下：

V_s＝(V_s1,V_s2,…,V_si,…,V_sN)；其中，N为模型中地层的层数，V_si为第i层横波速度，在本例中V_s＝(V_s1,V_s2,V_s3)。

步骤405，利用所述拉梅常数、剪切模量和地层质量密度为初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

该多分量数据处理方法即为多分量数据双反演方法。

步骤406，通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙走向模型，其中，所述裂隙走向与正北方向的夹角记为θ；

本实施例中，优先采用纵波方位AVO(AmplitudeVersusOffset,振幅随偏移距的变化)法获取裂隙裂缝走向。

步骤407，利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和正演模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

本实施例中，初始裂隙裂缝密度模型取值如图2中虚线所示；模拟地震波场：

地震波场正演模拟是利用HTI介质中弹性波动方程和伪谱法进行模拟；其中，震源函数采用零相位雷克子波。

步骤408，计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，输出反演结果或进入步骤409。

当E<δ时，输出反演结果，并停止反演，否则，进入步骤409；δ取值范围为：1.0e^-4～1.0e^-6；

残差数据根据以下公式进行计算：

Δd_z＝d_zobs-d_zcal；

Δd_x＝d_xobs-d_xcal；

Δd_y＝d_yobs-d_ycal；

其中Δd＝(Δd_x，Δd_y，Δd_z)^t为残差数据，d_obs＝(d_xobs，d_yobs，d_zobs)^t为测量的三分量炮集数据，d_cal＝(d_xcal，d_ycal，d_zcal)^t为模拟的三分量炮集数据，上标t表示转置；

所述误差能量根据以下公式进行计算：

$E=\frac{1}{2}{\mathrm{\&Delta;d}}^t{C}_D^{-1}\mathrm{\&Delta;}d;$

其中，

其中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵，X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，在二维情况，取p≥0.5，在三维情况，取p≥1，本实施例中，取p＝0.6。

在本实施例中，给定δ＝1.0e^-6，当E<δ时，输出反演结果，并停止反演，否则，进入步骤409。

步骤409，将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；

解所述弹性波动方程，按反时间方向进行波场传播，具体为：

将计算的残差数据作为输入数据，用HTI介质中弹性波动方程和伪谱法计算剩余地震波场，并按反时间方向传播。

步骤410，利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，计算反演公式中的系数张量C_pjkl(x),下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向；

具体计算方法如下：

$C_{1111}\left(x\right)=a{\cos }^4\mathrm{\&theta;}+b{\sin }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

$C_{2222}\left(x\right)=a{\sin }^4\mathrm{\&theta;}+b{\cos }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

C₃₃₃₃(x)＝a；

C₂₃₂₃(x)＝C₂₃₃₂(x)＝C₃₂₂₃(x)＝C₃₂₃₂(x)＝ccos²θ；

C₁₃₁₃(x)＝C₁₃₃₁(x)＝C₃₁₁₃(x)＝C₃₁₃₁(x)＝csin²θ；

$C_{1212}\left(x\right)=C_{1221}\left(x\right)=C_{2112}\left(x\right)=C_{2121}\left(x\right)=\frac{1}{4}(a+b-2d){\sin }^{2}2\mathrm{\&theta;}+c{\cos }^{2}2\mathrm{\&theta;};$

$C_{1122}\left(x\right)=C_{2211}\left(x\right)=\frac{1}{4}(a+b-4c){\sin }^{2}2\mathrm{\&theta;}+d({\sin }^{4}2\mathrm{\&theta;}+{\cos }^{4}2\mathrm{\&theta;});$

C₁₁₃₃(x)＝C₃₃₁₁(x)＝gcos²θ+dsin²θ；

C₂₂₃₃(x)＝C₂₃₃₂₂(x)＝gsin²θ+dcos²θ；

$C_{1112}\left(x\right)=C_{1121}\left(x\right)=C_{1211}\left(x\right)=C_{2111}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{2212}\left(x\right)=C_{2221}\left(x\right)=C_{1222}\left(x\right)=C_{2122}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\sin }^2\mathrm{\&theta;}-b{\cos }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{3312}\left(x\right)=C_{3321}\left(x\right)=C_{1233}\left(x\right)=C_{2133}\left(x\right)=\frac{1}{2}(g-d)sin2\mathrm{\&theta;};$

$\begin{array}{c}{C}_{2331}\left(x\right)=C_{2313}\left(x\right)=C_{3231}\left(x\right)=C_{3213}\left(x\right)=C_{1332}\left(x\right)=C_{3132}\left(x\right)=C_{1323}\left(x\right)=C_{3123}\left(x\right)\\ =\frac{1}{2}{\sin }_2\left(x\right)\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;\end{array};$

其中a、b、c、d、g为系数，θ为裂隙裂缝走向；除上述各元素外，张量C_pjkl(x)的其它元素值为零；

其中，系数a、b、c、d、g通过以下公式进行计算：

$a=\frac{{\mathrm{\&lambda;}}^2}{(1+E_N)(\mathrm{\&lambda;}+2\mathrm{\&mu;})}(\frac{1}{(1+E_N)}-1)S;$

$b=-\frac{(\mathrm{\&lambda;}+2\mathrm{\&mu;})}{{(1+E_N)}^2}S;$

$c=-\frac{\mathrm{\&mu;}}{{(1+E_T)}^2}T;$

$d=-\frac{\mathrm{\&lambda;}}{{(1+E_N)}^2}S;$

$g=\frac{-{\mathrm{\&lambda;}}^2(1+E_N)}{{(1+E_N)}^2(\mathrm{\&lambda;}+2\mathrm{\&mu;})}S;$

其中，S、T、E_N、E_T通过以下公式进行计算：

γ＝μ/(λ+2μ)；

$S=\frac{16}{3(3-2\mathrm{\&gamma;})};$

$T=\frac{4}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

$E_N=\frac{16e}{3(3-2\mathrm{\&gamma;})};$

$E_T=\frac{4e}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

其中，λ为拉梅常数、μ为剪切模量，e为裂隙裂缝密度，α为裂隙/缝纵横比，k′为填充流体的体积模量。

步骤411，将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向；T为地震记录长度，本实施例中，T＝500ms；N_s为用于反演的炮的总数,本实施例中，N_s＝15；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……，Ns；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

步骤412，迭代修改裂隙密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为所述第n次迭代修改后的裂隙密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙密度梯度；

其中，β_n是利用抛物线法进行计算得到的。

步骤413，将修改后的裂隙密度模型作为初始裂隙密度模型，进入步骤407。

根据图2虚线所示，为本实施例所给出的初始数据，经20次迭反演后，结果如图2和图3所示。图3的点划线为反演的裂隙密度沿裂隙层某深度变化曲线，与实线所示的实际裂隙密度比较，反演结果很好地预测了地层裂隙密度值。图2为断层左侧裂隙密度抽道显示，图中点划线为反演结果，实线为实际值，比较可见，除界面处外，预测值与实际基本一致。

实施例二

为一五层地质模型，模型大小和网格大小同实施例一。第三层和第四层发育有裂隙，第三层的裂隙密度为0.12，裂隙方向与x₁轴夹角为-35°，第四层的裂隙密度为0.06，裂隙方向与x₁轴夹角为-75°，实施例2的模型如图7所示，其相应参数见表2。

层号	厚度(m)	p(g/cm2)	Vp(m/s)	Vs(m/s)	e	θ(度)
							1	400	2400	3640	1750	0	0
2	70	2300	3157	1579	0	0
							3	50	2450	3513	1887	0.12	-35
4	80	2500	3726	2053	0.06	-75
							5	600	2400	3640	1750	0	0

表2

该模型上采集了30炮三分量地震数据，图8为其中某炮的数据，进行反演，具体实现步骤与实施例1基本相同。不同之处有点：

(1)模型上采集了30炮,每炮101个检波点，检波点间距10m，最小炮检距0m，最大炮检距1000m，时间采样率1ms，采集长度600ms。

(2)实施例一步骤106中改用纵波时差速度椭圆法获取裂隙裂缝走向。

(3)实施例一步骤107中初始裂隙密度模型取值如图9中①线所示。

图9中还包括20次迭反演后反演的裂隙密度图，③表示真实数据，②为反演结果。比较两者可以看出反演结果除了一些高频抖动外，两层反演的裂缝密度与真实数据非常相近，表明本发明的裂隙裂缝密度反演方法很好地获得了地层的真实密度。

实施例三

如图10所示，为本发明一种地震地层裂隙裂缝密度反演系统的结构图，包括以下单元：

预处理单元1001，用于采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

对所述三分量地震数据进行预处理包括：两水平分量旋转，三分量静校正、地表一致性振幅补偿和叠前去除噪音。

深度域纵波速度及地层质量密度计算单元1002，用于在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并计算地层的质量密度；

包括：

时间域纵波速度计算子单元，用于根据以下公式计算所述时间域纵波速度：

$V_{pi}^2=\frac{V_{r,i}^2{t}_{pi}-V_{r,i-1}^2{t}_{p(i-1)}}{t_{pi}-t_{p(i-1)}};$

其中t_pi为第i层底纵波双程旅行时，V_pi为第i层的纵波速度，V_r,i为第i层纵波均方根速度；

需要说明的是，每一个层都有顶和底，因而第i层也会有顶和底。

纵波时间-深度转换子单元，用于根据以下公式将时间域的纵波速度转换成深度域纵波速度：

$h_i=\frac{1}{2}(t_{pi}-t_{p(i-1)})V_{pi};$

其中t_pi为第i层底纵波双程旅行时，h_i为第i层的厚度，V_pi为第i层纵波速度；

地层的质量密度计算子单元，用于根据以下公式计算所述地层的质量密度：

$\mathrm{\&rho;}=0.31V_p^{0.25};$

其中ρ为质量密度，V_p为纵波速度；

所述深度域纵波速度为所有地层纵波速度的集合：

V_p＝(V_p1,V_p2,…,V_pi,…,V_pN)，其中，N为模型中地层的层数，V_pi为第i层纵波速度。

深度域横波速度计算单元1003，用于在所述三分量炮集数据中径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度；

包括：

时间域横波速度计算子单元，用于根据以下公式计算所述时间域横波速度：

$V_{si}=\frac{V_{pi}}{r_i};$

其中r_i为所述转换波速度分析获得的第i层纵横波速度比，V_pi为第i层纵波速度，V_si为第i层横波层速度；

横波时间-深度转换子单元，用于根据以下公式将时间域的横波层速度转换成深度域横波速度：

h_i＝(t_si-t_s(i-1))V_si；

其中，t_si为第i层底横波旅行时，h_i为第i层的厚度，V_si为第i层横波层速度；

需要说明的是，每一个层都有顶和底，因而第i层也会有顶和底。

第i层底横波旅行时计算子单元，用于根据以下公式计算所述第i层底横波旅行时：

$t_{si}=t_{ci}-\frac{1}{2}{t}_{pi};$

其中，t_ci为第i层转换波旅行时，t_pi为第i层纵波双程旅行时。

初始模型计算单元1004，用根据所述深度域纵波速度、所述深度域横波速度以及所述 地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

包括：

拉梅常数计算子单元，用于根据以下公式计算所述拉梅常数：

$\mathrm{\&lambda;}=\mathrm{\&rho;}(V_p^2-2V_s^2);$

其中λ为拉梅常数，ρ为质量密度，V_p为纵波速度，V_s为横波速度；

剪切模量计算子单元，用于根据以下公式计算所述剪切模量：

μ＝ρV_s ²；

其中μ为剪切模量，ρ为质量密度，Vs为横波速度，为所有地层横波速度的集合：

V_s＝(V_s1,V_s2,…,V_si,…,V_sN)；其中，N为模型中地层的层数，V_si为第i层横波速度。

模型计算单元1005，用于利用所述拉梅常数、剪切模量和地层质量密度的初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

裂隙走向模型计算单元1006，用于通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙裂缝走向模型，其中，所述裂隙裂缝走向与正北方向的夹角记为θ；

正演模拟单元1007，用于利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

所述地震波场正演模拟具体为：利用HTI介质中弹性波动方程，通过伪谱法进行模拟；

其中，震源函数采用零相位雷克子波。

反演结果输出单元1008，用于计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，当E<δ时，输出反演结果，并停止反演，否则，进入逆时传播单元；δ取值范围为：1.0e^-4～1.0e^-6；

包括：

残差数据计算子单元，用于据根据以下公式计算残差数据：

Δd_z＝d_zobs-d_zcal；

Δd_x＝d_xobs-d_xcal；

Δd_y＝d_yobs-d_ycal；

其中Δd＝(Δd_x，Δd_y，Δd_z)^t为残差数据，d_obs＝(d_xobs，d_yobs，d_zobs)^t为测量的三分量炮 集数据，d_cal＝(d_xcal，d_ycal，d_zcal)^t为模拟的三分量炮集数据，上标t表示转置；

误差能量计算子单元，用于根据以下公式计算误差能量：

$E=\frac{1}{2}{\mathrm{\&Delta;d}}^t{C}_D^{-1}\mathrm{\&Delta;}d;$

其中，

其中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵，X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，在二维情况，取p≥0.5，在三维情况，取p≥1。

逆时传播单元1009，用于将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；即：解所述弹性波动方程，按反时间方向进行波场传播。

系数张量计算单元1010，用于利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，计算反演公式中的系数张量C_pjkl(x)(其中下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向)；

系数张量计算子单元，用于根据以下公式计算反演公式中的系数张量C_pjkl(x)：

$C_{1111}\left(x\right)=a{\cos }^4\mathrm{\&theta;}+b{\sin }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

$C_{2222}\left(x\right)=a{\sin }^4\mathrm{\&theta;}+b{\cos }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

C₃₃₃₃(x)＝a；

C₂₃₂₃(x)＝C₂₃₃₂(x)＝C₃₂₂₃(x)＝C₃₂₃₂(x)＝ccos²θ；

C₁₃₁₃(x)＝C₁₃₃₁(x)＝C₃₁₁₃(x)＝C₃₁₃₁(x)＝csin²θ；

$C_{1212}\left(x\right)=C_{1221}\left(x\right)=C_{2112}\left(x\right)=C_{2121}\left(x\right)=\frac{1}{4}(a+b-2d){\sin }^{2}2\mathrm{\&theta;}+c{\cos }^{2}2\mathrm{\&theta;};$

$C_{1122}\left(x\right)=C_{2211}\left(x\right)=\frac{1}{4}(a+b-4c){\sin }^{2}2\mathrm{\&theta;}+d({\sin }^{4}2\mathrm{\&theta;}+{\cos }^{4}2\mathrm{\&theta;});$

C₁₁₃₃(x)＝C₃₃₁₁(x)＝gcos²θ+dsin²θ；

C₂₂₃₃(x)＝C₂₃₃₂₂(x)＝gsin²θ+dcos²θ；

$C_{1112}\left(x\right)=C_{1121}\left(x\right)=C_{1211}\left(x\right)=C_{2111}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{2212}\left(x\right)=C_{2221}\left(x\right)=C_{1222}\left(x\right)=C_{2122}\left(x\right)=\frac{1}{2}\sin 2\mathrm{\&theta;}\&lsqb;a{\sin }^2\mathrm{\&theta;}-b{\cos }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{3312}\left(x\right)=C_{3321}\left(x\right)=C_{1233}\left(x\right)=C_{2133}\left(x\right)=\frac{1}{2}(g-d)sin2\mathrm{\&theta;};$

$\begin{array}{c}{C}_{2331}\left(x\right)=C_{2313}\left(x\right)=C_{3231}\left(x\right)=C_{3213}\left(x\right)=C_{1332}\left(x\right)=C_{3132}\left(x\right)=C_{1323}\left(x\right)=C_{3123}\left(x\right)\\ =\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;\end{array};$

其中a、b、c、d、g为系数，θ为裂隙裂缝走向；除上述各元素外，张量C_pjkl(x)的其它元素值为零；

系数张量系数计算子单元，用于通过以下公式计算系数a、b、c、d、g：

$a=\frac{{\mathrm{\&lambda;}}^2}{(1+E_N)(\mathrm{\&lambda;}+2\mathrm{\&mu;})}(\frac{1}{(1+E_N)}-1)S;$

$b=-\frac{(\mathrm{\&lambda;}+2\mathrm{\&mu;})}{{(1+E_N)}^2}S;$

$c=-\frac{\mathrm{\&mu;}}{{(1+E_T)}^2}T;$

$d=-\frac{\mathrm{\&lambda;}}{{(1+E_N)}^2}S;$

$g=\frac{-{\mathrm{\&lambda;}}^2(1+E_N)}{{(1+E_N)}^2(\mathrm{\&lambda;}+2\mathrm{\&mu;})}S;$

其中，S、T、E_N、E_T通过以下公式进行计算：

γ＝μ/(λ+2μ)；

$S=\frac{16}{3(3-2\mathrm{\&gamma;})};$

$T=\frac{4}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

$E_N=\frac{16e}{3(3-2\mathrm{\&gamma;})};$

$E_T=\frac{4e}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

其中，λ为拉梅常数、μ为剪切模量，e为裂隙裂缝密度，α为裂隙/缝纵横比，k′为填充流体的体积模量。

裂隙裂缝密度梯度值计算单元1011，用于将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向；T为地震记录长度，Ns为用于反演的炮的总数；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……，Ns；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

迭代修改单元1012，用于根据所述裂隙裂缝密度梯度，按以下公式迭代修改裂隙裂缝密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为第n次迭代修改后的裂隙裂缝密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙裂缝密度梯度；

循环单元1013，用于将修改后的裂隙裂缝密度模型作为初始裂隙裂缝密度模型，即取裂隙裂缝密度e₀＝e_n+1,进入正演模拟单元。

通过以上实施例的说明，本发明能够达到以下有益效果：全波形反演充分利用了采集到的地震波信息，并自动考虑了地震波在地下的传播行为，能给出高精度地震反演结果。本发明通过HTI介质中弹性全波形反演理论，提出了全波形裂隙裂缝密度反演方法，实现了地震地下地层裂隙裂缝密度的定量预测。

本领域一般技术人员在此设计思想之下所做任何不具有创造性的改造，均应视为在本发明的保护范围之内。

Claims (16)

1.一种地震地层裂隙裂缝密度反演方法，其特征在于，包括以下步骤：

步骤一，采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

步骤二，在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并计算地层的质量密度；

步骤三，在所述三分量炮集数据中径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度，其中，所述时间域横波速度根据以下公式进行：

$V_{si}=\frac{V_{pi}}{r_i};$

其中r_i为所述转换波速度分析获得的第i层纵横波速度比，V_pi为第i层纵波速度，V_si为第i层横波层速度；

所述横波时间-深度转换具体为：根据以下公式将时间域的横波层速度转换成深度域横波速度：

h_i＝(t_si-t_s(i-1))V_si；

其中，t_si为第i层底横波旅行时，h_i为第i层的厚度，V_si为第i层横波层速度；

所述第i层底横波旅行时根据以下公式计算：

$t_{si}=t_{ci}-\frac{1}{2}{t}_{pi};$

其中，t_ci为第i层转换波旅行时，t_pi为第i层纵波双程旅行时；

步骤四，根据所述深度域纵波速度、所述深度域横波速度以及所述地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

步骤五，利用所述拉梅常数、剪切模量和地层质量密度为初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

步骤六，通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙裂缝走向模型，其中，所述裂隙裂缝走向与正北方向的夹角记为θ；

步骤七，利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

步骤八，计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，当E<δ时，输出反演结果，并停止反演，否则，进入步骤九；δ取值范围为：1.0e^-4～1.0e^-6；

步骤九，将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；

步骤十，利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，计算反演公式中的系数张量C_pjkl(x),下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向，其中，计算反演公式中的系数张量C_pjkl(x)的方法具体如下：

$C_{1111}\left(x\right)=a{\cos }^4\mathrm{\&theta;}+b{\sin }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

$C_{2222}\left(x\right)=a{\sin }^4\mathrm{\&theta;}+b{\cos }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

C₃₃₃₃(x)＝a；

C₂₃₂₃(x)＝C₂₃₃₂(x)＝C₃₂₂₃(x)＝C₃₂₃₂(x)＝c cos²θ；

C₁₃₁₃(x)＝C₁₃₃₁(x)＝C₃₁₁₃(x)＝C₃₁₃₁(x)＝c sin²θ；

$C_{1212}\left(x\right)=C_{1221}\left(x\right)=C_{2112}\left(x\right)=C_{2121}\left(x\right)=\frac{1}{4}(a+b-2d){\sin }^{2}2\mathrm{\&theta;}+c{\cos }^{2}2\mathrm{\&theta;};$

$C_{1122}\left(x\right)=C_{2211}\left(x\right)=\frac{1}{4}(a+b-4c){\sin }^{2}2\mathrm{\&theta;}+d({\sin }^{4}2\mathrm{\&theta;}+{\cos }^{4}2\mathrm{\&theta;});$

C₁₁₃₃(x)＝C₃₃₁₁(x)＝g cos²θ+d sin²θ；

C₂₂₃₃(x)＝C₂₃₃₂₂(x)＝g sin²θ+d cos²θ；

$C_{1112}\left(x\right)=C_{1121}\left(x\right)=C_{1211}\left(x\right)=C_{2111}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{2212}\left(x\right)=C_{2221}\left(x\right)=C_{1222}\left(x\right)=C_{2122}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\sin }^2\mathrm{\&theta;}-b{\cos }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{3312}\left(x\right)=C_{3321}\left(x\right)=C_{1233}\left(x\right)=C_{2133}\left(x\right)=\frac{1}{2}(g-d)sin2\mathrm{\&theta;};$

$\begin{array}{c}{C}_{2331}\left(x\right)=C_{2313}\left(x\right)=C_{3231}\left(x\right)=C_{3213}\left(x\right)=C_{1332}\left(x\right)=C_{3132}\left(x\right)=C_{1323}\left(x\right)=C_{3123}\left(x\right)\\ =\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;\end{array};$

其中a、b、c、d、g为系数，θ为裂隙裂缝走向；除上述各元素外，系数张量C_pjkl(x)的其它元素值为零；

步骤十一，将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以步骤10获得的系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向；T为地震记录长度，Ns为用于反演的炮的总数；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……，N_s；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

步骤十二，根据所述裂隙裂缝密度梯度，按以下公式迭代修改裂隙裂缝密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为第n次迭代修改后的裂隙裂缝密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙裂缝密度梯度；

步骤十三，将修改后的裂隙裂缝密度模型作为初始裂隙裂缝密度模型，即取裂隙裂缝密度e₀＝e_n+1,进入步骤七。

2.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于：

所述对所述三分量地震数据进行预处理包括：两水平分量旋转，三分量静校正、地表一致性振幅补偿和叠前去除噪音。

3.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于：

所述时间域纵波速度计算根据以下公式进行：

$V_{pi}^2=\frac{V_{r,i}^2{t}_{pi}-V_{r,i-1}^2{t}_{p(i-1)}}{t_{pi}-t_{p(i-1)}};$

其中t_pi为第i层底纵波双程旅行时，V_pi为第i层的纵波速度，V_r,i为第i层纵波均方根速度；

所述纵波时间-深度转换具体为：根据以下公式将时间域的纵波速度转换成深度域纵波速度：

$h_i=\frac{1}{2}(t_{pi}-t_{p(i-1)})V_{pi};$

其中t_pi为第i层底纵波双程旅行时，h_i为第i层的厚度，V_pi为第i层纵波速度；

所述计算地层的质量密度，根据以下公式计算：

$\mathrm{\&rho;}=0.31V_p^{0.25};$

其中ρ为质量密度，V_p为纵波速度；

所述深度域纵波速度为所有地层纵波速度的集合：

V_p＝(V_p1,V_p2,…,V_pi,…,V_pN)，其中，N为模型中地层的层数，V_pi为第i层纵波速度。

4.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于：

所述拉梅常数根据以下公式计算：

$\mathrm{\&lambda;}=\mathrm{\&rho;}(V_p^2-2V_s^2);$

其中λ为拉梅常数，ρ为质量密度，V_p为纵波速度，V_s为横波速度；

所述剪切模量根据以下公式计算：

$\mathrm{\&mu;}={\mathrm{\&rho;V}}_s^2;$

其中μ为剪切模量，ρ为质量密度，Vs为横波速度，为所有地层横波速度的集合，公式表述如下：

V_s＝(V_s1,V_s2,…,V_si,…,V_sN)；其中，N为模型中地层的层数，V_si为第i层横波速度。

5.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于，

所述地震波场正演模拟具体为：利用水平横向各向同性HTI介质中弹性波动方程，通过伪谱法进行模拟；

其中，震源函数采用零相位雷克子波。

6.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于：

所述残差数据根据以下公式进行计算：

Δd_z＝d_zobs-d_zcal；

Δd_x＝d_xobs-d_xcal；

Δd_y＝d_yobs-d_ycal；

其中Δd＝(Δd_x，Δd_y，Δd_z)^t为残差数据，d_obs＝(d_xobs，d_yobs，d_zobs)^t为测量的三分量炮集数据，d_cal＝(d_xcal，d_ycal，d_zcal)^t为模拟的三分量炮集数据，上标t表示转置；

所述误差能量根据以下公式进行计算：

$E=\frac{1}{2}{\mathrm{\&Delta;d}}^t{C}_D^{-1}\mathrm{\&Delta;}d;$

其中，

其中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵，X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，在二维情况，取p≥0.5，在三维情况，取p≥1。

7.如权利要求5所述的地震地层裂隙裂缝密度反演方法，其特征在于，

所述逆时传播具体为：解所述弹性波动方程，按反时间方向进行波场传播。

8.如权利要求1所述的地震地层裂隙裂缝密度反演方法，其特征在于，所述系数a、b、c、d、g通过以下公式进行计算：

$a=\frac{{\mathrm{\&lambda;}}^2}{(1+E_N)(\mathrm{\&lambda;}+2\mathrm{\&mu;})}(\frac{1}{(1+E_N)}-1)S;$

$b=-\frac{(\mathrm{\&lambda;}+2\mathrm{\&mu;})}{{(1+E_N)}^2}S;$

$c=-\frac{\mathrm{\&mu;}}{{(1+E_T)}^2}T;$

$d=-\frac{\mathrm{\&lambda;}}{{(1+E_N)}^2}S;$

$g=\frac{-{\mathrm{\&lambda;}}^2(1+E_N)}{{(1+E_N)}^2(\mathrm{\&lambda;}+2\mathrm{\&mu;})}S;$

其中，S、T、E_N、E_T通过以下公式进行计算：

γ＝μ/(λ+2μ)；

$S=\frac{16}{3(3-2\mathrm{\&gamma;})};$

$T=\frac{4}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

$E_N=\frac{16e}{3(3-2\mathrm{\&gamma;})};$

$E_T=\frac{4e}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

其中，λ为拉梅常数、μ为剪切模量，e为裂隙裂缝密度，α为裂隙/缝纵横比，k′为填充流体的体积模量。

9.一种地震地层裂隙裂缝密度反演系统，其特征在于，包括以下单元：

预处理单元，用于采集三分量地震数据，并对所述三分量地震数据进行预处理，获得测量三分量炮集数据；

深度域纵波速度及地层质量密度计算单元，用于在所述三分量炮集数据中的垂直分量炮集数据中抽取共中心点道集数据，进行纵波速度分析、时间域纵波速度计算和纵波时间-深度转换，获得深度域纵波速度，并计算地层的质量密度；

深度域横波速度计算单元，用于在所述三分量炮集数据中径向水平分量炮集数据中抽取共转换点道集数据，进行转换波速度分析和横波速度计算，将转换波时间域横波速度通过横波时间-深度转换，转换到深度域，得到深度域横波速度，其中，所述深度域横波速度计算单元包括：

时间域横波速度计算子单元，用于根据以下公式计算所述时间域横波速度：

$V_{si}=\frac{V_{pi}}{r_i};$

其中r_i为所述转换波速度分析获得的第i层纵横波速度比，V_pi为第i层纵波速度，V_si为第i层横波层速度；

横波时间-深度转换子单元，用于根据以下公式将时间域的横波层速度转换成深度域横波速度：

h_i＝(t_si-t_s(i-1))V_si；

其中，t_si为第i层底横波旅行时，h_i为第i层的厚度，V_si为第i层横波层速度；

第i层底横波旅行时计算子单元，用于根据以下公式计算所述第i层底横波旅行时：

$t_{si}=t_{ci}-\frac{1}{2}{t}_{pi};$

其中，t_ci为第i层转换波旅行时，t_pi为第i层纵波双程旅行时；

初始模型计算单元，用根据所述深度域纵波速度、所述深度域横波速度以及所述地层质量密度计算拉梅常数和剪切模量，获得所述地层质量密度、拉梅常数和剪切模量的初始模型；

模型计算单元，用于利用所述拉梅常数、剪切模量和地层质量密度为初始模型，通过多分量数据处理方法获取地层拉梅常数、剪切模量和质量密度模型；

裂隙走向模型计算单元，用于通过纵波方位速度椭圆法或纵波方位振幅法获得裂隙裂缝走向模型，其中，所述裂隙裂缝走向与正北方向的夹角记为θ；

正演模拟单元，用于利用所述拉梅常数、剪切模量和地层质量密度模型、所述裂隙裂缝走向模型以及初始裂隙裂缝密度模型，进行地震波场正演模拟，获得正演模拟三分量地震波场数据和模拟三分量炮集数据，其中，所述初始裂隙裂缝密度模型的裂隙裂缝密度取e₀＝0.01；

反演结果输出单元，用于计算所述模拟三分量炮集数据与所述测量三分量炮集数据的差，即残差数据，并计算误差能量E，当E<δ时，输出反演结果，并停止反演，否则，进入逆时传播单元；δ取值范围为：1.0e^-4～1.0e^-6；

逆时传播单元，用于将所述残差数据作为源数据，通过与所述正演模拟过程相反的逆时传播，获得所述残差数据的剩余地震波场；

系数张量计算单元，用于利用所述地层质量密度、拉梅常数、剪切模量、裂隙裂缝走向数据和所述初始裂隙裂缝密度模型，根据以下公式计算反演公式中的系数张量C_pjkl(x)(其中下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，这里1、2和3分别代表x、y和z三个方向)：

$C_{1111}\left(x\right)=a{\cos }^4\mathrm{\&theta;}+b{\sin }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

$C_{2222}\left(x\right)=a{\sin }^4\mathrm{\&theta;}+b{\cos }^4\mathrm{\&theta;}+(\frac{1}{2}d+c){\sin }^{2}2\mathrm{\&theta;};$

C₃₃₃₃(x)＝a；

C₂₃₂₃(x)＝C₂₃₃₂(x)＝C₃₂₂₃(x)＝C₃₂₃₂(x)＝c cos²θ；

C₁₃₁₃(x)＝C₁₃₃₁(x)＝C₃₁₁₃(x)＝C₃₁₃₁(x)＝c sin²θ；

$C_{1212}\left(x\right)=C_{1221}\left(x\right)=C_{2112}\left(x\right)=C_{2121}\left(x\right)=\frac{1}{4}(a+b-2d){\sin }^{2}2\mathrm{\&theta;}+c{\cos }^{2}2\mathrm{\&theta;};$

$C_{1122}\left(x\right)=C_{2211}\left(x\right)=\frac{1}{4}(a+b-4c){\sin }^{2}2\mathrm{\&theta;}+d({\sin }^{4}2\mathrm{\&theta;}+{\cos }^{4}2\mathrm{\&theta;});$

C₁₁₃₃(x)＝C₃₃₁₁(x)＝g cos²θ+d sin²θ；

C₂₂₃₃(x)＝C₂₃₃₂₂(x)＝g sin²θ+d cos²θ；

$C_{1112}\left(x\right)=C_{1121}\left(x\right)=C_{1211}\left(x\right)=C_{2111}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{2212}\left(x\right)=C_{2221}\left(x\right)=C_{1222}\left(x\right)=C_{2122}\left(x\right)=\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\sin }^2\mathrm{\&theta;}-b{\cos }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;;$

$C_{3312}\left(x\right)=C_{3321}\left(x\right)=C_{1233}\left(x\right)=C_{2133}\left(x\right)=\frac{1}{2}(g-d)\sin 2\mathrm{\&theta;};$

$\begin{array}{c}{C}_{2331}\left(x\right)=C_{2313}\left(x\right)=C_{3231}\left(x\right)=C_{3213}\left(x\right)=C_{1332}\left(x\right)=C_{3132}\left(x\right)=C_{1323}\left(x\right)=C_{3123}\left(x\right)\\ =\frac{1}{2}sin2\mathrm{\&theta;}\&lsqb;a{\cos }^2\mathrm{\&theta;}-b{\sin }^2\mathrm{\&theta;}-(d+2c)cos2\mathrm{\&theta;}\&rsqb;\end{array};$

其中a、b、c、d、g为系数，θ为裂隙裂缝走向；除上述各元素外，系数张量C_pjkl(x)的其它元素值为零；

裂隙裂缝密度梯度值计算单元，用于将所述正演模拟三分量地震波场的三分量数据和所述剩余地震波场的三分量数据分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，再乘以系数张量，得到裂隙裂缝密度梯度值：

$\mathrm{\&delta;}e\left(x\right)=-C_{pjkl}\left(x\right)\underset{s=1}{\overset{N_s}{\&Sigma;}}{\&Integral;}_0^{T}dt\frac{\&part;{\mathrm{\&psi;}}_j}{\&part;x_p}(x,t;x_s)\frac{\&part;u_k}{\&part;x_l}(x,t;x_s);$

其中u_k为正演波场在k方向的分量,ψ_j为回传的剩余波场在j方向的分量；下标p＝(1,2,3)，j＝(1,2,3)，k＝(1,2,3)，l＝(1,2,3)，其中1、2和3分别表示x、y和z三个方向；T为地震记录长度，Ns为用于反演的炮的总数；向量x表示地下某点空间位置，x_s为震源点空间位置，t为时间；C_pjkl(x)为系数张量；δe为裂隙裂缝密度梯度；s为炮号，s＝1，2，……，N_s；x_p表示p方向的坐标变量；x_l表示l方向的坐标变量；

迭代修改单元，用于根据所述裂隙裂缝密度梯度，按以下公式迭代修改裂隙裂缝密度值：

e_n+1＝e_n-β_nδe_n；

其中，e_n+1为第n+1次迭代修改后的裂隙裂缝密度，e_n为第n次迭代修改后的裂隙裂缝密度值，β_n为第n次迭代修改步长，δe_n为第n次迭代后的裂隙裂缝密度梯度；

循环单元，用于将修改后的裂隙裂缝密度模型作为初始裂隙裂缝密度模型，即取裂隙裂缝密度e₀＝e_n+1,进入正演模拟单元。

10.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于：

所述预处理单元中，对所述三分量地震数据进行预处理包括：两水平分量旋转，三分量静校正、地表一致性振幅补偿和叠前去除噪音。

11.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于，所述深度域纵波速度及地层质量密度计算单元包括：

时间域纵波速度计算子单元，用于根据以下公式计算所述时间域纵波速度：

$V_{pi}^2=\frac{V_{r,i}^2{t}_{pi}-V_{r,i-1}^2{t}_{p(i-1)}}{t_{pi}-t_{p(i-1)}};$

其中t_pi为第i层底纵波双程旅行时，V_pi为第i层的纵波速度，V_r,i为第i层纵波均方根速度；

纵波时间-深度转换子单元，用于根据以下公式将时间域的纵波速度转换成深度域纵波速度：

$h_i=\frac{1}{2}(t_{pi}-t_{p(i-1)})V_{pi};$

其中t_pi为第i层底纵波双程旅行时，h_i为第i层的厚度，V_pi为第i层纵波速度；

地层的质量密度计算子单元，用于根据以下公式计算所述地层的质量密度：

$\mathrm{\&rho;}=0.31V_p^{0.25};$

其中ρ为质量密度，V_p为纵波速度；

所述深度域纵波速度为所有地层纵波速度的集合：

V_p＝(V_p1,V_p2,…,V_pi,…,V_pN)，其中，N为模型中地层的层数，V_pi为第i层纵波速度。

12.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于，所述模型计算单元包括：

拉梅常数计算子单元，用于根据以下公式计算所述拉梅常数：

$\mathrm{\&lambda;}=\mathrm{\&rho;}(V_p^2-2V_s^2);$

其中λ为拉梅常数，ρ为质量密度，V_p为纵波速度，V_s为横波速度；

剪切模量计算子单元，用于根据以下公式计算所述剪切模量：

$\mathrm{\&mu;}={\mathrm{\&rho;V}}_s^2;$

其中μ为剪切模量，ρ为质量密度，Vs为横波速度，为所有地层横波速度的集合：

V_s＝(V_s1,V_s2,…,V_si,…,V_sN)；其中，N为模型中地层的层数，V_si为第i层横波速度。

13.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于，

所述地震波场正演模拟具体为：利用HTI介质中弹性波动方程，通过伪谱法进行模拟；

其中，震源函数采用零相位雷克子波。

14.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于，所述反演结果输出单元包括：

残差数据计算子单元，用于根据以下公式计算残差数据：

Δd_z＝d_zobs-d_zcal；

Δd_x＝d_xobs-d_xcal；

Δd_y＝d_yobs-d_ycal；

其中Δd＝(Δd_x，Δd_y，Δd_z)^t为残差数据，d_obs＝(d_xobs，d_yobs，d_zobs)^t为测量的三分量炮集数据，d_cal＝(d_xcal，d_ycal，d_zcal)^t为模拟的三分量炮集数据，上标t表示转置；

误差能量计算子单元，用于根据以下公式计算误差能量：

$E=\frac{1}{2}{\mathrm{\&Delta;d}}^t{C}_D^{-1}\mathrm{\&Delta;}d;$

其中，

其中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵，X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，在二维情况，取p≥0.5，在三维情况，取p≥1。

15.如权利要求13所述的地震地层裂隙裂缝密度反演系统，其特征在于，

所述逆时传播单元，用于以所述残差数据为源，利用HTI介质弹性波动方程，按反时间方向进行波场传播。

16.如权利要求9所述的地震地层裂隙裂缝密度反演系统，其特征在于，所述系数张量计算单元包括：

系数张量系数计算子单元，用于通过以下公式计算系数a、b、c、d、g：

$a=\frac{{\mathrm{\&lambda;}}^2}{(1+E_N)(\mathrm{\&lambda;}+2\mathrm{\&mu;})}(\frac{1}{(1+E_N)}-1)S;$

$b=-\frac{(\mathrm{\&lambda;}+2\mathrm{\&mu;})}{{(1+E_N)}^2}S;$

$c=-\frac{\mathrm{\&mu;}}{{(1+E_T)}^2}T;$

$d=-\frac{\mathrm{\&lambda;}}{{(1+E_N)}^2}S;$

$g=\frac{-{\mathrm{\&lambda;}}^2(1+E_N)}{{(1+E_N)}^2(\mathrm{\&lambda;}+2\mathrm{\&mu;})}S;$

其中，S、T、E_N、E_T通过以下公式进行计算：

γ＝μ/(λ+2μ)；

$S=\frac{16}{3(3-2\mathrm{\&gamma;})};$

$T=\frac{4}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

$E_N=\frac{16e}{3(3-2\mathrm{\&gamma;})};$

$E_T=\frac{4e}{3\mathrm{\&gamma;}\&lsqb;1-\mathrm{\&gamma;}+k^{\&prime;}/\left(\mathrm{\&pi;}\mathrm{\&alpha;}\mathrm{\&mu;}\right)\&rsqb;};$

其中，λ为拉梅常数、μ为剪切模量，e为裂隙裂缝密度，α为裂隙/缝纵横比，k′为填充流体的体积模量。