[发明专利]一种高分辨率的多尺度波动方程反演方法

权利要求书

1.一种高分辨率的多尺度波动方程反演方法，其特征在于包括：

步骤101，采集炮集数据；

步骤102，输入地层层位数据H(x)；

步骤103，建立约束模型m_p(x,z)＝(K_p,ρ_p),其中，输入的模型是所需反演的储层段的体积模量K_p和密度ρ_p；

步骤104，记每个炮集数据为一个二维数据P_o(t,x)，其中，x和z代表空间坐标，t代表时间变量，设炮集个数为n_s，对炮集数据P_o(t,x)进行多尺度分解，得到在不同主频尺度下的地震炮集P_o(t,x,ω_i)，ω_i＝i·Δω,多尺度分解公式可表示为：

G(ω)=12πσe-ω22σ2]]>

其中，i表示虚数，ω表示高斯函数的中心频率，Δω表示频率采样率；σ表示高斯函数的宽度，多尺度下地震炮集可表示为

P_o(t,x,ω_i)＝FFT^-1[FFT(P(t,x))·G(ω_i)]

对每个尺度的炮集数据P_o(t,x,ω_i)，重复下述步骤105至步骤111；

步骤105，设初始模型的输入m₀(x,z)＝(K₀,ρ₀)，K₀为所需反演储层段的初始体积模量，ρ₀为所需反演储层段的初始密度；

步骤106，针对频率ω_i，以主频为ω_i的零相位雷克子波，根据声波波动方程及吸收边界条件计算模拟炮集数据P_s(t,x,ω_i)；

步骤107，计算接收到的地震炮集数据P_o(t,x,ω_i)与模拟炮集数据P_s(t,x,ω_i)的残差d(t,x)，记为

d(t,x)＝P_o-P_s(4)

其中，P_o为P_o(t,x,ω_i)，P_s为P_s(t,x,ω_i)；

步骤108，以-d(-t,x)为震源，替换声波波动方程的震源项，模拟残差记录的反向炮集数据P_s'；

步骤109，根据上述计算得到的炮集数据和反向炮集数据构造反演模型变化的梯度：

g=∂∂mE(p)=<δρ,∂∂tPs·∂∂tPs′>+<δK,▿Ps·▿Ps′>---(5)]]>

E(p)为炮集数据能量误差，上述式(5)的反演模型参数的共轭修改量可表示为：

δρiω,js=1ρ2(x)∫0T∂∂tPs(r,t;xs;ωi)∂∂tPs′(r,t;xs;ωi)dt---(6)]]>

δKiω,js=1K2(x)∫0T▿Ps(r,t;xs;ωi)·▿Ps′(r,t;xs;ωi)dt---(7)]]>

其中，r为接收点的位置，t为时间变量，x_s为震源的位置，表示梯度算子，T为炮集记录时间长度，K,ρ为反演模型参数，为反演模型参数修改量，i_ω表示第i个频率，j_s表示第j炮；对该频率下的炮集中所有炮重复步骤106至步骤109，累加计算整个模型修改量；

步骤110，利用步骤109计算的所有模型修改量，生成新的反演模型，

所述新的反演模型为：

m_n+1＝m_n-α_ng_n(8)

其中，m＝(K,ρ)为反演模型，n为迭代次数，α_n为第n次迭代步长，g_n＝(δK,δρ)为第n次迭代的反演模型变化梯度；

步骤111，按照下式计算能量误差，

E=12ΔdtCD-1Δd---(9)]]>

式中E为能量误差，Δd^t为Δd的转置，Δd为上述步骤107中的残差，C_D为数据协方差矩阵；

当E＜δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为本频段内下次反演的初始模型，进入重复步骤105；其中i_ω为第i个频率，n_ω为频率采样个数；

当E＞δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为下一频段反演的初始模型，进入重复步骤104；

当E＜δ，i_ω＝n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为最终的反演结果，其中，δ为任意小的实数。