本發(fā)明涉及地球物理勘查技術領域，主要用于提高三維全波形反演縱波速度場的精度。

背景技術：

速度是描述地下介質情況的重要參數，地球物理勘探的重點在于如何恢復地下介質所有尺度的速度信息。速度分析的方法有很多，但是始終達不到理想的效果。射線層析、波動方程層析等層析成像的方法可恢復低波數的速度信息，偏移的方法可提供高波數的反射率信息，然而這些方法都不能同時恢復所有波數的速度參數。全波形反演利用疊前地震數據的全部信息，可提供地下介質的高分辨率成像。常規(guī)全波形反演是非線性梯度類最優(yōu)化方法，以觀測數據與模擬數據的最小二乘目標函數值最小為標準更新地下速度模型。三維全波形反演主要存在以下難點：(1)巨大的計算量問題：三維全波形反演基于三維波動方程正演模擬的迭代反演方法，在一次迭代過程中，需要計算多炮的波動方程正演模擬，誤差波場反向傳播，這兩部分每個時間步都需要進行差分運算，伴隨著巨大的計算量。(2)海量存儲的問題：基于伴隨狀態(tài)梯度算子的全波形反演方法要求存儲整個正演波場，三維情況下的存儲量是不可容忍的，可以說大部分的時間消耗在數據的I/O讀寫上。計算量大可以通過多種并行算法進行加速，而I/O讀寫速度只能靠提高硬件性能來改善。(3)梯度算子的求取效率與精度：常規(guī)全波形反演梯度算子是由震源正傳波場的時間二階偏導數與殘差波場反傳波場互相關所得，震源正傳波場包含著球面波傳播過程中幾何擴散的能量損失，震源波場與殘差反傳波場互相關就導致了梯度深、淺層能量更加不均衡，會引起反演速度深層精度不足。(4)海森矩陣的求?。簩μ荻阮A處理同樣是全波形反演中的一個重要環(huán)節(jié)，可有效的提高反演的精度。然而海森矩陣的求取伴隨著巨大的存儲量和計算量，如何避開海森矩陣的求取或是對海森矩陣進一步的近似都是需要考慮的問題。如何在提高效率、降低耗費且保證精度是三維全波形反演的首要問題。因此需要發(fā)展一種三維全波形反演能量加權梯度預處理方法。

技術實現要素：

本發(fā)明的目的正是為了解決現有技術存在的問題，提供一種高效高精度的三維全波形反演能量加權梯度預處理方法。

本發(fā)明的技術方案為：

一種三維全波形反演能量加權梯度預處理方法，具體步驟如下：

(1)三維高精度有限差分正演模擬

第一，根據初始速度的道頭信息確定三維正演模擬的觀測系統(tǒng)，主要確定炮檢波點的位置；第二，根據初始速度的最大最小值求取滿足有限差分數值模擬差分穩(wěn)定性和頻散關系的三維正演模擬參數；第三，引入完全匹配層邊界條件，確定其需要的參數，用于消除正演模擬的邊界反射干擾；第四，用時間二階空間十階有限差分方法進行三維正演模擬，存儲邊界波場用于波場重建，并計算存儲到達每個速度網格點振幅最大值。

(2)源波場重建及誤差波場逆時傳播求取梯度

正演模擬得到的模擬炮集與實際觀測炮集對應做差求取殘差波場。讀取步驟(1)中存儲的邊界波場作為邊界條件用時間二階空間十階有限差分方法進行源波場重建。源波場重建的同時進行殘差波場逆時傳播。根據式(1)的伴隨狀態(tài)法進行梯度計算，將源波場重建的波場與殘差波場逆時傳播得到的波場對應時刻進行零延遲互相關得到梯度算子。

其中E為目標函數，表示梯度，m為模型參數，v表示各網格點的速度值，x_s表示震源點的位置，x表示各網格點的位置，t表示每個時間步，T表示最大時間步，為源波場關于時間的二階偏導數，q為以殘差波場為震源的逆時反傳波場。

(3)能量加權梯度預處理

取存儲的波場到達每個速度網格點能量的最大值，即進而求取每點能量的最大值為波場到達每個速度網格點能量的最大值即初至波的能量值，表征的是波傳播球面波幾何擴散的過程，用此對梯度進行預處理，得到能量加權梯度算子。

其中，表示能量加權梯度算子。

(4)求取合適步長迭代更新速度

首先給一個試探步長，再用Armijo條件的一維線搜索方法求取合適的步長作用與能量加權梯度算子對速度進行迭代更新。

本發(fā)明的技術效果體現在：

常規(guī)全波形反演梯度算子是由震源正傳波場的時間二階偏導數與殘差波場反傳波場互相關所得，震源正傳波場包含著球面波傳播過程中幾何擴散的能量損失，震源波場與殘差反傳波場互相關就導致了梯度深、淺層能量更加不均衡，會引起反演速度深層精度不足。利用海森矩陣對梯度進行預處理可有效的均衡梯度的能量，但是三維情況下的海森矩陣構建困難，計算量與存儲量巨大，因此我們只提取海森矩陣中球面波傳播幾何擴散的信息對梯度進行預處理達到均衡梯度能量的效果，提高速度深層的反演精度。

該方法基于GPU集群搭建多級異構并行算法解決三維全波形反演的運算效率問題；利用源波場重建策略解決三維情況下的海量存儲問題；在保證運算效率的情況下提高梯度算子構建精度。

附圖說明

圖1為截取的部分二維Marmousi模型并重采樣作為真實速度

圖2為對圖1的真實速度平滑作為初始速度

圖3為不同全波形反演算法迭代60次速度的對比，其中圖a為常規(guī)算法反演速度；圖b為本文算法反演速度

圖4為常規(guī)全波形算法迭代60次與本文算法迭代60次誤差下降曲線的對比

圖5為不同全波形反演算法迭代60次單道速度的對比，其中圖a為常規(guī)算法反演速度；圖b為本文算法反演速度

圖6為全波形反演能量加權梯度預處理方法流程圖

圖7為三維SEG/EAGE推覆體速度模型，其中圖a為真實速度；圖b為初始速度

圖8為三維SEG/EAGE推覆體模型不同迭代次數全波形反演速度，其中圖a為第1次迭代的三維速度體與三個不同方向的剖面；圖b為第10次迭代的三維速度體與三個不同方向的剖面；圖c為第40次迭代的三維速度體與三個不同方向的剖面

圖9為三維SEG/EAGE推覆體模型不同迭代次數的誤差下降曲線

圖10為三維SEG/EAGE推覆體模型最終反演速度與真實速度以及初始速度的對比，其中圖a為真實速度；圖b為初始速度；圖c為反演速度

具體實施方式

下面結合附圖說明本發(fā)明的具體實施方式：

通過模型測試說明具體的技術方案：

第一步：三維高精度有限差分正演模擬

第一，根據初始速度的道頭信息確定三維正演模擬的觀測系統(tǒng)，主要確定炮檢波點的位置；第二，根據初始速度的最大最小值求取滿足有限差分數值模擬差分穩(wěn)定性和頻散關系的三維正演模擬參數；第三，引入完全匹配層邊界條件，確定邊界厚度以及衰減吸收系數，波場能量到達邊界處逐漸衰減為零，從而消除三維正演模擬的邊界反射干擾；第四，用時間二階空間十階有限差分方法對速度模型進行三維正演模擬，存儲邊界波場用于波場重建，并計算波場到達每個速度網格點振幅最大值。

第二步：源波場重建及誤差波場逆時傳播求取梯度

第一步中三維正演模擬得到的模擬炮集與實際觀測炮集對應做差求取殘差波場。讀取第一步中存儲的邊界波場作為邊界條件，用時間二階空間十階有限差分方法進行源波場重建。源波場重建的同時進行殘差波場逆時傳播。根據下式(1)的伴隨狀態(tài)法進行梯度計算，將源波場重建的波場與殘差波場逆時傳播得到的波場對應時刻進行零延遲互相關得到梯度算子。

式中E為目標函數，表示梯度，m為模型參數，v表示各網格點的速度值，x_s表示震源點的位置，x表示各網格點的位置，t表示每個時間步，T表示最大時間步，為源波場關于時間的二階偏導數，q為以殘差波場為震源的逆時反傳波場。

第三步：能量加權梯度預處理

讀取存儲的源波場到達每個速度網格點振幅的最大值，即進而求取每點能量的最大值為(前一個公式是最大振幅值，后一個公式是最大能量值)。源波場到達每個速度網格點能量的最大值即初至波的能量值，表征的是波傳播球面波幾何擴散的過程，用此對梯度進行預處理，得到能量加權梯度算子，

其中，表示能量加權梯度算子。

第四步：求取合適步長迭代更新速度

每次迭代，首先給一個試探步長，再用一維線搜索方法求取新的步長，此步長滿足Armijo條件使誤差下降，然后將求得的步長作用于能量加權梯度算子對速度進行迭代更新。圖3到圖5展示的是二維Marmousi模型反演測試的效果圖。從圖3可以看出圖b中的速度成像淺層與深層的能量較圖a更加均衡。而從圖4誤差曲線下降圖的也可以看出文中基于能量加權梯度算子的全波形反演算法誤差收斂更快，且收斂的誤差值更小。從圖5為單道速度反演結果的對比可以看出基于能量加權梯度算子的全波形反演算法深部速度更新更加準確。

圖8到圖10展示的是按照圖6流程進行三維SEG/EAGE推覆體模型反演測試的效果圖。圖8展示的是不同迭代次數三維SEG/EAGE推覆體模型全波形反演的速度結果。其中圖a為第1次迭代的三維速度體與三個不同方向的剖面；圖b為第10次迭代的三維速度體與三個不同方向的剖面；圖c為第40次迭代的三維速度體與三個不同方向的剖面。三個不同方向分別為橫縱測線方向8000m處的剖面與深度3000m處的切片。從圖7中可以看出隨著迭代次數的增加，構造細節(jié)成像越來越清晰，特別對于3000m處目的層切片的河道處的細節(jié)恢復明顯。圖9是迭代過程中的誤差下降曲線，誤差做過歸一化。從圖中可以看出，17次之前誤差下降迅速，17次之后逐漸趨向收斂。圖10是三維SEG/EAGE推覆體模型全波形反演迭代第40次的最終反演速度。圖a是初始速度，圖b是真實速度，圖c是全波形反演速度。從圖中可以看出速度深淺層成像效果均衡，并可以看出全波形反演在構造細節(jié)上優(yōu)勢明顯，可以對類似于河道等的構造細節(jié)精細刻畫。