[发明专利]一种基于辅助坐标系的起伏地表波形反演方法

权利要求书

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表示早至波时窗。