[發明專利]基于時頻域目標函數的鹽下構造最小二乘逆時偏移方法

有效

• 申請號: 202010430122.X
• 申請日: 2020-05-20
• 公開（公告）號: CN111999764B
• 公開（公告）日: 2021-04-13
• 發明（設計）人: 胡勇;許永忠;胡明順;楊磊
• 申請（專利權）人: 中國礦業大學
• 主分類號: G01V1/28
• 分類號: G01V1/28;G01V1/32;G01V1/34
• 代理公司: 江蘇圣典律師事務所 32237
• 代理人: 程化銘

1.一種基于時頻域目標函數的鹽下構造最小二乘逆時偏移方法，其步驟是：

步驟1.使用MATLAB軟件及其地震數據處理軟件Crews工具包對地震數據進行預處理，并將預處理后的地震數據作為時頻域最小二乘逆時偏移的數據輸入；

步驟2.構建初始速度模型，作為波動方程正演模擬的模型輸入；

步驟3.定義觀測系統，提取震源子波；利用初始速度模型進行波動方程正演模擬，存儲正傳波場，并與觀測數據的反傳波場互相關，獲得逆時偏移成像結果；在速度模型上正演模擬對應的聲波方程為：

其中x,z為坐標軸，s表示慢度場，u為地震波場，f為震源子波，t為時間；

步驟4.計算擾動介質產生的擾動地震波場，擾動介質表示為：

將慢度場的平方s²表示為背景慢度的平方與擾動慢度平方Δs²的疊加；對應的總地震波場同樣分解為背景地震波場u₀和擾動地震波場u_s的疊加，即：

u＝u₀+u_s

其中背景地震波場u₀與總地震波場u都滿足波動方程：

當擾動地震波場u_s＜＜u₀時，用背景地震波場代替總地震波場u₀≈u₀+u_s，即：

步驟5.利用Gabor變換獲得觀測數據和模擬數據的時頻域振幅相位信息；其中觀測數據和模擬數據Gabor變換為：

其中h(τ-t)表示高斯窗函數；t表示時間；τ為中間變量；ω為角頻率，u(t)和d(t)分別表示時間域模擬數據和觀測數據；和表示時頻域模擬數據和觀測數據，F_h[·]表示對地震數據作用的Gabor變換算子；

步驟6.構建基于時頻域振幅相位的最小二乘逆時偏移目標函數：

其中i表示虛數；ns為震源數目；nr為檢波器數目；模型擾動和表示模擬數據和觀測數據在時頻域的相位信息；和表示觀測數據和模擬數據在時頻域振幅信息；ε∈[0,1]表示權重因子，用來控制目標函數中振幅和相位的比例；當ε＝0時，目標函數為純相位目標函數，完全忽略了波形振幅的影響；當ε＝1時，目標函數將退化為常規最小二乘目標函數；求目標函數關于模型參數的偏導數，獲得目標函數對模型擾動的更新量：

步驟7.利用Gabor反變換將時頻域的梯度表達式轉換到時間域中，其中為了對公式做進一步的簡化，地震觀測數據和模擬數據Gabor逆變換表示為：

其中表示Gabor逆變換算子；因此，目標函數的梯度進一步在時空域簡化為：

步驟8.定義伴隨震源，并將其反傳至模型空間；對應的伴隨震源為：

步驟9.利用零延遲互相關方法獲得最小二乘逆時偏移的梯度；根據波動方程的表達式，

Born正演模擬用矩陣的形式來表示：

A_su_s＝A₀u₀

其中計算其關于模型擾動的偏導數：

同樣將梯度寫成矩陣的形式為：

步驟10.利用L-BFGS優化算法計算模型擾動的更新方向；迭代公式如下：

m_k+1＝m_k-α_kH_kg_k

其中m_k為第k步模型擾動的更新量，α_k為步長，H_k為近似海森矩陣的逆，為模型擾動的梯度；在L-BFGS優化更新中，只需要保存少數的向量對，用于更新Hessian矩陣，其更新公式如下：

H_k+1＝V_k^TH_kV_k+ρ_ks_ks_k^T

V_k＝I-ρ_ky_ks_k^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中H_k+1是根據向量對{s_k,y_k}和H_k計算得到；H_kg_k的乘積通過梯度g_k與向量對{s_k,y_k}之間一系列向量的內積與向量的和來獲得；其中的近似Hessian矩陣的逆矩陣H_k需滿足以下更新公式：

步驟11.判斷是否滿足終止條件，若滿足則輸出基于時頻域目標函數的最小二乘逆時偏移結果；若不滿足終止條件，將當前的成像結果繼續作為下一次循環的初始擾動模型，直到滿足終止條件。

2.根據權利要求1所述的基于時頻域目標函數的鹽下構造最小二乘逆時偏移方法，其特征是：采用有限差分法進行波動方程正演模擬。