反演速度模型及其建立方法和地下结构的像的获得方法与流程

本发明涉及地球物理勘探领域，特别是油气勘探领域的反演速度模型及其建立方法和地下结构的像的获得方法。

背景技术：

为在油气勘探领域找到大油田，需要对地下结构精确成像。为获得地下结构的像，首先需要知道地下结构的速度信息。自从tarantola于1984年提出基于广义最小二乘反演理论的时间域全波形反演方法(以基于地震波的位移和压力建立的反演速度模型作为全波形反演的初始速度模型)以来，地球物理学研究者似乎认为找到了一种可以高精度重建地下速度的方法。几十年来，地球物理学研究者发表了上千篇相关文章，主要探讨如何应用地震全波形获得地下结构信息。伴随着计算机集群计算能力的大幅度提高，更多的全波形反演方法被提出，其中最有影响力的方法之一是pratt(1995)提出的基于多尺度的频率域全波形反演方法，该方法从低频信息开始迭代，逐渐走向高频，最终获得高精度的地下速度结构。基于多尺度的频率域全波形反演方法通过将问题分解为不同尺度求解以保证反演过程稳定收敛。

不过，目前所有的研究者都由于遇到以下两个问题而使关于如何高精度重建地下速度的研究难以继续前进：其中之一，是观测的地震数据中没有低频信息，如在油气勘探领域的地震数据的主频一般在15hz以上，而batean于2013年证明，地震数据中须含有低于2hz以下的频率数据，才能够通过波形反演获得高精度地下速度结构。其中之二是，在全波形反演的初始速度模型合成的地震记录与地震观测记录的差大于半个波长的情形下，反演迭代将不会收敛，因此为确保反演迭代过程的收敛，要求初始速度模型必须与真实速度模型很接近。

因此，本发明需要一种新的方案来解决上述问题。

技术实现要素：

技术问题

有鉴于此，本发明要解决的技术问题是，如何获得与真实速度模型更为接近的、作为地震波全波形反演的初始速度模型的反演速度模型。

解决方案

为了解决上述技术问题，根据本发明的一个实施例，提供了一种反演速度模型的建立方法，该方法包括：

将地震观测记录转换为第一中间数据，且该第一中间数据中包含频率不高于第一预设频率的数据；

根据设定的初始速度模型计算得出地震合成记录，并将地震合成记录转换为第二中间数据，且该第二中间数据中包含频率不高于第一预设频率的数据；

基于所述第一中间数据和所述第二中间数据，建立目标函数；

根据所述目标函数，获取地震波的速度的梯度δ(m)；

计算δ(m)的迭代步长α；

根据所述δ(m)以及所述α，确定反演速度模型。

在上述建立方法的优选技术方案中，所述第一中间数据和所述第二中间数据均为能够反映地震波强度的数据。

在上述建立方法的优选技术方案中，所述第一中间数据为地震观测记录的平方，所述第二中间数据为地震合成记录的平方。

在上述建立方法的优选技术方案中，“基于所述第一中间数据和所述第二中间数据，建立目标函数”的步骤进一步包括：对所述第一中间数据和所述第二中间数据分别进行傅里叶变换，保留所述第一中间数据和所述第二中间数据中频率不高于第一预设频率的部分，基于所述不高于第一预设频率的部分建立目标函数。

在上述建立方法的优选技术方案中，所述目标函数具体为：

公式(1)中，m为地震波的速度，p0是地震观测记录，是p0的复共轭，是第一中间数据，pc是地震合成记录，pc^*是pc的复共轭，是第二中间数据。

在上述建立方法的优选技术方案中，“根据所述目标函数，计算地震波的速度的梯度δ(m)”的步骤进一步包括：

计算目标函数对地震波的速度的导数，应用泰勒展开，求取速度的梯度δ(m)；

公式(2)中，δp²＝(pcpc^*-p0p0^*)，是pc对m的一阶导数。

在上述建立方法的优选技术方案中，“根据所述δ(m)以及所述α，确定反演速度模型”的步骤进一步包括：

计算并比较本次迭代和上次迭代的目标函数值：

如果本次迭代的目标函数值大于或等于上次迭代的目标函数值，则迭代终止，根据上次迭代的结果确定反演速度模型；

如果本次迭代的目标函数值小于上次迭代的目标函数值，则判断本次迭代的目标函数的值是否小于等于预设的目标函数的收敛误差β；

若本次迭代的目标函数的值是小于等于预设的目标函数的收敛误差β，则迭代终止，根据本次迭代的结果确定反演速度模型；

若本次迭代的目标函数的值大于预设的目标函数的收敛误差β，则通过以下的公式(3)确定出更新后的初始速度模型，继续迭代，直到目标函数的值小于收敛误差β；

m＝m0+αδ(m)(3)

公式(3)中，m0为地震波的初速度。

在上述建立方法的优选技术方案中，所述第一预设频率为2-5hz的量。

为了解决上述技术问题，根据本发明的另一个实施例，提供了一种反演速度模型，该反演速度模型是根据前述的建立方法构建。

为了解决上述技术问题，根据本发明的第三个实施例，提供了一种地下结构的像的获得方法，该方法包括：

以前述的反演速度模型作为初始速度模型，采用地震波全波形反演方法获得地下速度模型；

基于该地下速度模型，获得地下结构的速度；

根据所述地下结构的速度，获得地下结构的像。

有益效果

本发明的地震波强度全波形反演(fullintensitywaveforminversion，fiwi)方法在建立反演速度模型的过程中，通过将很难获得低频信息的地震数据直接转换为包含低频信息(尤其是极低频信息)的如能够反映地震波强度的中间数据，并基于中间数据构建目标函数，从低频开始迭代，获得对设定的初始速度模型修正后的反演速度模型。然后将该反演速度模型作为传统的地震波全波形反演(fullwaveforminversion，fwi)方法的初始速度模型，通过一系列的迭代获得更高精度的地下速度结构。由于中间数据中包含低频信息，因此反演速度模型的建立不再要求地震数据中含有低频信息；由于基于低频信息的反演迭代更容易收敛，因此不再要求设定的初始反演速度模型与真实速度模型非常接近。

根据下面参考附图对示例性实施例的详细说明，本发明的其它特征及方面将变得清楚。

附图说明

包含在说明书中并且构成说明书的一部分的附图与说明书一起示出了本发明的示例性实施例、特征和方面，并且用于解释本发明的原理。

图1示出本发明一个实施例的反演速度模型的建立方法的流程示意图。

图2示出波场反传播的原理示意图。

图3示出本发明的整个反演流程的逻辑框图。

图4a示出第一种地震记录(互井数据)中的单道数据的示意图；图4b示出该单道数据的频谱示意图；图4c示出该单道数据的地震波强度数据的示意图；图4d示出该地震波强度数据的频谱示意图。

图5a示出真实速度模型(与互井数据对应)的示意图。图5b示出以5490米/秒的常速度模型作为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图；图5c示出以5490米/秒的常速度模型作为初始速度模型，采用本发明的fiwi方法迭代后的反演结果示意图；图5d示出将图5c的反演结果作为初始速度模型，采用传统的fwi方法迭代100次后的示意图。

仍为图5a所示的真实速度模型，图6a示出以5000米/秒的常速度模型作为初始速度模型，采用传统的fwi方法进行多尺度迭代后的结果示意图；图6b示出以5000米/秒的常速度模型作为初始速度模型，采用本发明的fiwi方法进行多尺度迭代后的反演结果示意图；图6c示出以图6b所示的反演结果作为初始速度模型，采用传统的fwi方法迭代30次后的结果示意图；图6d示出以图6b所示的结果作为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图。

图7a示出第二种地震记录(地表数据)中的单道数据的波形示意图，图7b示出该单道数据的频谱示意图；图7c示出将该单道数据的波形转换为地震波强度数据的结果示意图；图7d示出该地震波强度数据的频谱示意图。

图8a示出第二种真实速度模型(与地表数据对应)的示意图；图8b是反演采用的初始速度模型的示意图；图8c示出以图8b作为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图；图8d示出以图8b作为初始速度模型，采用本发明的fiwi方法进行迭代后的反演结果示意图；图8e示出将图8d的反演结果作为初始速度模型，采用传统的fwi方法迭代100次后的示意图。

图9a示出采用传统的fwi方法时不同的速度模型的对比图；图9b示出采用本发明的fiwi方法时不同的速度模型的对比图。

具体实施方式

以下将参考附图详细说明本发明的各种示例性实施例、特征和方面。附图中相同的附图标记表示功能相同或相似的元件。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。

另外，为了更好地说明本发明，在下文的具体实施方式中给出了众多的具体细节。本领域技术人员应当理解，没有某些具体细节，本发明同样可以实施。在一些实例中，对于本领域技术人员熟知的方法和手段未作详细描述，以便于凸显本发明的主旨。

图1示出本发明的一种反演速度模型的建立方法的流程示意图。如图1所示，该反演速度模型的建立方法主要包括：

s100、将地震观测记录转换为包含频率不高于第一预设频率(低频信息)的第一中间数据；其中，第一预设频率可以为不大于2-5hz的量，具体的设定值可以根据实际需要灵活调整。如本实施例中以第一预设频率为2hz为例来阐述本发明的原理。即低频信息为不高于2hz的信息。

s200、根据设定的初始速度模型计算得出地震合成记录，并将地震合成记录转换为包含频率不高于第一预设频率的低频信息的第二中间数据；

s300、基于第一中间数据和第一中间数据，建立目标函数；

s400、计算地震波的速度的迭代步长；以及根据目标函数，获取地震波的速度的梯度；

s500、根据地震波的初速度、地震波的速度的梯度以及地震波的速度的迭代步长，确定反演速度模型。

可以看出，通过将不包含低频信息的地震数据转换为必然包含低频信息的中间数据，使反演速度模型的建立不再要求地震数据中含有低频信息。而且由于基于低频信息的反演迭代更容易收敛，因此不再要求作为反演速度模型与真实速度模型非常接近。

地震波的强度的定义为在单位时间和单位面积通过地震波能流垂直面的能量(bornandwolf，1999)，地震波的强度信息能够反映地震波在地下能量的传播情况，而含有低频甚至极低频信息是地震波的强度信息中的固有特征。鉴于此，与目前直接利用地震数据建立反演速度模型不同，本发明在建立反演速度模型的过程中，首先将作为地震数据的地震观测记录和地震合成记录分别转换成与地震波的强度相关(能够反映地震波的强度)的第一中间数据和第二中间数据。由于地震波的强度信息中必然含有低频部分和高频部分，因此第一中间数据和第二中间数据也必然含有低频部分和高频部分。在去掉第一中间数据和第二中间数据中的高频部分之后，用低频部分构建目标函数，从低频开始迭代，获得一个反演速度模型。将该反演速度模型作为fwi方法初始速度模型，经过一系列的迭代获得更高精度的地下速度结构。由于转换为中间数据之后必然含有低频成分，因此本发明不再要求地震数据本身必须含有低频信息，这样反演迭代显然更容易收敛，因此本发明不再要求初始反演速度模型与真实速度模型十分接近。

其中，地震波的强度信息是与地震数据的平方相关，应当理解为二者的变化趋势大致相同，如可以用一个系数来实现地震波强度信息与地震数据之间的转换。具体而言，地震波的强度信息是正比于地震数据的平方的量，因此可以直接将地震波的强度信息或者与地震波的强度信息的趋势等同的量(如地震数据的平方)作为中间数据来建立目标函数。下面以地震观测记录的平方和地震合成记录的平方分别作为第一中间数据(以下称为第一地震波强度数据)和第二中间数据(以下称为第二地震波强度数据)为例，来阐述本发明的原理和建立反演速度模型的过程。

反演速度模型的建立过程具体包括如下步骤：

s1)将所有炮的地震观测记录的平方作为第一地震波强度数据。对第一地震波强度数据进行傅里叶变换，获得第一地震波强度数据的分频信息，即将第一地震波强度数据分成高频部分与低频部分，去掉其中的高频部分，仅保留低频部分。再通过逆傅里叶变换，回到地震数据的时间域，即可获得只含有低频信息的地震波强度数据。

s2)设置一个初始速度模型(可以是任意模型，如常速度模型或者线性速度模型)，以butterworth或者ricker为震源子波，采用有限差分计算出地震波动方程，具体如公式(1)所示：

其中，(x,y,z)为任意网格点的地震合成记录对应的地下三维坐标，p为地震波的波场，m0为地震波的初速度，t为地震波的移动时间。

根据地震波动方程(1)进一步建立地震合成记录，并将地震合成记录的平方作为第二地震波强度数据。对第二地震波强度数据进行傅里叶变换，获得第二地震波强度数据的分频信息，即将第二地震波强度数据分成低频部分与高频部分，去掉其中的高频部分，仅保留低频部分。再通过逆傅里叶变换，回到地震数据的时间域，即可获得只含有低频信息的地震波强度数据。

s3)建立第一地震波强度数据与第二地震波强度数据的2范数的平方的目标函数，具体如公式(2)所示：

其中，m为地震波的速度，p0是地震观测记录，是p0的复共轭，是作为第一中间数据的第一地震波强度数据；pc是地震合成记录，即地震波动方程(1)的解，pc^*是pc的复共轭，是作为第二中间数据的第二地震波强度数据。

s4)由于不能保证设置出的初始速度模型与真实速度模型接近，因此需要更新初始速度模型。更新模型的参数包括地震波的速度的梯度δ(m)和δ(m)的迭代步长α；其中：

s41)基于(第一、第二)地震波强度数据，应用泰勒展开，计算目标函数对速度的导数，获得地震波的速度的梯度方程(3)。从低频开始迭代，求取要更新的地震波的速度的梯度δ(m)；

其中，是pc对速度m的梯度，也叫frechét导数。其中，δ(m)的计算过程具体为：

s411)将地震合成记录的地震波场乘以第二地震波强度数据与第一地震波强度数据的差的积(即pcδp²)作为复共轭震源；

s412)图2示出波场反传的原理示意图。如图2所示，对从s点正向传播的震源，通过双程波波动方程计算出在地下任何一个位置的正向传播的地震波场的二阶时间导数ps；

s413)继续参照图2，对从r点反向传播的共轭震源pcδp²，通过双程波波动方程，计算出在地下任何一个位置的反向传播的地震波场pr；

s414)在地下任何一个网格点，对ps与pr作互相关，进而参照公式(3)计算出m的梯度δ(m)。

s42)通过抛物线法搜索δ(m)的迭代步长α；

s5)设定一个迭代收敛误差，经过一系列的迭代，当δ(m)的修正量小于该设定的迭代收敛误差，则终止迭代，并根据地震波的初速度m0、m的梯度δ(m)和δ(m)的迭代步长α获得一个δ(m)修改的反演速度模型，即：

m＝m0+αδ(m)(4)

参照公式(2)示出的目标函数可以看出，在初始速度模型与真实速度模型完全吻合的情形下，目标函数的值为0，且目标函数的最小值为0。因此，在理想状态下，若终止迭代时的目标函数取值为0，则通过fiwi方法获得的反演速度模型恰好为真实速度模型。

以s5)获得的反演速度模型作为fwi方法的初始速度模型，经过一系列迭代(如80次、100次等)，获得最终的地下速度模型。根据该地下速度模型，获得地下结构的速度。根据地下结构的速度，最终获得地下结构的像。基于该地下结构的像，即可在油气勘探领域确定大油田地下结构，或者进行关于地球全球构造的研究。

可以看出，由于作为中间数据的地震波强度数据中必然包含有低频信息，因此不再要求地震数据中包含低频信息。以地震波强度数据为基础建立目标函数，最终速度的梯度修正后的反演速度模型作为fwi方法的初始速度模型来重建真实速度模型，因此本发明不再要求初始速度模型与真实速度模型十分接近。

图3示出整个反演流程的逻辑框图。如图3所示，整个反演流程主要包括：1)通过初始模型(即设定的初始速度模型，如常速度模型或者线性模型)正演计算出正演数据，即模拟得到地震合成记录；2)将观测数据(即地震观测记录)和正演数据分别转换为能够数据强度(即能够反映地震波强度的量)。具体而言，将观测数据和正演数据分别转换为第一地震波强度数据和第二地震波强度数据；3)对数据强度进行多尺度低通滤波，从低频到高频，依次建立对应的目标函数。具体而言，将(第一、第二)地震波强度数据分解成不同频率的尺度，应用多尺度反演方法计算目标函数对于速度的修正量：如果目标函数的值不小于上一次迭代的值，则终止fiwi方法的迭代流程；如果目标函数下降(即目标函数的值小于上一次迭代的值)，则计算搜索步长和更新的速度的梯度(速度修正量)，并基于搜索步长和速度的梯度更新初始模型，以更新后的初始模型再次正演计算出新的正演数据，以观测数据和新的正演数据为基础进一步建立目标函数，以此类推，经过多次迭代后，直到目标函数值小于预设的收敛误差β(如直至目标函数不再下降)，即终止fiwi方法的迭代流程。迭代流程终止时获得的反演速度模型作为传统的fwi方法的初始速度模型，经过一系列的反演迭代，最终获得最终的高精度地下速度模型。可以看出，与传统的依据地震数据(如地震波的位移或者压力场)建立目标函数的方法完全不同，本发明是应用由地震数据直接得出的中间数据(如地震波强度数据)来建立目标函数。仍以中间数据为地震波强度数据为例，由于地震波强度数据中包含有低频信息，因此不再需要地震数据含有低频信息，并且可以克服地震反演过程中存在的周波跳跃问题。

下面通过两个实施例来验证本发明的fiwi方法的可行性和有效性。

实施例1

图5a示出第一种真实速度模型(与互井数据对应)的示意图，图5b示出以5490米/秒的常速度模型作为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图，图5c示出以5490米/秒的常速度模型作为初始速度模型，采用本发明的fiwi方法进行迭代后的反演结果示意图(修正后的初始速度模型)，图5d示出以图5c的反演结果作为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图。如图5a所示，本发明的第一种真实速度模型是一个互井模型，该模型在水平方向(distance)的距离为120米，在深度方向(depth)的距离为300米，水平方向与深度方向之间的网格间距为1米。在模型的左侧布置76炮，炮间距为3米，起始炮点为(x＝10米，z＝50米)，在模型的右侧布置76个检波器，道间距为3米，起始道为(x＝110米，z＝50米)。模型的速度变化范围为4390-6160米/秒。利用空间四阶、时间二阶有限差分计算地震合成记录。

为了进一步限制地震合成记录的数据频带，采用butterworth作为震源子波，震源子波的频带为150-200-400-450hz。图4a示出地震合成记录中(互井数据)的单道数据(即单道井间数据)的示意图，图4b示出单道井间数据的频谱(振幅谱)示意图。需要说明的是，由于ricker子波在低频段存在渐变衰减的缺陷，而butterworth子波不存在低频信息，因此以butterworth作为震源子波能够使振幅谱在通带范围内尽量保持平稳。图4c示出单道井间数据的地震波强度数据(地震波的振幅的平方)的示意图，图4d示出单道井间数据的地震波强度数据的频谱(振幅谱)的示意图。对比图4b和图4d可以看出，相比地震数据，地震波强度数据包含更多的低频信息，而且地震波强度数据明显地将频谱分为低频部分和高频部分。这样一来，只需要去掉其中的高频部分、利用保留的低频部分即可实现地震波强度的全波形反演。

下面采用两个不同大小的初始常速度模型作为设定的初始速度模型来测试本发明根据地震波强度数据建立的反演速度模型对初始速度模型的依赖。

测试例11：

以5490米/秒的常速度模型作为设定的初始速度模型的情形下，参照图5a可以看出，与垂直方向(velocity)给出的调色棒对比，真实速度模型的速度绝大部分都分布在5490米/秒附近或者高于此速度。也就是说，由于该初始速度模型相比真实速度模型的差异较小，因此采用传统的fwi方法迭代100次后，迭代收敛，迭代结果如图5b所示。仍以5490米/秒的常速度模型作为初始速度模型，采用本发明的fiwi方法进行多尺度反演迭代，在反演迭代的过程中，用butterworth滤波器对地震波强度数据进行分频(分为3个频段)，并从低频到高频进行迭代(0-0-50-100hz迭代1次，0-0-100-150hz迭代5次，0-0-150-200hz迭代1次)，迭代结果如图5c所示。进一步以fiwi方法获得的反演速度模型作为初始速度模型，用传统的fwi方法迭代100次，结果如图5d所示。

结合图5a、5b和5d可以看出，在初始速度模型为5490米/秒的情形下，通过传统的fwi方法和本发明的fiwi方法都可以获得一个较好的速度反演结果。也就是说，在初始速度模型与真实速度模型十分接近的情形下，采用传统的fwi方法和本发明的fiwi方法均能够使真实速度模型的细节得到较好的重建。

测试例12：

以5000米/秒的常速度模型作为设定的初始速度模型的情形下，继续参照图5a可以看出，与垂直方向(velocity)给出的调色棒对比，真实速度模型的速度只有很少的部分接近5000米/秒。也就是说，由于该设定的初始速度模型相比真实速度模型的差异较大，因此反演过程存在更加严重的周波跳跃(cycle-skipping)问题。具体地，采用传统的fwi方法迭代了9次之后便不再收敛，迭代结果如图6a所示。也就是说，在以5000米/秒的常速度模型作为设定的初始速度模型的情形下，采用传统的fwi方法不能重建真实速度模型的真实速度。仍以5000米/秒的常速度模型作为设定的初始速度模型，首先采用本发明的fiwi方法进行多尺度反演迭代，在反演迭代的过程中，用butterworth滤波器将地震波强度数据进行分频(分为3个频段)，并从低频到高频进行迭代(其中，0-0-20-50hz迭代4次，0-0-40-80hz迭代2次，0-0-60-120hz迭代3次)，反演结果如图6b所示。以此反演结果作为初始速度模型，采用传统的fwi方法迭代30次和100次的结果分别如图6c和图6d所示。可以看出，真实速度模型的细节在两种不同迭代次数的情形下均能够得到较好的重建。

可以看出，在设定的初始速度模型与真实速度模型的误差较大的情形下，通过传统的fwi方法由于存在严重的周波跳跃问题而使迭代不收敛，从而导致无法重建出竖井模型的真实速度。而本发明的fiwi方法则可以较大程度地克服周波跳跃问题，因此即使在设定的初始速度模型与真实速度模型的误差较大的情形下，采用本发明的fiwi方法仍然能够获得一个较好的速度反演结果，即仍然能够较好地重建出竖井模型的真实速度。

对比上述测试例11(以5490米/秒的常速度模型作为设定的初始速度模型)和测试例12(以5000米/秒的常速度模型作为设定的初始速度模型)可以看出，与传统的fwi方法相比，本发明的fiwi方法在保证能够较好地重建出真实速度模型的前提下，对设定的初始速度模型的依赖更低，即不再要求设定的初始速度模型与真实速度模型十分接近。

实施例2

图8a示出第二种真实速度模型(与地表数据对应)的示意图，图8b是设定的初始速度模型的示意图，图8c示出以图8b为初始速度模型，采用传统的fwi方法迭代100次后的结果示意图；图8d示出以图8b为初始速度模型，采用本发明的fiwi方法进行迭代后的反演结果示意图(反演速度模型)；图8e示出将图8d的反演结果作为初始速度模型，采用传统的fwi方法迭代100次后的示意图。

如图8a所示，本发明的第二种真实速度模型是2d-marmousi模型，该模型是从原始的marmousi模型重新采样得到的。将该模型在水平方向(distance)上离散化为661个网格，在垂直方向(depth)上离散化为201个网格，网格间距在水平方向与垂直方向上均为10米。如图8b所示，反演所采用的设定的初始速度模型线性速度模型，具体而言，该设定的初始速度模型为速度从1500米/秒到3500米/秒线性增加的模型。地震波动方程正演有111炮，炮间距为60米，均匀分布在地表。每炮对应661个检波器，检波器之间的间距为10米，地震子波是主频为10hz的ricker子波。

图7a示出第二种地震数据(地表数据)中的单道数据的波形示意图，将地震数据采用5-7-9-12hz的butterworth滤波器进行滤波，得到如图7b所示的单道数据的频谱(振幅谱)示意图。可以看到，频谱中没有5hz以下的低频信息。将地震数据转换为地震波强度数据的结果如图7c所示，地震波强度数据的频谱(振幅谱)如图7d所示。可以看到，相比于图7b中的地震数据，地震波强度数据的频谱中含有非常低的频率。

下面通过以下两种反演流程来验证fiwi方法在重建真实速度模型的过程中的作用：

第一种反演流程：基于设定的初始速度模型，采用传统的fwi方法进行反演迭代；

第二种反演流程：基于设定的初始速度模型，采用fiwi方法进行反演迭代，获得反演速度模型；以该反演速度模型作为初始速度模型，采用传统的fwi方法进行反演迭代(以下称联合fiwi和fwi的方法)。

测试例21：

在第一种反演流程中，首先以图8b所示的线性速度模型作为初始速度模型，采用传统的fwi方法进行多尺度反演。具体地，依次按照5-7hz、5-9hz、5-11hz和5-12hz总共迭代100次，最终得到的反演结果如图8c所示。

可以看出，虽然图8c的反演结果提供了一定的精细结构，但与图8a示出的真实的2d-marmousi模型的总体趋势不一致，特别是与较深部分(depth较大的部分)的高速度不匹配。这主要是由于地震数据缺乏频率低于5hz的低频信息而导致初始速度模型与真实速度模型在较深部分具有较大误差，从而使反演陷入了局部极值。

测试例22：

在第二种反演流程中，首先将地震波强度数据采用butterworth滤波器进行滤波，仅保留0-2hz的低频部分，然后以图8b所示的线性速度模型作为初始速度模型，采用传统的fwi方法进行多尺度反演迭代，迭代95次后的反演结果如图8d所示。

可以看出，由于地震波强度数据中含有低频信息，因此采用本发明的fiwi方法进行多尺度反演时会产生一个低分辨率但具有正确的长波长分量的中间模型(参照8d)。然后以该中间模型作为初始速度模型，采用传统的fwi方法进行反演迭代，反演过程与第一种反演流程相同，得到的反演结果如图8e所示。可以看出，与传统的fwi方法相比较，采用联合fiwi和fwi的方法的反演结果与真实速度模型更加吻合，即能够更为准确地重建真实速度模型。

进一步地，图9a和图9b示出了在采用fwi和fiwi方法时实施例2中的不同速度模型的对比图。其中，图9a和图9b是在相同的水平位置处(3.3km)抽取的两个垂直剖面，如图9a所示，实线为真实速度模型，虚线为初始速度模型，点线为采用传统的fwi方法进行多尺度反演的结果。如图9b所示，实线为真实速度模型，虚线为采用本发明的fiwi方法得到的中间模型，点线为采用联合fiwi和fwi的方法得到的最终反演结果。

结合图8a-8e和图9a-9b可以看出，联合fiwi和fwi的方法可以得到较好的反演结果，即能够较好地重建真实速度模型。具体地，进一步参照图9b，通过使用更精确的中间模型作为fwi方法的初始速度模型，最终为真实速度模型的浅层和较深部分均提供了更加均衡的反演速度。而且，采用接近真实速度模型的初始速度模型可以较好地反演出真实速度模型的精细结构。举例而言，真实速度模型(参照图8a)中的顶部1.0千米处存在分层结构，联合fiwi和fwi的方法(参照图8e、图9a和图9b)适当地再现了这些层的结构，而传统的fwi方法(参照图8c)却未能产生这些浅层的结构。

可以看出，在本发明的联合fiwi和fwi的方法中，由于地震波的强度场包含有低频信息，因此地震波强度全波形反演方法可以从超低频数据开始，通过应用多尺度递归迭代为地震波全波形反演准备一个较好的初始速度模型。这样一来，联合fiwi和fwi的方法对反演开始时设定的初始速度模型没有要求，如初始速度模型可以是常速度或者线性渐进速度，而且本发明的联合fiwi和fwi的方法克服了传统的fwi方法在地震合成记录与地震观测记录超过半个周期时无法获得收敛解的问题。如与传统的fwi方法相比，本发明在“地震数据中无低频信息、设定的初始速度模型为与真实速度模型差异较大的常速度模型”的情形下仍然能够实现反演迭代的收敛，即仍然能够较好地重建真实速度模型。

很明显地，通过本发明中联合fiwi和fwi的反演方案能够将地震波的反演研究向前推进一大步，这样的反演研究能够被广泛地应用于油气勘探领域中。

需要说明的是，尽管以包含低频信息的中间数据为地震数据的平方作为示例介绍了联合fiwi和fwi的全波形反演方案如上，但本领域技术人员能够理解，本发明应不限于此。事实上，完全可根据实际情况，灵活地调整中间数据的具体形式，如对地震数据进行处理之后再求平方，或者平方之后再赋予相应的系数等，只要满足由地震数据直接得出的中间数据由于与地震波强度数据相关因此必然包含低频信息的条件即可。

至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。