地震波是照亮地球内部的一盏明灯,随着地震观测和数字信号处理技术的发展,地震波不仅可以提供地球内部的结构信息,其波速变化还能反映地球内部的介质变化信息(陈颙,朱日祥,2005; 陈颙等,2007a,b; 丘学林等,2007; Chen et al,2008,2010)。而其中,基于背景噪声格林函数的波速变化和基于精确定位的相似地震研究是了解介质性质变化的重要方法之一(王伟涛,王宝善,2011; Xu,Song,2009)。

为了通过地震波了解地壳介质性质变化,在云南、甘肃省和新疆地区建立了大型陆地水体气枪震源观测基地(陈颙等,2017; 王伟涛等,2017),并积累了多年的观测数据。如何精确测量地震波速变化就成为了解地壳介质性质变化及其过程的关键问题之一(王彬等,2012; 王宝善等,2011,2016; Chen et al,2014; Liu et al,2019; Su et al,2019; Wang,Chen,2019)。

地震波速变化测量本质上是要解决数字信号处理的时延估计问题,关键是要提高时延估计的精度(刘自凤等,2015; 叶泵等,2017; 栾奕等,2016; 向涯等,2017; 林建民等,2010; Knapp,Carter,2003)。对于确定性的周期或非周期信号、或平稳随机噪声,增加取样长度有利于提高时延估计的精度; 但是对于非平稳随机信号,增加取样长度将带来新的误差。尽管大容量气枪震源具有较高的重复性,而实际地震记录信号由于存在随机干扰和能量衰减,具有一定的非平稳性。同时气枪震源能量有限,传播距离较远时信号信噪比低,而传播距离较近的台站震相波列发育较短,有效信号较短,这些因素决定了进行相关计算时只能进行较短取样(蒋生淼等,2017; 李孝宾等,2016,2017; 徐逸鹤等,2016)。

鉴于互相关的时延检测方法受相关窗口长度和信号信噪比影响较大,本文将改进的线性调频Z变换(Modified Chirp Z Transform,简称MCZT)谱细化和相关峰精确插值(Fine Interpolation of Correlation Peak,简称FICP)引入到陆地大容量气枪主动源观测中的波速干涉测量中,介绍了一种能够获得高精度走时差的短时窗时延估计方法。

陆地大型水体大容量气枪的能量释放由压力脉冲和气泡脉冲组成,震源持续时间约为2 s,不是典型的脉冲源,同时存在水位和水深等影响因素。为了消除或减弱这些影响因素,通常对接收台站记录与源附近的台站(参考台)记录做反褶积运算(王宝善等,2012; 翟秋实等,2016; 陈佳等,2017),得到接收台近似于脉冲震源的记录,即格林函数。

目前陆地大型水体大容量气枪的波速测量,多数采用处理流程和方法如下(王宝善等,2012; Wang et al,2019):

(1)参考格林函数选取:将接收台站的格林函数叠加得到参考格林函数或者选取某一次激发的台站格林函数作为参考函数。

(2)当前格林函数选取:将每一次激发的格林函数或者相邻几次激发的格林函数叠加形成信噪比更高的格林函数,作为该段时间(内)的格林函数。

(3)震相选取:根据震源距和传播特征选取具有明确物理意义的震相,用于干涉测量,通常选取P波和S波震相。

(4)干涉测量:将当前格林函数与参考格林函数的震相信号做互相关进行时延估计,得到时延量(即走时差),然后将走时差换算成波速变化率。也可用伸缩法(行鸿彦,唐娟,2008)直接得到波速变化率。

(5)波速变化率:根据震相走时将走时差转换为波速变化率。

地震波速变化测量常用方法有互相关法、伸缩法和尾波干涉法(王宝善等,2016; Alexandre,2014)。这些方法都是基于互相关运算,相关运算可以在时域或频域进行,本文主要讨论在频域计算中如何提高计算精度。而在频域进行互相关计算时,需要较精确的频谱估计。

标准快速傅里叶变换(Fast Fourier Transform,简称FFT)的计算速度快,但对于实数序列的FFT计算存在2个显著的缺陷:一是在计算结果中,序列左半边代表信号的频谱,右半边是左半边的共扼对称,只有一半是有用信息,浪费了另一半; 二是FFT计算出的信号的频谱相当于对频谱的离散抽样,存在栅栏效应。谱分辨率受取样长度t_p制约,Δf=1/t_p。在满足信号平稳性的前提下,为了减小栅栏效应的影响,可以增加取样长度,但这也增加了计算量。另外,标准FFT计算的相关函数,是从N点序列到N点序列的计算,相关函数序列的点的时间间隔与采样周期一致,会降低相关波形的分辨率(间隔增加),难以得到精确的时延。针对此问题,传统的方法是对相关峰进行插值计算找峰值点,因插值结果不能完全反映相关峰的形状,因此,插值计算会带来新的误差,且增加了运算量(杨亦春等,2002)。综上认为,精确的互谱估计和互谱插值对时延的精确估计很重要。

在大容量气枪源观测中,地震仪的采样率一般为100 Hz,实际时延量远小于一个采样间隔,互相关得到的时延量只有一个采样间隔的分辨率。为了提高分辨率,对互相关函数的相关峰通过余弦插值等方法进行插值,可以得到更精确的时延估计值,取得较好的结果(Wang et al,2019)。但余弦插值计算过程与信号的信息无关,即与互谱无关,插值结果只能近似反映相关峰的形状。因此要达到计算精度更高而计算量小的目的,有必要寻求同时提高互谱和相关峰计算精度的方法。为此,本文探讨利用MCZT谱细化和FCIP在陆地大容量气枪主动源观测中的波速干涉测量的应用,并与余弦插值方法进行了对比。

频谱的细化计算常见的方法有ZOOM算法、频率抽取法、降采样法、ZFFT(降采样)、cFFT(级联FF)等(杨亦春等,2002; 侯朝焕等,1990),这几种方法的共同特点是计算量大。MCZT算法可在不增加计算量的前提下提高频谱计算精度(马少春等,2014; 杨亦春等,2003; 唐娟,行鸿彦,2007; 杜娟,程擂,2010; 王莉等,2011; 徐世友,2003; Benesty et al,2004)。

对长度为N的时域序列x₁(n)和x₂(n),其相关函数为:

R(m)= E[x₁(n)·x₂(n+m)](1)

式中:R表示相关函数; m表示时延; E表示数学期望; n表示时间。

相关函数的主峰对应的时间就是信号x₁(n)和x₂(n)之间的时延量, R(n)的频谱R(k)与x₁(n)和x₂(n)的频谱X₁(k)和X₂(k)的关系为:

R(k)=X₁(k)·X^*₂(k)(2)

故频谱计算精度直接影响相关函数的精度。

MCZT定义为:

X(k)=MCZT[x(n)]=∑^N-1_n=0x(n)e^{-(j2π)/(N₁)kn}

(k= 0,1,…,N-1)(3)

定义MCZT的逆变换为IMCZT:

x'(n)=IMCZT[X(k)]=∑^N-1_n=0X(k)e^{-(j2π)/(N₁)kn}

(n=0,1,…,N-1)(4)

由式(3)(4)可看出,x'(n)不等于x(n),IMCZT与MCZT的计算结果在数值上是共轭关系,MCZT是CZT(线性调频Z变换)的一个特例。频谱的间隔由N₁决定,Δf=f_s/N,f_s为信号采样率,当N₁=f_s时,频谱的分辨率就达到1 Hz,N₁的选择不受N的限制,可以采用FFT实现快速计算。采用MCZT计算出的信号频谱,是从0 Hz开始到f_max=f_s*(N/N₁)Hz,即是截取了一段低频谱,为保证不影响相关精度,应该覆盖信号的主要频率范围(Sorensen,Burrus,1993),该方法适用于对低频信号分析。

主动源台阵接收到的主要信号频率为2～7 Hz,相对于100 Hz的采样率,MCZT不仅能计算出了被FFT丢失的谱,而且谱峰的高度也更加精确。另外由于MCZT只计算了有用信息所在的频段,因此还具有抑制宽带噪声的作用。

FICP最基本的思想是对2个信号的互谱补零,从而提高相关函数的分辨率。从采样定理可知,采样信号的频谱相当于原始信号频谱的周期重复,且重复间隔等于采样频率。离散傅里叶变换(Discrete Fourier Transfrom,简称DFT)计算出的信号频谱,其右半部分是左半部分的共扼对称,总长度对应采样率数值。提高采样率并增加序列点数(t_p不增加)不会改变频谱的形态,只会拉开左右两段谱之间的间隔。因此,如果把计算出的频谱中间加零,拉开左右距离,逆变换计算出的时域波形的采样率就会提高(侯朝焕等,1990)。

对于2个信号x₁(n), x₂(n)(n=0,1,…,N-1),加窗函数W(n)(n=0,1,…,N-1)后,用MCZT分别计算其细化频谱,根据相关定理计算得到2个信号的互谱:

R(k)=X₁(k)·X₂(k)(k=0,1,…,N-1)(5a)

式中,R(k)只是谱的前一部分,通过谱的共轭对称性,对R₁(x)补零扩展为N₁点R(k):

R(k)={R₁(k)(k=0,1,…,N-1)

0(k=N,N+1,…,N₁-n)

R^*₁(N₁-1-k)

(k=N₁-N+1,N₁-N+2,…,N₁-1)(5b)

通常信号的时延处在有限范围内,相关函数的主峰处在零值附近,因此只需要计算相关波形的左边(后段)有限点和右边(前段)有限点。实现这一算法的方法可以采用有限点时域(或频域)序列计算在有限点频域序列(或时域)的算法(Sorensen,Burrus,1993)该方法虽然考虑到了减少计算量的问题,但是不能对变换结果进行插值运算。在计算时延时,只需要计算相关峰附近的点,为此可以对逆变换计算进行简化。

按照式(5b)将互谱拉开,设定互谱的长度为N₂, 对于n=0,1,…,N-1的时域范围,相关函数前段为:

R_n1(n)=1/(N₂)∑^N-1_k=0R(k)e^{(j2π)/(N₂)kn}+1/(N₂)∑^N₂-1_k=N2-NR(k)e ^{j(2π)/(N₂)kn}

(6)

将式(6)两项分别计算可以发现,这两项分别为IMCZT计算过程:

R_n1(n)=1/(N₂)[A(n)+e^{-j(2π)/(N₂)nN}·B(n)](7)

式中:

A(n)=IMCZT[R₁(k)]=∑^N-1_k=0R₁(k)e ^{j(2π)/(N₂)kn}(8)

B(n)=IMCZT[R₃(k)](9)

式中:

R₃(k)={R^*₁(N-k)(k=1,2,…,N-1)

0(k=0)(10)

同理,对于n=N₂-N,N₂-N+1,…,N₂-1的时域范围, 取n=N₂-N+n₁, 对应的n₁=0,1,…,N-1相关函数为:

R_n2(n)=(E(n)+f(n))/(N₂)(n=0,1,…,N-1)(11)

式中:

E(n)=∑^N-1_k=0R₁(k)e ^{j(2π)/(N₂)k(N₂-N+n₁)}=IMCZT[H(k)]

(n=N₂-N+n₁; n₁=0,1,…,N-1; k=0,1; …,N-1)(12)

式中:

H(k)=R₁(k)e ^{j(2π)/(N₂)kN}(k=0,1,…,N-1)(13)

F(n)=G(n)·IMCZT[M(k)](14)

G(n)=e ^{j(2π)/(N₂)(N2-nN)}(15)

M(k)=R₃(k)e ^{j(2π)/(N₂)kN}(n = 0,1,…,N-1; k = 0,1,…,N-1)(16)

由此可组合出相关波的峰为:

R_n(n)={R_n2(N+n)(n=-N,…,-1)

R_n1(n)(n=0,1,…,N-1)(17)

通过设置N₂,可以使相关函数的峰值分辨率提高N₂/N₁倍。

一方面,提高互谱的频谱分辨率可以使相关波形更光滑,相关峰顶更圆,且频谱的细化不受取样长度的影响; 另一方面,提高插值的倍数,可以使相关峰顶分辨率提高。可见采用FICP算法可以提高时延估计的精度,特别是对于低采样的信号,其优点更明显。

本文根据地下介质波速变化测量精度的需求和宾川气枪主动源数据特点,介绍了一种能够在短时窗互相关计算中获得高精度时延估计的方法,通过实验对比分析,得到以下结论:

(1)基于MCZT和FICP的时延估计方法可以得到较高精度的时延估计值,通过模拟实验表明,该方法较余弦插值法有一定的提升,明显优于伸缩法。

(2)本文方法在短时窗互相关时差计算中优势较大,且在走时差计算过程中,通过使用MCZT,只计算有用信号所在频段,实现对噪声干扰的抑制,提高了地震信号时延估计的精度。

(3)利用FICP相关峰精确插值方法,克服了传统方法插值结果不能完全反映相关峰形状、插值计算而带来信号误差、增加运算量等问题,提高了相关峰顶的分辨率,特别是对低采样信号优势明显。

(4)实际数据计算表明,本文方法略优于余弦插值方法,也具有一定的抗噪能力。