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

背景技術(shù)：

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

技術(shù)實(shí)現(xiàn)要素：

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

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

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

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

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

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

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

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

(3)能量加權(quán)梯度預(yù)處理

取存儲的波場到達(dá)每個速度網(wǎng)格點(diǎn)能量的最大值，即進(jìn)而求取每點(diǎn)能量的最大值為波場到達(dá)每個速度網(wǎng)格點(diǎn)能量的最大值即初至波的能量值，表征的是波傳播球面波幾何擴(kuò)散的過程，用此對梯度進(jìn)行預(yù)處理，得到能量加權(quán)梯度算子。

其中，表示能量加權(quán)梯度算子。

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

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

本發(fā)明的技術(shù)效果體現(xiàn)在：

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

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

附圖說明

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

圖2為對圖1的真實(shí)速度平滑作為初始速度

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

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

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

圖6為全波形反演能量加權(quán)梯度預(yù)處理方法流程圖

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

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

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

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

具體實(shí)施方式

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

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

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

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

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

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

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

第三步：能量加權(quán)梯度預(yù)處理

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

其中，表示能量加權(quán)梯度算子。

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

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

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