CN103207409A - 一种频率域全波形反演地震速度建模方法 - Google Patents

Abstract

本发明涉及一种频率域全波形反演地震速度建模方法，包括以下步骤：1）获取原始地震炮集记录、震源子波信息和反演所用的初始模型；2）分析步骤1）获得的信息，基于正演算法和优化算法，确定基本的反演参数和从低频到高频的全波形反演框架；3）对不同的频率，计算获得最合适的频域正反演模型网格；4）在低频反演时采用主成分分析方法，压缩参与反演的数据维度；5）判断对应不同频率的投影矩阵维度是否满足阈值转换标准，当满足该转换标准时，进入下一步，否则回到步骤4）；6）引入震源编码方法，并利用随机相位编码方法压制串扰噪声；7）判定是否满足迭代截止条件，如果满足迭代截止条件，则进入下一步；否则回到步骤6）；8）如果没有完成所有频率的反演，则回到步骤3），直至完成所有频率的反演，获得最后的速度模型，输出速度模型。

Description

一种频率域全波形反演地震速度建模方法

技术领域

本发明涉及一种建模方法，特别是关于一种频率域全波形反演地震速度建模方法。

背景技术

由于全波形反演的实现是一个多次迭代优化的过程，从而需要进行大量基于波动方程的波场正演模拟运算，因此计算成本过大是制约该方法应用和发展的重要因素，如何压缩计算成本和提高算法效率一直是国内外关于全波形反演研究的前沿性课题。提高全波形反演计算效率主要有两个途径，一是高效率的优化算法，二是对原始数据的压缩。其中，使用收敛速度较快的优化算法虽然能够减少反演中的迭代次数，例如高斯牛顿法或牛顿法，但是这些方法由于在计算海森矩阵或近似海森矩阵时要投入较大的计算成本，因此对整体的计算效率提升并不明显。而对原始数据空间有效的压缩不仅可以减少方法对存储量的要求，还可以降低迭代中波场正演模拟的计算成本，因此该手段是提高反演效率较为直接有效的途径。特别是将全波形反演问题引入到频率域中，利用几个离散频率从低到高依次反演，不仅在很大程度上压缩了原始数据（时间域）的大小，而且频率域波场正演模拟在多震源情况下能够体现更为高效的计算性能。因此，频率域全波形反演被大多数人认可和研究，然而在反演过程中对应不同的频率尺度所反映出的特性不同，所以该方法还有很多潜力尚未开发。

关于如何降低原始数据空间的大小，目前常用的方法有炮集抽稀方法、频率域正反演以及震源编码技术。对原始数据的压缩必然会导致反演所用信息量的减少，从而降低了反演结果的精度。因此在计算效率和计算精度中是否具有较好的折衷是评价各个方法在应用中的关键因素。对原始资料进行炮集抽稀虽然可以较为直接地压缩原始数据，但是原始信息的损失比例和效率的提升比是相同的，从而导致了反演结果分辨率的降低。因此，在使用该方法之前需要进行分辨率估计或者地下照明度预测等，除非原始资料的冗余度较大（炮间距或道间距非常小），不然该方法并不适合用于全波形反演。

将震源编码技术应用到全波形反演中，通过将多个独立炮集编码组合为超级炮集，降低了原始数据的大小，同时也减少了参与正演的震源个数。但该方法受串扰噪声影响较大，反演结果的分辨率不足，另外该方法受环境噪声的影响也较大。

发明内容

针对上述问题，本发明的目的是提供一种能够联合各种方法，高效率的频率域全波形反演地震速度建模方法。

为实现上述目的，本发明采取以下技术方案：1、一种频率域全波形反演地震速度建模方法，包括以下步骤：1）获取原始地震炮集记录、震源子波信息和反演所用的初始模型；2）分析步骤1）获得的信息，基于正演算法和优化算法，确定基本的反演参数和从低频到高频的全波形反演框架；3）对不同的频率，计算获得最合适的频域正反演模型网格；4）在低频反演时采用主成分分析方法，压缩参与反演的数据维度；5）判断对应不同频率的投影矩阵维度是否满足阈值转换标准，当满足该转换标准时，进入下一步，否则回到步骤4）；6）引入震源编码方法，并利用随机相位编码方法压制串扰噪声；7）判定是否满足迭代截止条件，如果满足迭代截止条件，则进入下一步；否则回到步骤6）；8）如果没有完成所有频率的反演，则回到步骤3），直至完成所有频率的反演，获得最后的速度模型，输出速度模型。

其中步骤2）中，所述正演算法采用基于混合网格的有限差分方法，所述优化算法采用共轭梯度法。

其中步骤3）中，对不同的频率采用下面公式求其所对应的最适网格间距Δh：

$\left\{\begin{array}{c}\mathrm{\Δh}\≤\frac{v_{\min }}{M_{\min }{f}_{\mathrm{ifreq}},}\\ \mathrm{\Δh}=2^N\·\mathrm{\ΔS}\end{array}\right.---\left(1\right)$

式中，v_min为介质最小速度，M_min为一个波长所包含网格的最少个数，由正演算法的精度控制，f_ifreq为反演对应频率，ifreq为反演频率的索引，Δs为炮间距，N为整数。

所述步骤4）中，在使用主成分分析时，首先对频率域波场残差矩阵进行奇异值分解，然后通过积累能量分析和主能量成分提取确定一个投影矩阵，进而将原始数据投影至低维空间上，实现对原始数据的压缩；主成分分析方法的表达式为：

$\left\{\begin{array}{c}\mathrm{\δD}={\mathrm{LXR}}^H\\ \frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_i^2}{\underset{i}{\mathrm{\Σ}}{\mathrm{\σ}}_i^2}\≥\mathrm{\χ}\\ D_{\mathrm{obs}}^k=D_{\mathrm{obs}}{D}^k,S^k={\mathrm{SR}}^k\end{array}\right.,---\left(2\right)$

式中，δD为波场残差矩阵，L和R分别表示左、右奇异值向量矩阵，R^H为R的共轭转置，X为对角奇异值（σ_i）矩阵，i为奇异值矩阵的维度索引，χ为控制投影矩阵维度k的阈值参数，R^k为投影矩阵，D_obs和S分别表示观测数据以及原始震源矩阵。

其中步骤5）中，所述阈值转换标准的表达式为：

k(f_ifreq)＜μM,         （3）

式中，k(f_ifreq)表示对应频率f_ifreq得到的投影矩阵维度，M表示原始震源矩阵维度，μ为预设的最小压缩比的倒数，即控制转换标准的阈值；

其中步骤6中，所述随机相位震源编码方法是通过将多个独立震源编码组合成为一个“超级炮”，同时在迭代过程中随机生成编码序列来压制震源组合所产生的串扰噪声，震源编码的表达式为：

$\left\{\begin{array}{c}\stackrel{\‾}{s}=\underset{j=1}{\overset{n_s}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}S{\hat{e}}_j=\underset{j=1}{\overset{n_s}{\mathrm{\Σ}}}\exp \left({\mathrm{i\φ}}_j\right)S\widehat{e_j}\\ {\mathrm{\φ}}_j\∈\left[\mathrm{0,2}\right]\mathrm{\π}\end{array}\right.,---\left(4\right)$

式中，表示“超级炮”，

表示单位向量，ns表示参与编码的震源个数，j＝1,2...,n_s表示参与编码的震源索引，α_j表示相位编码序列。

其中步骤7）中，所述迭代截止条件公式为：

$\left\{\begin{array}{c}\frac{{\mathrm{\ψ}}_k-{\mathrm{\ψ}}_{k-1}}{{\mathrm{\ψ}}_{k-1}}\≤\mathrm{\ϵ},\mathrm{\ϵ}>0\\ n_k\left(f_{\mathrm{ifreq}}\right)\≥N_{\max }\end{array}\right.,---\left(5\right)$

式中，ψ_k表示对应第k次迭代的目标函数，ε为预设的迭代停止阈值，n_k(f_ifreq)表示对应频率f_ifreq的迭代次数，N_max为预设的最大迭代次数。

本发明由于采取以上技术方案，其具有以下优点：本发明方法充分开发了在不同频率尺度下各种方法在全波形反演中的应用潜力，既大幅度地提高计算效率，又克服了各种方法的一些缺点，进而保证了全波形反演算法的稳定性。因此，在使用本发明方法进行地震速度建模时不仅能够得到高精度、高分辨率的速度模型，同时也满足工业生产中对计算成本的要求。这也是本发明方法能够代替或联合其他现有速度建模方法应用于地震资料处理中的必要条件和优势所在。本发明可以广泛用于深度域偏移成像、时深转换等。

附图说明

图1是本发明流程示意图

图2Marmousi模型应用实例-联合高效率全波形反演

图3Marmousi模型应用实例-炮集抽稀方法

图4Marmousi模型应用实例-震源编码方法

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明一种频率域全波形反演地震速度建模方法，包括以下步骤：

1）获取原始地震炮集记录、震源子波信息和反演所用的初始模型；其中初始模型是指一个粗略的初始速度模型；

2）分析由步骤1）获得的信息，基于现有的正演算法和优化算法，进而确定基本的反演参数和从低频到高频的全波形反演框架；其中正演算法采用基于混合网格的有限差分方法进行波场模拟，优化算法采用共轭梯度法；

3）对不同的频率，计算获得最合适的频域正反演模型网格；其中对不同频率，采用下面公式求其所对应的最适合的网格间距Δh。低频时网格间距较大，高频时网格间距较小。

$\left\{\begin{array}{c}\mathrm{\Δh}\≤\frac{v_{\min }}{M_{\min }{f}_{\mathrm{ifreq}}}\\ \mathrm{\Δh}=2^N\·\mathrm{\ΔS}\end{array}\right.,---\left(1\right)$

式中，v_min为介质最小速度，M_min为一个波长所包含网格的最少个数，由正演算法的精度控制，f_ifreq为反演对应频率，ifreq为反演频率的索引，Δs为炮间距，N为整数。这样，在低频反演时，由于采用的网格间距较大，因此在正反演中对应模型的网格维度较小，从而大幅度降低了在低频时的计算成本。

4）在低频反演时采用主成分分析方法，压缩参与反演的数据维度；主成分分析方法随着反演频率的增大，提速比逐渐降低。首先对频率域波场残差矩阵进行奇异值分解，然后通过积累能量分析和主能量成分提取确定一个投影矩阵，进而将原始数据投影至低维空间上，实现对原始数据的压缩；

$\left\{\begin{array}{c}\mathrm{\δD}={\mathrm{LXR}}^H\\ \frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_i^2}{\underset{i}{\mathrm{\Σ}}{\mathrm{\σ}}_i^2}\≥\mathrm{\χ}\\ D_{\mathrm{obs}}^k=D_{\mathrm{obs}}{D}^k,S^k={\mathrm{SR}}^k\end{array}\right.,---\left(2\right)$

式中，δD为波场残差矩阵，L和R分别表示左、右奇异值向量矩阵，R^H为R的共轭转置，X为对角奇异值（σ_i）矩阵，i＝1,2...,k为奇异值矩阵的维度索引，χ为控制投影矩阵维度k的阈值参数，R^k为投影矩阵，D_obs和S分别表示观测数据以及原始震源矩阵。

5）判断对应不同频率的投影矩阵维度是否满足下面的阈值转换标准，如果满足该转换标准，则进入步骤6）；否则回到步骤4）；

k(f_ifreq)＜μM,         （3）

其中，k(f_ifreq)表示对应频率f_ifreq得到的投影矩阵维度，M表示原始震源矩阵维度（列数），μ为预设的最小压缩比的倒数，即控制转换标准的阈值。

6）引入震源编码方法，并利用随机相位编码方法压制串扰噪声；随机相位震源编码方法是通过将多个独立震源编码组合成为一个“超级炮”，同时在迭代过程中随机生成编码序列来压制震源组合所产生的串扰噪声，震源编码的表达式为：；

$\left\{\begin{array}{c}\stackrel{\‾}{s}=\underset{j=1}{\overset{n_s}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}S{\hat{e}}_j=\underset{j=1}{\overset{n_s}{\mathrm{\Σ}}}\exp \left({\mathrm{i\φ}}_j\right)S\widehat{e_j}\\ {\mathrm{\φ}}_j\∈\left[\mathrm{0,2}\right]\mathrm{\π}\end{array}\right.,---\left(4\right)$

式中，

表示“超级炮”，

表示单位向量，n_s表示参与编码的震源个数，j＝1,2...,n_s表示参与编码的震源索引，α_j表示相位编码序列。

7）判定是否满足迭代截止条件，如果满足迭代截止条件，则进入下一步；否则回到步骤6）；其中迭代截止条件公式为：

$\left\{\begin{array}{c}\frac{{\mathrm{\ψ}}_k-{\mathrm{\ψ}}_{k-1}}{{\mathrm{\ψ}}_{k-1}}\≤\mathrm{\ϵ},\mathrm{\ϵ}>0\\ n_k\left(f_{\mathrm{ifreq}}\right)\≥N_{\max }\end{array}\right.,---\left(5\right)$

式中，ψ_k表示对应第k次迭代的目标函数，ε为预设的迭代停止阈值，n_k(f_ifreq)表示对应频率f_ifreq的迭代次数，N_max为预设的最大迭代次数。

8）如果没有完成所有频率的反演，则回到步骤3），直至完成所有频率的反演，获得最后的速度模型，输出速度模型。

为了更好的说明本发明，下面列举一实施例：

如图2所示，a为一个标准的二维Marmousi模型，试验中采用的原始网格大小为121*365，网格间距为25m，在1～21Hz范围内选择15个频率进行反演，每个频率迭代20次；b为反演采用的初始速度模型，c为采用常规全波形频率域方法得到的速度模型，d为采用本发明方法得到的速度模型。通过对比可以看出，采用本发明方法能够保证反演结果的精度要求，同时与常规全波形反演相比整体的计算效率提高了8.75倍，充分体现了本发明方法的优势。下面列出了对应不同频率所采用的具体方法以及对效率的提升比（如表1所示）

附表1本发明方法在Marmousi模型试验中的实现及效率分析

从表中可以看出，本发明方法在低频时所获得的效率提高明显优于高频时，这主要取决于变网格法的特性。

全波形反演是建立高精度速度模型的重要工具，采用频率域全波形反演相对于时间域全波形反演能够提高计算效率，但计算成本仍然较大。如何在频率域全波形反演的基础上进一步提高计算效率，降低计算成本是需要解决的问题。现有技术主要有炮集抽稀方法和震源编码方法。

这里仍然以二维Marmousi模型测试来说明，常规的频率域全波形反演方法虽然可以得到较高分辨率的速度模型（如图2中c所示），但其计算成本太大，限制了其应用。

炮集抽稀方法较为直接地提高了计算的效率，但由于反演中所用到的信息少，所以反演得到的速度模型的精度较低，如图3所示，a为抽取20%炮集数据；b为抽取40%炮集数据。对比可以看出，当抽取比例较小时，大量有效信息没有参与每次迭代反演中，因此反演的精度和分辨率受到严重的影响（如表2所示）：

附表2Marmousi模型应用各方法效率和精度数值分析

在震源编码方法中采用随机相位编码技术虽然可以在一定程度上压制串扰噪声，但是当参与编码的震源个数较多时（即获得较高的效率提升比），串扰噪声对反演结果的干扰仍较为明显，尤其是在低频时（如图4所示），每个超级炮由五个独立震源编码生成，每个频率迭代a为20次；b为80次。从图中可以看出，使用随机相位震源编码方法时，虽然可以通过增加迭代次数压制串扰噪声的影响，获得较为精确的反演模型，但是这样却增加了计算量（如表2所示），降低了该方法的效率提升比。

另外，对于不同频率，炮集抽稀方法和震源编码方法对效率的提升比均是保持一致的。也就是说，这两种方法并没有体现和发挥出频率域全波形反演在不同频率尺度下的特点和潜力。本发明方法采用与上述两种提高效率的方法不同，本发明的频率域全波形反演方法是在频率域反演框架下，利用变网格法、主成分分析和震源编码等三种方法对应不同频率所表现的优缺点，动态地选择或结合各个方法，在保证了速度建模精度的同时大幅度地提高了计算效率。

本发明方法同抽稀炮方法和震源编码方法相对比（如图2、图3、图4所示），可以发现：首先，炮集抽稀方法所得速度模型分辨率明显较低；其次，震源编码方法反演所得速度模型受串扰噪声干扰较严重，且界面连续性较差；最后，本发明方法充分利用了变网格法等三种方法在不同频率尺度下所表现的特点，相比其他两种方法，本发明所建速度模型更加接近真实速度模型，精度较高。

上述表2列出了本发明方法同炮集抽稀方法、震源编码方法的效率和反演精度对比。其中，迭代次数表示每个频率执行的迭代次数，效率提升比是指各方法相对于常规方法节省耗时的倍数，反演精度是描述反演后模型同实际模型的相似程度。对比可以看出，首先炮集抽稀方法受限于反演精度的要求，所能提供的效率提升比较低；其次，震源编码方法在参与编码震源较多时受串扰噪声影响较大，需要较多的迭代次数对此进行压制，从而对效率的提升亦有限；最后，本发明方法同其他两种方法不同，它有机地结合变网格法等三种方法对效率的多层次提高，对效率的提升比明显优于其他两种方法。

上述各实施例仅用于说明本发明，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (8)

1.一种频率域全波形反演地震速度建模方法，包括以下步骤： 

1）获取原始地震炮集记录、震源子波信息和反演所用的初始模型； 

2）分析步骤1）获得的信息，基于正演算法和优化算法，确定基本的反演参数和从低频到高频的全波形反演框架； 

3）对不同的频率，计算获得最合适的频域正反演模型网格； 

4）在低频反演时采用主成分分析方法，压缩参与反演的数据维度； 

5）判断对应不同频率的投影矩阵维度是否满足阈值转换标准，当满足该转换标准时，进入下一步，否则回到步骤4）； 

6）引入震源编码方法，并利用随机相位编码方法压制串扰噪声； 

7）判定是否满足迭代截止条件，如果满足迭代截止条件，则进入下一步；否则回到步骤6）； 

8）如果没有完成所有频率的反演，则回到步骤3），直至完成所有频率的反演，获得最后的速度模型，输出速度模型。 

2.如权利要求1所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤2）中，所述正演算法采用基于混合网格的有限差分方法，所述优化算法采用共轭梯度法。 

3.如权利要求1所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤3）中，对不同的频率采用下面公式求其所对应的最适网格间距Δh： 

式中，v_min为介质最小速度，M_min为一个波长所包含网格的最少个数，由正演算法的精度控制，f_ifreq为反演对应频率，ifreq为反演频率的索引，Δs为炮间距，N为整数。 

4.如权利要求2所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤3）中，对不同的频率采用下面公式求其所对应的最适网格间距Δh： 

式中，v_min为介质最小速度，M_min为一个波长所包含网格的最少个数，由正演算法的精度控制，f_ifreq为反演对应频率，ifreq为反演频率的索引，Δs为炮间距，N为整数。 

5.如权利要求1或2或3或4所述的一种频率域全波形反演地震速度建模方法，其特征在于：所述步骤4）中，在使用主成分分析时，首先对频率域波场残差矩阵进 行奇异值分解，然后通过积累能量分析和主能量成分提取确定一个投影矩阵，进而将原始数据投影至低维空间上，实现对原始数据的压缩；主成分分析方法的表达式为： 

式中，δD为波场残差矩阵，L和R分别表示左、右奇异值向量矩阵，R^H为R的共轭转置，X为对角奇异值（σ_i）矩阵，i为奇异值矩阵的维度索引，χ为控制投影矩阵维度k的阈值参数，R^k为投影矩阵，D_obs和S分别表示观测数据以及原始震源矩阵。 

6.如权利要求1或2或3或4或5所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤5）中，所述阈值转换标准的表达式为： 

k(f_ifreq)＜μM,         （3） 

式中，k(f_ifreq)表示对应频率f_ifreq得到的投影矩阵维度，M表示原始震源矩阵维度，μ为预设的最小压缩比的倒数，即控制转换标准的阈值。 

7.如权利要求1或2或3或4或5所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤6中，所述随机相位震源编码方法是通过将多个独立震源编码组合成为一个“超级炮”，同时在迭代过程中随机生成编码序列来压制震源组合所产生的串扰噪声，震源编码的表达式为： 

式中，

表示“超级炮”，表示单位向量，n_s表示参与编码的震源个数，j＝1,2...,n_s表示参与编码的震源索引，α_j表示相位编码序列。 

8.如权利要求1或2或3或4或5或6所述的一种频率域全波形反演地震速度建模方法，其特征在于：其中步骤7）中，所述迭代截止条件公式为： 

式中，ψ_k表示对应第k次迭代的目标函数，ε为预设的迭代停止阈值，n_k(f_ifreq)表示对应频率f_ifreq的迭代次数，N_max为预设的最大迭代次数。 

Images