[發明專利]一種基于聲波方程的線性優化隱式時空域有限差分數值模擬方法

有效

• 申請號: 201711142912.2
• 申請日: 2017-11-17
• 公開（公告）號: CN107976710B
• 公開（公告）日: 2019-05-28
• 發明（設計）人: 王恩江;巴晶;吳春芳;戚慧;郭夢秋;周健
• 申請（專利權）人: 河海大學
• 主分類號: G01V1/28
• 分類號: G01V1/28;G01V1/36
• 代理公司: 南京縱橫知識產權代理有限公司 32224
• 代理人: 董建林;姚蘭蘭

1.一種基于聲波方程的線性優化隱式時空域有限差分數值模擬方法，其特征在于，包括以下幾個步驟：

(1)讀取參數，所述參數包括正演模擬所需要的速度模型參數文件、有限差分算子長度、子波函數及子波主頻、正演所采用的時間與空間步長及地震記錄時長；

(2)基于菱形差分算子和二階時間中心差分，得到時間導數的時間高階離散格式；

(3)求解二階空間導數隱式差分格式，在顯式高階差分的基礎上，分母中引入二階中心差分格式，構造空間隱式差分格式求解二階空間導數，用于減小有限差分算子長度；進一步，基于最小二乘法，以L2范數目標函數極小化空間域頻散關系，求解優化的隱式差分系數；

(4)基于所述時間高階離散格式和空間隱式離散格式，帶入聲波方程，構建同時具有時間高階和空間隱式的差分遞推格式，得到時空域頻散關系；

(5)固定空間差分系數，采用優化策略極小化整個遞推格式頻散關系，建立關于時間差分格式中差分系數的L2范數目標函數，利用線性優化算法對時間差分系數進行求解；

(6)采用混合吸收邊界條件對邊界反射能量進行吸收，按照時間高階和空間隱式差分格式離散波動方程并進行遞推，得到任意時刻的波場及整個地震記錄；

(7)記錄波場快照，輸出地震記錄并結束；

步驟(2)中，求解二階時間導數時采用中心差分，形式如下：

其中，為壓力場，τ為時間步長，t為時間；為了進一步提高時間離散精度，將菱形差分算子引入到公式(1)中，得出如下的時間高階離散公式：

其中，h為空間采樣步長，v為速度，a_m,n為有限差分系數；N為菱形算子中的差分算子長度；n和m為菱形算子中的坐標位置序列；

步驟(3)中，采用如下空間隱式離散格式對二階空間導數進行近似，形式如下：

其中，a′₀、a′_m和b為隱式差分系數，當b＝0時，公式(3)退化為顯式差分格式；h為空間步長，x為方向坐標，M為差分算子長度；公式(3)對應的空間頻散關系如下：

其中，β＝k_xh，k_x為水平方向波數；極小化公式(4)對應的空間頻散關系，目標函數為關于差分系數的線性凸函數，最優差分系數通過求解如下公式得到：

其中，n＝1,2,L,M+1，d＝[d₁,d₂,L,d_M,d_M+1]^T＝[a′₁,a′₂,L,a′_M-1,a′_M,b]^T，q為介于0和1之間的自然數，控制著積分范圍，

步驟(4)中，在求解時間導數時，采用公式(2)計算；求解空間導數時，采用公式(3)計算，將公式(2)和公式(3)帶入到聲波方程并進行平面波理論分析，過程如下：

均勻模型中，聲波方程壓力P滿足如下平面波方程，

其中，ω是角頻率，是虛數單位，k_z為垂直方向波數；m、n和l為位置和時間序列；

將公式(7)帶入公式(2)并進行化簡得到如下公式：

同樣，對空間導數公式進行分析，得到如下頻散關系：

將上述(8)、(9)和(10)帶入聲波方程并簡化和重組得到如下時空域頻散關系，

其中，ω為角頻率，k_x和k_z分別為沿x和z方向的波數，r為庫朗數，r＝vτ/h；

步驟(5)中，首先，差分系數a′_m和b由公式(5)計算得到，在計算差分系數a_m,0和a_m,n時保持固定，此時優化差分系數a_m,0和a_m,n是線性優化問題，直接線性求解，定義如下目標函數：

其中，d＝[d₁,L,d_N+(N-1)N/2]^T＝[a_1,0,L,a_N,0,a_1,1,L,a_1,N-1,L,a_N-1,1]^T，表示差分系數；θ為傳播角度，q為介于0和1間的實數，波數k_x和k_z與波數k之間滿足如下關系：

k_x＝kcos(θ),k_z＝ksin(θ), (13)ψ(β,θ,r)及定義如下：

目標函數是關于差分系數a_m,0和a_m,n的線性函數，差分系數a_m,0和a_m,n通過求解如下公式得到，

2.根據權利要求1所述的基于聲波方程的線性優化隱式時空域有限差分數值模擬方法，其特征在于，步驟(6)中，整個計算區域被劃分為內部區域、邊界區域和過渡區域，具體的方法如下：

(6-1)利用公式(2)和公式(3)離散時間和空間導數，采用公式(5)和公式(16)計算差分系數，進行雙程波方程求解，求解聲波方程，計算所有區域處的雙程波波場；

(6-2)求解單程波方程，計算邊界區域和過渡區域處的單程波波場，計算公式及方法如下：

針對二階聲波方程，以左邊界和左上角為例，左邊界單程波方程為：

其中，x和z代表坐標；

左上角單程波方程為：

其他邊界處及角點處單程波可以類似得到，并進行對應求解；采用上述吸收邊界條件策略計算各位置處任意時刻波場值；

(6-3)計算最終波場：

邊界處波場為單程波場值，內部區域處波場為雙程波場值，中間過渡區域處波場通過如下線性組合得到：

P＝B_iP^two+(1-B_i)P^one, (19)

其中，P^two為雙程波場，P^one為單程波場，B_i為權重系數，B_i＝(i-1)/N_b，i＝2,L,N_b，N_b為過渡區域吸收邊界層數，過渡區域從內到外，B_i由(N_b-1)/N_b變化到1/N_b；如果N_b＝1，則此時混合吸收邊界條件退化為單程波吸收邊界條件。