卡尔曼滤波作为一种数值估计优化方法,与应用领域的背景结合性很强.因此在应用卡尔曼滤波解决实际问题时, 最重要的是根据所应用领域知识建立相应的数学模型,再从这个模型出发进行滤波器的设计与实现工作.在反射波地震勘探领域中,地震反射波时距曲线一般近似为双曲线,其运动方程满足双曲时距方程.考虑均匀水平层状介质,假设地震波在反射界面之上介质中的传播速度为v, 炮点的自激自收时间为t₀,各检波器的炮检距分别为x₁,x₂,…,x_k,x_k+1,…,x_N,各检波道的波达时间分别为t₁,t₂,…,t_k,t_k+1,…,t_N,则由地震反射波双曲时距方程有

考虑到不可能获得地震波反射的精确模型以及许多不可预测的因素,在上述的状态方程(2)中引入了过程噪声U_k,过程噪声刻画地震波实际在地下传播过程中的运动规律与我们所建立的状态方程模型之间的误差, 例如：我们所建立的状态方程是基于地震波在单一介质中匀速传播以及反射界面为水平的前提条件下的, 而实际地下局部区域介质属性变化会引起地震波传播速度的微小变化,局部反射界面非水平会引起地震波传播方向的微小变化进而影响波达时间等,这些都可看做过程噪声(或称模型误差)来建模. 另外考虑在地震波的接收和测量过程中也会不可避免地存在噪声,在量测方程(3)中引入量测噪声W_k+1,例如：检波器的热噪声、检波器附近风吹草动引起的噪声等都可看做量测噪声(或称观测误差)来建模. 通常假设过程噪声和量测噪声是均值为零的高斯白噪声序列, 同时假定不同检波道的过程噪声、量测噪声是各自相互独立的, 且过程噪声和量测噪声序列之间互不相关.

由于状态向量两分量之间具有一定的相关性，进一步假设过程噪声U_k可用Γ_ku_k表示, 其中过程噪声分布矩阵

假设第k个检波道的状态估计为_k/k,状态估计的协方差为P_k/k,利用最小均方误差准则可逐一推导上节所建立的跟踪地震记录同相轴的卡尔曼滤波系统模型下的滤波过程. 与常规卡尔曼滤波模型相同,这里首先对状态向量、状态协方差和量测向量进行预测,它们的一步预测分别为

可见经过上述滤波过程不但拾取了同相轴位置(时间)信息,而且同时估计出了地震波在地下介质中的传播速度.

由上节的滤波模型可知, 运用卡尔曼滤波的一个重要前提条件是初始化状态估计和状态估计的协方差. 假设量测噪声W_k是均值为0、方差为σ_w²的高斯白噪声, 即,过程噪声U_k=Γ_ku_k, u_k是均值为0、方差为σ_u²的高斯白噪声, 即,且假设量测噪声与过程噪声相互独立. 此时可采用两点差分法进行初始化, 即利用最先得到的两个量测值Z₁和Z₂来初始化. 初始状态估计为

综上所述, 根据反射地震波传播时距方程, 以同相轴各子波到达时间和地震波传播速度的估计量为状态向量, 以观测到的各道地震子波到达时间为量测, 我们建立了跟踪反射同相轴的系统模型和滤波模型, 该模型融入了波达时间、炮检距、层速度等地震波传播过程中的重要物理量, 利用可观测到的各检波道地震子波的到达时间和炮检距逐次进行卡 尔曼滤波迭代, 从而得到滤波后的同相轴估计信息.

上述用于跟踪地震记录同相轴的卡尔曼滤波系统模型的输入数据是不同的炮检距和在不同炮检距道集上检测到的地震子波的波达时间. 检测波达时间是相关的另外一个研究要点, 限于篇幅在本文的仿真实验中, 仅通过简单地取各道固定窗口内地震信号幅度最大值点所对应的时间作为量测到的各道波达时间, 这在单同相轴的情况下是可行的, 多同相轴情况下的检测可通过门限方法获取.然后将各道波达时间经过平方作为卡尔曼滤波系统的量测数据, 进行滤波跟踪同相轴. 具体仿真步骤如下:

(1)根据给定仿真条件, 合成N道含噪声地震观测记录；

(2)对每道地震信号, 记录下给定窗口内幅度最大值点所对应的时间, 并将该时间平方得到量测数据Z_k(k=1,2,…,N)；

(3)估计量测噪声W_k的方差σ_w²和到达时的平方的误差的方差σ_u²，按照式(11)和(12)初始化

卡尔曼滤波器的状态估计_1/1

和状态估计的协方差P_1/1;

(4)按照2.2节中的滤波过程进行预测、滤波,逐步迭代得到滤波后状态向量序列 其中的即为经卡尔曼滤波后得到的同相轴第k道的波达时间和速度信息.

为了评价所提出方法跟踪滤波同相轴的性能, 考察如下三个评估指标(何友等，2009): 同相轴精度、发散度和收敛速度.

(1)同相轴精度

同相轴精度包括波达时间精度和速度精度. 波达时间精度定义为滤波稳定后波达时间估计的均方根误差, 速度精度定义为滤波稳定后速度估计的均方根误差, 计算公式分别为

(2)发散度η

由于量测噪声以及实际地震波传播过程与所建模型间过程噪声的存在,有时会导致滤波发散,即滤波值和真实值之间的差值随着迭代次数的增加而无限增长.有时即使该差值不无限增长,但BC明显大于检波器的量测误差, 在这种情况下，该滤波结果也是没有实用价值的.因此本文采用如下定义:若从滤波稳定时刻开始连续多道波达时间均方根误差RMSE_time大于发散检验阈值,则认为该次滤波发散,反之认为该次滤波收敛.对发散度的评估在实际应用中通常采用Monte Carlo仿真方法,即发散度定义为Monte Carlo仿真实验中滤波发散次数N₀与总仿真次数N之比,即.该值越小说明滤波器正确跟踪同相轴的概率越大.

(3)收敛速度k₀

由于滤波初始阶段数据量少以及噪声的存在, 经常导致滤波不稳定, 表现为状态向量有较强的震荡, 随着数据量增加, 经过多步迭代, 滤波将逐步稳定. 若连续多道波达时间方差小于给定的稳定检验阈值, 则认为滤波稳定. 在滤波发散的情况下考虑稳定的快慢是没有意义的, 因此定义收敛速度为在滤波收敛的条件下滤波达到稳定所需的迭代次数k₀,该次数越小则表明滤波器收敛速度越快.

本节首先给出两个利用所提出滤波方法跟踪地震反射记录同相轴的仿真示例，然后通过Monte-Carlo仿真实验较详细地分析不同参数对滤波性能的影响.

(1)简单同相轴

首先考虑简单同相轴的情形，即真实轴为具有单一速度的连续双曲轴.按照3.1节中仿真步骤, 首先依据表1所列参数按双曲时距方程产生各道到达时，合成如图1a所示无噪声地震记录, 其中地震子波采用Ricker子波(Sheriff,1994)w(t)=(1-(f_m为地震子波的主频), 前100道子波主频取为50 Hz，101-180道主频取为40 Hz，181-256道取为30 Hz. 然后叠加高斯白噪声后形成观测地震记录.图1b给出了一张含噪声记录,其 各道峰值信噪比(PSNR, Peak Signal to Noise Ratio) 为服从[0 dB, 20 dB]上均匀分布的随机数(平均 PSNR=10 dB). 为了计算简便, 这里定义PSNR为

图1

Fig. 1

图1 简单同相轴的卡尔曼滤波拾取实验示例 (a)无噪声地震记录； (b)含噪声地震记录(平均PSNR=10 dB)；(c)检测和跟踪结果； (d)图(c)的局部放大；(e)速度估计误差. Fig.1 A tracking example for a simple event by using Kalman filter (a) Noise-free seismic record；(b) Noisy record (mean PSNR=10 dB)；(c) Results of detecting and tracking； (d) Enlarged part of subplot (c)；(e) Errors of velocity estimates.

上例表明对较简单的同相轴本文方法能有效地抑制噪声，拾取出正确的同相轴.为考察本文方法对复杂同相轴情形的有效性，对图1a所示的无噪声 记录，改变地震道1~150的层速度为2000 m/s，令 第151道无有效地震子波，地震道152~256的层速度为 1400 m/s，其他合成参数如表1所示不变.图2a显示了所合成的复杂双曲同相轴的地震记录，从第152道起由于速度变慢同相轴变得更向下倾斜，而不是原来的双曲线形态.图2b叠加了高斯白噪声，图2c显示了对复杂双曲轴的检测和跟踪结果，其中实线表示真实同相轴，小圆圈表示由最大值法找到的代表反射地震子波的量测点，实心点表示通过卡尔曼滤波跟踪得到的同相轴点.由图2c可见，即使在150道附近同相轴的速度发生了突然变换，但跟踪效果依然良好，与真实同相轴基本吻合.图2d放大显示了第150道附近的跟踪情况，图2e显示了每道的速度估计误差，除了初始估计误差较大外，在速度发生变化的第150道左右速度估计误差也较大，但随着迭代的增加误差逐渐减小.这是由于当速度发生变化的时候，虽然状态向量的预测误差较大，但此时新息也将增大，从而修正了预测值，使之逐渐收敛到新的速度.

表1(Table 1)

表1 产生合成记录所用参数 Table 1 Parameters used to synthesize the seismic record 总道数 自激自收时间/s 层速度/(m·s-1) 地震子波主频/Hz 偏移距/m 道间距/m 采样间隔/s 256 1 2500 50,40,30 50 5 0.002

表1 产生合成记录所用参数 Table 1 Parameters used to synthesize the seismic record

图2

Fig. 2

图2 复杂同相轴的卡尔曼滤波拾取实验示例 (a) 无噪声地震记录； (b) 含噪声地震记录(平均PSNR=10 dB)；(c) 检测和跟踪结果；(d) 图(c)的局部放大；(e) 速度估计误差. Fig.2 A tracking example for a complex event by using Kalman filter (a) Noise-free seismic record；(b) Noisy record (mean PSNR=10 dB)；(c) Results of detecting and tracking；(d) Enlarged part of subplot (c)；(e) Errors of velocity estimates.

(3)Monte-Carlo仿真

上面两实验示例给出了卡尔曼滤波方法跟踪较简单和较复杂地震反射同相轴的结果, 由于噪声具有随机性,下面通过Monte-Carlo仿真实验计算3.2节所列各项性能评估指标, 分析不同因素对跟踪滤 波性能的影响, 其中所有Monte Carlo仿真次数均为1000.

a.地震记录信噪比

改变所叠加高斯白噪声的功率,表2列出了对所合成的具有不同PSNR地震记录滤波跟踪后的同相轴各项性能评估指标值.由表2可见,随着PSNR的提高,滤波后同相轴精度越来越高,发散度越来越小,收敛速度也缓慢加快.

表2(Table 2)

表2 不同PSNR下的性能评估 Table 2 Performances of the proposed method for the cases of different PSNRs PSNR/dB 8 9 10 11 12 13 14 15 RMSEtime/s 6.7×10-4 6.5×10-4 6.3×10-4 6.2×10-4 6.1×10-4 6.0×10-4 6.0×10-4 5.9×10-4 RMSEv/(m·s-1) 31 27 25 23 21 19 19 17 发散度η 0.047 0.024 0.008 0.003 0.002 0 0.001 0.001 收敛速度k0 37 34 33 33 31 30 30 30

表2 不同PSNR下的性能评估 Table 2 Performances of the proposed method for the cases of different PSNRs

b.层速度

一般说来，层速度越高同相轴的形态表现为越陡峭, 为了考察本文方法对不同形态同相轴的跟踪能力, 改变表1中层速度值, 其它条件不变, 并固定PSNR=10 dB, 得到表3所列不同层速度下滤波跟踪性能评估指标值. 可见层速度越大, 速度估计精度越差, 收敛速度越慢, 到达时精度和发散度受层速度影响较小. 造成速度估计精度随层速度增大而变差的原因部分上是由于采样间隔(0.02 s)较大的缘故.

表3(Table 3)

表3 不同层速度下的性能评估 Table 3 Performances of the proposed method for the cases of different velocities 层速度/m·s-1 500 1000 1500 2000 2500 3000 3500 RMSEtime/s 0.0275 6.4×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.4×10-4 RMSEv/(m·s-1) 3 4 6 13 24 43 69 发散度η 0.005 0.007 0.011 0.011 0.013 0.009 0.009 收敛速度k0 5 12 19 26 32 41 48

表3 不同层速度下的性能评估 Table 3 Performances of the proposed method for the cases of different velocities

c.量测噪声标准差参数

在卡尔曼滤波系统中, 一般要事先估计量测噪声标准差, 然后作为一个输入参数传递给滤波过程. 在PSNR=10 dB的条件下, 表4列出了不同量测噪声标准差参数对滤波器性能的影响.可见当该参数较小(如:0.001,0.002)或较大(如:0.05,0.1,1)时会导致各项指标变差, 如:发散度变得较大, 即正确跟踪真实同相轴的概率较低, 因此设定合适的量测噪声标准差参数是十分重要的.

表4(Table 4)

表4 不同量测噪声标准差参数下的性能评估 Table 4 Performances of the proposed method for the cases of different σw σw/w 0.001 0.002 0.003 0.004 0.005 0.01 0.05 0.1 1 RMSEtime/s 0.0286 0.0012 6.3×10-4 6.3×10-4 6.3×10-4 6.6×10-4 0.0027 0.0076 0.0082 RMSEv/(m·s-1) 127 28 24 24 24 31 201 331 346 发散度η 0.880 0.225 0.025 0.009 0.010 0.043 0.257 0.482 0.503 收敛速度k0 28 31 33 33 35 35 64 89 93

表4 不同量测噪声标准差参数下的性能评估 Table 4 Performances of the proposed method for the cases of different σw

d.门限参数γ

门限参数γ影响着当前量测点能否参与本次迭代的滤波更新, γ越大则量测点参与滤波更新的概率就越大.设置门限γ的作用是剔除一些误差较大的量测点对跟踪同相轴的影响.表5列出了相同条件下, 不同γ参数对滤波器性能的影响.可见该参数设置过小(如:3)或过大会导致较大的发散度.前者是由于门限较小会导致部分量测点不能参与滤波更新过程，使滤波仅沿着原始预测方向迭代的缘故；后者是由于门限较大导致过多的误差较大的量测点参与滤波更新的缘故.

表5(Table 5)

表5 不同门限参数下的性能评估 Table 5 Performances of the proposed method with different parameters γ γ 3 5 7 9 11 13 15 17 20 RMSEtime/s 8.2×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 6.3×10-4 RMSEv/(m·s-1) 25 24 23 24 24 24 24 23 24 发散度η 0.145 0.036 0.020 0.012 0.012 0.015 0.017 0.018 0.018 收敛速度k0 32 33 33 34 32 34 32 33 35

表5 不同门限参数下的性能评估 Table 5 Performances of the proposed method with different parameters γ

由上可见, 本文方法跟踪滤波性能除受层速度影响较小外,地震记录信噪比以及预设的量测噪声标准差参数和门限参数都会不同程度地影响跟踪性能，适宜参数的选取还需进一步研究.