反演速度模型及其建立方法和地下結構的像的獲得方法與流程

2023-06-25 07:20:51 1

本發明涉及地球物理勘探領域，特別是油氣勘探領域的反演速度模型及其建立方法和地下結構的像的獲得方法。

背景技術：

為在油氣勘探領域找到大油田，需要對地下結構精確成像。為獲得地下結構的像，首先需要知道地下結構的速度信息。自從tarantola於1984年提出基於廣義最小二乘反演理論的時間域全波形反演方法(以基於地震波的位移和壓力建立的反演速度模型作為全波形反演的初始速度模型)以來，地球物理學研究者似乎認為找到了一種可以高精度重建地下速度的方法。幾十年來，地球物理學研究者發表了上千篇相關文章，主要探討如何應用地震全波形獲得地下結構信息。伴隨著計算機集群計算能力的大幅度提高，更多的全波形反演方法被提出，其中最有影響力的方法之一是pratt(1995)提出的基於多尺度的頻率域全波形反演方法，該方法從低頻信息開始迭代，逐漸走向高頻，最終獲得高精度的地下速度結構。基於多尺度的頻率域全波形反演方法通過將問題分解為不同尺度求解以保證反演過程穩定收斂。

不過，目前所有的研究者都由於遇到以下兩個問題而使關於如何高精度重建地下速度的研究難以繼續前進：其中之一，是觀測的地震數據中沒有低頻信息，如在油氣勘探領域的地震數據的主頻一般在15hz以上，而batean於2013年證明，地震數據中須含有低於2hz以下的頻率數據，才能夠通過波形反演獲得高精度地下速度結構。其中之二是，在全波形反演的初始速度模型合成的地震記錄與地震觀測記錄的差大於半個波長的情形下，反演迭代將不會收斂，因此為確保反演迭代過程的收斂，要求初始速度模型必須與真實速度模型很接近。

因此，本發明需要一種新的方案來解決上述問題。

技術實現要素：

技術問題

有鑑於此，本發明要解決的技術問題是，如何獲得與真實速度模型更為接近的、作為地震波全波形反演的初始速度模型的反演速度模型。

解決方案

為了解決上述技術問題，根據本發明的一個實施例，提供了一種反演速度模型的建立方法，該方法包括：

將地震觀測記錄轉換為第一中間數據，且該第一中間數據中包含頻率不高於第一預設頻率的數據；

根據設定的初始速度模型計算得出地震合成記錄，並將地震合成記錄轉換為第二中間數據，且該第二中間數據中包含頻率不高於第一預設頻率的數據；

基於所述第一中間數據和所述第二中間數據，建立目標函數；

根據所述目標函數，獲取地震波的速度的梯度δ(m)；

計算δ(m)的迭代步長α；

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

在上述建立方法的優選技術方案中，所述第一中間數據和所述第二中間數據均為能夠反映地震波強度的數據。

在上述建立方法的優選技術方案中，所述第一中間數據為地震觀測記錄的平方，所述第二中間數據為地震合成記錄的平方。

在上述建立方法的優選技術方案中，「基於所述第一中間數據和所述第二中間數據，建立目標函數」的步驟進一步包括：對所述第一中間數據和所述第二中間數據分別進行傅立葉變換，保留所述第一中間數據和所述第二中間數據中頻率不高於第一預設頻率的部分，基於所述不高於第一預設頻率的部分建立目標函數。

在上述建立方法的優選技術方案中，所述目標函數具體為：

公式(1)中，m為地震波的速度，p0是地震觀測記錄，是p0的復共軛，是第一中間數據，pc是地震合成記錄，pc*是pc的復共軛，是第二中間數據。

在上述建立方法的優選技術方案中，「根據所述目標函數，計算地震波的速度的梯度δ(m)」的步驟進一步包括：

計算目標函數對地震波的速度的導數，應用泰勒展開，求取速度的梯度δ(m)；

公式(2)中，δp2＝(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δp2)作為復共軛震源；

s412)圖2示出波場反傳的原理示意圖。如圖2所示，對從s點正向傳播的震源，通過雙程波波動方程計算出在地下任何一個位置的正向傳播的地震波場的二階時間導數ps；

s413)繼續參照圖2，對從r點反向傳播的共軛震源pcδp2，通過雙程波波動方程，計算出在地下任何一個位置的反向傳播的地震波場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的全波形反演方案如上，但本領域技術人員能夠理解，本發明應不限於此。事實上，完全可根據實際情況，靈活地調整中間數據的具體形式，如對地震數據進行處理之後再求平方，或者平方之後再賦予相應的係數等，只要滿足由地震數據直接得出的中間數據由於與地震波強度數據相關因此必然包含低頻信息的條件即可。

至此，已經結合附圖所示的優選實施方式描述了本發明的技術方案，但是，本領域技術人員容易理解的是，本發明的保護範圍顯然不局限於這些具體實施方式。在不偏離本發明的原理的前提下，本領域技術人員可以對相關技術特徵做出等同的更改或替換，這些更改或替換之後的技術方案都將落入本發明的保護範圍之內。