[發明專利]一種基于輔助坐標系的起伏地表波形反演方法

有效

• 申請號: 201610037642.8
• 申請日: 2016-01-20
• 公開（公告）號: CN105549080B
• 公開（公告）日: 2017-08-25
• 發明（設計）人: 曲英銘;李振春;李金麗;黃建平
• 申請（專利權）人: 中國石油大學（華東）
• 主分類號: G01V1/28
• 分類號: G01V1/28
• 代理公司: 濟南舜源專利事務所有限公司37205
• 代理人: 王連君

1.一種基于輔助坐標系的起伏地表波形反演方法，其特征在于：按照如下步驟進行：

步驟1：輸入初始全局速度場、常規疊前炮記錄、起伏高程及震源子波，并建立觀測系統；

步驟2：根據初始全局速度場及起伏高程進行網格剖分，將近地表附近的網格剖分成曲網格，深層的網格剖分成矩形網格；

步驟3：將初始全局速度場變換到輔助坐標系下的矩形網格，采用下式所示的變換格式：

x(ξ,η)=ξz(ξ,η)=zi-1(ξ)-zi(ξ)ηi-1(ξ)-ηi(ξ)(η-ηi(ξ))+zi(ξ);]]>

其中，x和z表示笛卡爾坐標系下的橫縱坐標；ξ和η表示輔助坐標系下的橫縱坐標；z_i-1(ξ)和z_i(ξ)是笛卡爾坐標系下第i層頂界面、底界面的高程，定義最深層高程為零；η_i-1(ξ)和η_i(ξ)是輔助坐標系第i層頂界面、底界面的縱向采樣點數，定義最深層縱向采樣點數為零；

步驟4：對常規疊前炮記錄劃分時窗，在輔助坐標系下對近地表附近的曲網格區域速度場應用早至波波形反演，更新近地表的速度場，早至波波形反演的梯度方向如下：

g(v)=2v3ΣxsΣt=0Te∂p∂t·p*;]]>

其中，g為梯度；x_s為炮點坐標；v為介質速度；p為聲壓；p^*是殘差波場的反向傳播；t為時間；T_e表示早至波時窗；

步驟5：判斷使用近地表速度場正演模擬的炮記錄與步驟4中劃分時窗后的常規疊前炮記錄之差是否滿足誤差條件；

若：判斷結果是使用近地表速度場正演模擬的炮記錄與步驟4中劃分時窗后的常規疊前炮記錄之差滿足誤差條件，則近地表速度場更新完成，然后執行步驟6；

或判斷結果是使用近地表速度場正演模擬的炮記錄與步驟4中劃分時窗后的常規疊前炮記錄之差不滿足誤差條件，則執行步驟4；

步驟6：將步驟1中輸入的常規疊前炮記錄分解成不同主頻的多尺度炮記錄，采用如下式所示的分解公式：

Fw(ω)=Wt(ω)Wo*(ω)|Wo(ω)|2+ϵ2]]>

其中，F_w是維納濾波器；W_o表示初始炮記錄子波；W_t是生成的炮記錄子波；ω是角頻率；ε為一個小數；*表示共軛轉置；

步驟7：將更新過近地表速度的全局速度場作為初始速度場，在輔助坐標系下應用全波形反演從低頻到高頻更新全局速度場，全波形反演的梯度方向如下：

g(v)=2v3ΣxsΣt=0TmaxΣf=f1fmax∂p∂t·pf*]]>

其中，f表示主頻；f₁和f_max為多尺度分解的最低頻率和最高頻率；T_max表示常規疊前炮記錄的最大記錄時間；表示主頻是f的殘差波場反傳；

步驟8：判斷使用全局速度場正演模擬的炮記錄與常規疊前炮記錄之差是否滿足誤差條件；

若：判斷結果是使用全局速度場正演模擬的炮記錄與常規疊前炮記錄之差滿足誤差條件，則全局速度場更新完成，然后執行步驟9；

或判斷結果是使用全局速度場正演模擬的炮記錄與常規疊前炮記錄之差不滿足誤差條件，則執行步驟7；

步驟9：將更新完成的全局速度場反變換到笛卡爾坐標系下，反變換公式如下：

ξ(x,z)=xη(x,z)=ηi-1(ξ)-ηi(ξ)zi-1(ξ)-zi(ξ)(z-zi(ξ))+ηi(ξ);]]>

步驟10：輸出反演的速度場。

2.根據權利要求1所述的基于輔助坐標系的起伏地表波形反演方法，其特征在于：在步驟4中，具體包括

步驟4.1：定義目標函數：

E=12ΣxsΣxrΣt(Ru(t,xr,xs)-dobs(t,xr,xs))2---(1);]]>

其中，u(t,x_r,x_s)代表模擬波場u＝(u,w,p)^T，其中T表示轉置；R為限定算子；d_obs(t,x_r,x_s)是常規疊前炮記錄；x_s和x_r表示震源點和檢波點的位置坐標；t表示時間；E為目標函數值；

步驟4.2：將目標函數變分，得到變分表達式：

δE(v)=12δ<(Ru-dobs),(Ru-dobs)>=<δ(Ru-dobs),(Ru-dobs)>=<(Rδu),(Ru-dobs)>---(2);]]>

步驟4.3：定義變換格式：

x(ξ,η)=ξz(ξ,η)=zi-1(ξ)-zi(ξ)ηi-1(ξ)-ηi(ξ)(η-ηi)+zi(ξ)---(3);]]>

通過變換格式(3)與鎖鏈法則得到下面的映射公式：

∂ξ∂x=1∂ξ∂z=0∂η∂z=ηi-1-ηizi-1(ξ)-zi(ξ)∂η∂x=ηi-1-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηizi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ---(4);]]>

通過映射公式(4)，得到輔助坐標系下的一階方程：

ρ∂u∂t-∂p∂ξ-∂p∂η(ηi-1(ξ)-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηi(ξ)zi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)=0ρ∂w∂t-∂p∂η(ηi-1(ξ)-ηi(ξ)zi-1(ξ)-zi(ξ))=01v2∂p∂t-ρ[∂u∂ξ+∂u∂η(ηi-1(ξ)-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηi(ξ)zi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)+∂w∂η(ηi-1(ξ)-ηi(ξ)zi-1(ξ)-zi(ξ))]=s---(5);]]>

其中，p是聲壓；u和w分別是水平方向和垂直方向的質點速度；v是介質速度；s表示震源項；ρ為密度；

將v+δv→u+δu代入一階方程(5)后與一階方程(5)相減得到如下方程：

ρ∂δu∂t-∂δp∂ξ-∂δp∂η(ηi-1-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηizi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)=0ρ∂δw∂t-∂δp∂η(ηi-1-ηizi-1(ξ)-zi(ξ))=01v2∂δp∂t-ρ[∂δu∂ξ+∂δu∂η(ηi-1-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηizi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)+∂δw∂η(ηi-1-ηizi-1(ξ)-zi(ξ))]=2δvv3∂p∂t---(6);]]>

進一步可得：

δu=L(2δvv3∂p∂t)---(7);]]>

其中，L代表正演算子；

步驟4.4：將方程(7)代入變分表達式(2)，可得：

δE(v)=<R(L(2δvv3∂p∂t)),(Ru-dobs)>=<(2δvv3∂p∂t),L*R*(Ru-dobs)>---(8);]]>

其中，L^*R^*(Ru-d_obs)表示波場的逆時傳播；

利用伴隨狀態法，L^*R^*(Ru-d_obs)可由下式求得：

ρ∂u*∂t-∂p*∂ξ-∂p*∂η(ηi-1-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηizi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)=0ρ∂w*∂t-∂p*∂η(ηi-1-ηizi-1(ξ)-zi(ξ))=01v2∂p*∂t-ρ[∂u*∂ξ+∂u*∂η(ηi-1-ηzi-1(ξ)-zi(ξ)·∂zi(ξ)∂ξ+η-ηizi-1(ξ)-zi(ξ)·∂zi-1(ξ)∂ξ)+∂w*∂η(ηi-1-ηizi-1(ξ)-zi(ξ))]=p-pobs---(9);]]>

則L^*R^*(Ru-d_obs)＝p^*；

步驟4.5：求取目標函數對速度模型的梯度，得到早至波波形反演的梯度方向：

g(v)=δE(v)δv=2v3ΣxsΣt=0Te∂p∂t·p*---(10);]]>

其中，g為梯度；x_s為炮點坐標；v為介質速度；p為聲壓；p^*是殘差波場的反向傳播；t為時間；T_e表示早至波時窗。