CN103901470B - 一种海底沉积层反射强度及时延估计方法 - Google Patents

Abstract

本发明属于信号处理领域，具体涉及一种对白噪声具有较强的抑制作用的海底沉积层反射强度及时延估计方法。本发明包括：水听器输出信号首先做包括滤波、放大和采集自动增益控制的信号调理；经调理的信号分别与信号下边频的正、余弦信号相乘并通过带通滤波器滤除直流偏置；根据确知信号参数对滤波输出信号做最优分数阶Fourier变换，同时获取时间量纲化过程引入的时延偏移量；对最优分数阶Fourier变换输出信号做时间量纲化变换；由量纲化变换结果得到沉积层反射强度序列，同时通过峰值检测获取各沉积层的时延估计信息。本发明计算量小且鲁棒性高，可以工程实时实现，对白噪声具有很强的抑制作用。

Description

一种海底沉积层反射强度及时延估计方法

技术领域

本发明属于信号处理领域，具体涉及一种对白噪声具有较强的抑制作用的海底沉积层反射强度及时延估计方法。

背景技术

随着海洋工程、海洋勘查和海洋科学研究的不断深入，海底浅地层探测技术逐渐得到了广泛应用。为了进行海底浅地层的剖面成像、对沉积层结构进行探测，通常采用反卷积的方法对海底回波信号进行处理，得到海底浅地层的反射强度和时延信息。首先是快速傅立叶变换被用于反卷积技术，然而当回波信号信噪比降低时，该方法变得非常不稳定。在以后的研究过程中，不断有时域(t域)和频域(f域)反卷积方法被提出，如时域的共轭梯度法、频域的AR预测滤波法、以及维纳滤波与AR预测滤波联合处理方法等。但是上述方法都存在一定的问题：共轭梯度法若迭代次数选取过多则会产生虚假的振荡结果，若迭代次数选取过少又得不到很好的结果；频域AR预测滤波法频率窗应选取信号信噪比较高的部分，从而扩展后的频谱才能比较平坦，否则算法不稳定，由于噪声的存在影响特征值的选取，导致迭代结果出现虚假信息。另外，从反卷积方法的发展历程来看，几乎所有的反卷积方法都需要进行迭代运算，导致运算量增大，并且对信噪比都有较高的要求，在低信噪比下算法都存在不稳定的问题。

随着浅地层探测设备商业化和产业化的发展，对海底浅地层探测技术的实时性和鲁棒性提出了更高的要求，因此，提出具有较强实时性和噪声抑制能力的可实时实现的浅地层信号处理方法意义重大。

分数阶Fourier变换(Fractional Fourier Transform，FrFT)是傅立叶变换的广义形式，特别适合于处理LFM信号，最优阶FrFT会使LFM信号在分数阶Fourier域（u域）出现能量聚集，形成冲激响应信号。由互易原理可知，u域LFM信号变换到t域可得到冲激响应信号。由于沉积层的冲激响应反映了海底沉积层的反射强度和各沉积层的时延信息，因此求得沉积层冲激响应看实现沉积层的反射强度及时延估计。另外，由于FrFT以LFM信号为正交基，对白噪声具有很强的抑制作用，且主要有FFT实现，计算量小。因此，基于FrFT技术的海底沉积层反射强度及时延估计可具有较强的实时性。

发明内容

本发明的目的在于提供一种提高信号处理的鲁棒性，且计算量更小得海底沉积层反射强度及时延估计方法。

本发明的目的是这样实现的：

（1）水听器输出信号首先做包括滤波、放大和采集自动增益控制的信号调理；

（2）经调理的信号分别与信号下边频的正、余弦信号相乘并通过带通滤波器滤除直流偏置；

（3）根据确知信号参数对滤波输出信号做最优分数阶Fourier变换，同时获取时间量纲化过程引入的时延偏移量；

（4）对最优分数阶Fourier变换输出信号做时间量纲化变换；

（5）由量纲化变换结果得到沉积层反射强度序列，同时通过峰值检测获取各沉积层的时延估计信息。

步骤(3)中的时延偏移量为：

${\mathrm{\ξ}}_{\mathrm{offest}}=\left|\right[\mathrm{\Δ}{f}_L-\frac{\mathrm{\ΔT}}{2}.\frac{1-\sin \left({\mathrm{\α}}_c\right)}{\cos \left({\mathrm{\α}}_c\right)}].\frac{S}{\mathrm{\Δk}}|$

式中f_L表示LFM信号的下边频，T为截取待处理信号的长度，α_c=arctan(k)为LFM信号时频分布与时间轴的夹角，k为LFM信号的频率变化率，Δ表示对应参量的量纲归一化值，S=sqrt(T/f_s)为量纲归一化时间因子。

步骤(4)中的时间量纲化处理过程中，无量纲的u域样本点序列U=1，2，...N变换为以时间为量纲的时域样本点序列Num_j（j=N_offset，N_offset+1，...N_offset+N）时，进行变换：

${\mathrm{Num}}_j=\left\{\right[\mathrm{\Δ}{f}_L-\frac{(U_i-N/2)}{S{f}_s\cos \left({\mathrm{\α}}_c\right)}+\frac{\mathrm{\ΔT}}{2}\tan \left({\mathrm{\α}}_c\right)]/\mathrm{\Δk}+\mathrm{\Δ}{\mathrm{\ξ}}_{\mathrm{offest}}\}S.f_{s}i=\mathrm{1,2}...N$

其中，N为处理信号样本的点数，N_offset＝ζ_offset·f_s为时间量纲化过程引入偏移量所对应的点数。

本发明的有益效果在于：计算量小且鲁棒性高，可以工程实时实现，对白噪声具有很强的抑制作用。本发明具有广泛的应用前景，虽然针对海底浅地层探测提出，但可以推广应用于其它任何利用LFM信号作为探测信号的目标探测领域。

附图说明

图1是海底沉积层反射强度及时延估计方法实施流程图；

图2是LFM信号在t-f平面及最优u-v平面上的分布；

图3是u域无量纲量变换为t域量纲时存在的偏移量；

图4是时延LFM信号及其FrFT；

图5是海底沉积层冲激响应及回被信号波形；

图6是海底沉积层反射强度序列；

图7是SNR为6dB时回波信号及沉积层反射强度序列；

图8是SNR为3dB时本发明方法与AR预测滤波方法的对比。

具体实施方式

下面结合附图对本发明的具体实施方式作更详细的描述：

水听器输出信号首先做包括滤波、放大和采集自动增益控制的信号调理；

经调理的信号分别与信号下边频的正、余弦信号相乘并通过带通滤波器滤除直流偏置；

根据确知信号参数对滤波输出信号做最优分数阶Fourier变换，同时获取时间量纲化过程引入的时延偏移量；

对最优分数阶Fourier变换输出信号做时间量纲化变换；

由量纲化变换结果得到沉积层反射强度序列，同时通过峰值检测获取各沉积层的时延估计信息。

为了实现本发明的目的，将LFM信号作为探测信号，充分利用FrFT技术对LFM信号的优越处理性能，提出无量纲u域至时间量纲t域变换的量纲化过程，实现基于FrFT的海底沉积层反射强度及时延实时估计。

本发明涉及一种针对LFM信号的海底沉积层反射强度及时延估计方法，该方法实施的流程如图1所示，整个流程分解为以下5个步骤完成：

步骤(1)：水听器输出信号首先做包括滤波、放大和采集自动增益控制的信号调理；

步骤(2)：经调理的信号分别与信号下边频的正、余弦信号相乘并通过带通滤波器滤除直流偏置；

步骤(3)：根据确知信号参数对滤波输出信号做最优分数阶Fourier变换，同时获取时间量纲化过程引入的时延偏移量；

步骤(4)：对最优分数阶Fourier变换输出信号做时间量纲化变换；

步骤(5)：由量纲化变换结果得到沉积层反射强度序列，同时通过峰值检测获取各沉积层的时延估计信息。

为详细说明本发明的机理，以具体实例加以说明：

本发明中采用LFM信号作为探测信号，信号可表示为：

其中，A为信号幅度，f₀为中心频率，k为频率变化率，T为信号脉冲宽度。

由射线声学原理可知，发射信号经不同海底浅地层反射后到达接收点，最终接收信号是各反射信号叠加的结果。海底沉积层（信道）的冲激响应为：

$h\left(t\right)=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{A}_i\mathrm{\δ}(t-{\mathrm{\ι}}_i)---\left(2\right)$

其中，A_i为第i个沉积层反射回波信号的幅度，τ_i为第i个沉积层的时延，N为海底沉积层的个数。从而海底回波信号可表示为：

y(t)＝s(t)*h(t) (3)

信号x(t)的FrFT为：

$X_{\mathrm{\α}}\left(u\right)=\left\{F^{\mathrm{\α}}\right[x\left(t\right)\left]\right\}\left(u\right)={\∫}_{-\∞}^{\∞}x\left(t\right)K_{\mathrm{\α}}(t,u)\mathrm{dt}---\left(4\right)$

式中α＝pπ2是FrFT运算对应的旋转角度。FrFT的变换核K_α(t,u)表示为：

$K_p(t,u)=\left\{\begin{array}{cc}\sqrt{\frac{1-j\cos \mathrm{t\α}}{2\mathrm{\π}}}\exp (j\frac{t^2+u^2}{2}\cot \mathrm{\α}-\mathrm{jtu}\csc ),& \mathrm{\α}\≠\mathrm{n\π}\\ \mathrm{\δ}(t-u),& \mathrm{\α}=2\mathrm{n\π}\\ \mathrm{\δ}(t+u),& \mathrm{\α}=(2n\±1)\mathrm{\π}\end{array}\right.---\left(2\right)$

而p阶FrFT的逆变换被定义为：

$x\left(t\right)={\∫}_{-\∞}^{+\∞}{X}_p\left(u\right)K_{-p}(t,u)\mathrm{du}---\left(3\right)$

与LFM信号FrFT最大幅度或WD直线积分最大值(或最小带宽——分数阶域带宽)相对应的阶数与旋转角度表示为：

其中Cr为归一化频率变化率(Normalized LFM rate)，表示为：

C_r＝ΔB/ΔT (6)

式中ΔB为信号的归一化频带宽度(对采样频率的尺度归一化)，ΔT为信号归一化脉宽(对观测时间窗长度的尺度归一化)。

FrFT是广义的Fourier变换，而LFM信号在最优的u域上具有很强的聚集性，这一特点有利于FrFT分析与处理LFM信号。FrFT的处理结果取决于参量α，它可以被看作信号逆时针旋转任意角度α时的信号表示。根据FrFT的定义以及性质，可得旋转角度为α_opt时LFM信号在时频平面上分布的几何关系如图2所示，U₀即为LFM信号在u域的分布位置。

由量纲归一化和LFM信号在时频域分布位置的几何关系，推得任意时延τ时LFM信号 在u域分布位置ΔU_τ的表达式为：

$\mathrm{\&Delta;}{U}_{\mathrm{\&iota;}}=[\mathrm{\&Delta;}{f}_L-\mathrm{\&Delta;\&iota;}.\mathrm{\&Delta;k}].\cos \left({\mathrm{\&alpha;}}_c\right)+\frac{\mathrm{\&Delta;T}}{2}.\sin \left({\mathrm{\&alpha;}}_c\right)\&le;\mathrm{\&Delta;\&iota;}<\mathrm{\&Delta;T}---\left(7\right)$

式中f_L表示LFM信号的下限频率，τ为LFM信号分量的时延，k为LFM信号的频率变化率，T为截取待处理信号的长度，f_s为采样频率，α_c=arctan(k)为LFM信号时频分布与时间轴的夹角，式中“Δ”表示对应参量的量纲归一化值。表示成离散序列的形式，LFM信号在u域分布位置为U_τ=ΔU_τSf_s，S=sqrt(T/f_s)为量纲归一化时间因子。

在海底浅地层探测中，为了直观地生成海底浅地层剖面图，总是希望得到t域海底浅地层反射强度序列，因此将u域信息转换至t域是一个时关键性的问题。为了描述u域上任意样本点对应t域信号的时延量，从而将无量纲的u域参量变换至t域，需要得知u域样本点与t域对应样本点间存在的内在关系，而这一过程即为无量纲u域至t域的量纲化过程。由FrFT的特性可知，在时间量纲化过程中，u域与t域相应采样点之间会存在一个偏移量。由量纲归一化以及FrFT的定义推得，该偏移量与u域信号最后一个样本点映射到t域坐标的位置相对应，推得t域上的该偏移量ζ_offset表示为：

${\mathrm{\&xi;}}_{\mathrm{offest}}=\left|\right[\mathrm{\&Delta;}{f}_L-\frac{\mathrm{\&Delta;T}}{2}.\frac{1-\sin \left({\mathrm{\&alpha;}}_c\right)}{\cos \left({\mathrm{\&alpha;}}_c\right)}].\frac{S}{\mathrm{\&Delta;k}}|---\left(8\right)$

t-f平面上u域信息变换至t域时存在偏移量的示意图如图3所示。在计算得到偏移量ζ_offset的前提下即可将信号的u域信息量纲化变换到t域，最终得到与之对应的具有时间量纲的信号信息。

无量纲的u域样本点序列U=1,2,...N变换为以时间为量纲的t域样本点序列Num_j（j=N_offset,N_offset+1,...N_offset+N）时存在以下关系：

${\mathrm{Num}}_j=\left\{\right[\mathrm{\&Delta;}{f}_L-\frac{(U_i-N/2)}{S{f}_s\cos \left({\mathrm{\&alpha;}}_c\right)}+\frac{\mathrm{\&Delta;T}}{2}\tan \left({\mathrm{\&alpha;}}_c\right)]/\mathrm{\&Delta;k}+\mathrm{\&Delta;}{\mathrm{\&xi;}}_{\mathrm{offest}}\}S.f_{s}i=\mathrm{1,2}...N---\left(9\right)$

其中，N为处理信号样本点数，N_offset为时间量纲化过程引入偏移量所对应的点数。上式即实现了u域无量纲量至t域时间量纲量的量纲化变换过程。从而可直接得到海底沉积层反射强度序列，通过峰值检测容易的得到各沉积层的时延信息。

下面以具体仿真实例给出本发明方法的整个实现过程：

LFM信号时延为5ms，脉宽为6ms，频带范围为1～6kHz，采样频率为40kHz时的波形及其最优阶FrFT如图4所示，截取待处理信号长度为36ms。利用公式(7)求得LFM信号在u域的分布位置为341.25，与图4的仿真结果相吻合，存在的小数偏差是由于FrFT输出信号的横坐标取值为整数引起的。

图5依次给出了依据公式（2）和公式（3）得到的海底沉积层的冲激响应及其回波信号，其中：发射信号参数与上述参数相同，各沉积层反射界面对应回波时刻分别为：3、7、9.5、12、14、16、18、21、23、25.5、28ms。首先，对海底回波信号进行最优阶FrFT变换；然后对u域至t域变换时的偏移量进行计算，用于时间量纲化变换时对时延信息进行修正；最后，采用本发明提出的时间量纲化变换，将回波信号的u域表示转换至时间域，从而得到海底沉积层的反射强度序列，如图6所示。由图6可知，在不存在噪声的情况下可以准确的对海底沉积层各反射层的反射强度和时延信息进行估计。

当引入噪声时，对SNR=6dB的回波信号进行处理，回波信号的t域波形及利用本发明方法处理得到的海底沉积层反射强度序列如图7所示。由于受到噪声干扰的影响，在时频域上已不能清晰地判断出各回波信号分量所处的位置（如图7（a）所示），但利用本发明提出的基于FrFT的海底沉积层反射强度及时延估计方法计算得到的沉积层反射强度序列如图7(b)所示，与图6相比，除极个别小反射强度沉积层信息被噪声掩盖外，仍能精确估计出绝大多数主要海底沉积层界面的反射强度及时延信息。

为了进一步说明本发明的处理性能，将其与AR预测滤波处理方法进行对比。AR预测滤波方法已被Satchidananda Panda成功应用于海底浅地层探测，并得到了较好的处理效果。图8给出了SNR=3dB时采用两种方法的处理结果，其中：图8(a)为本发明方法的处理结果，图8(b)为采用AR预测滤波方法处理得到的海底沉积层反射强度。通过对两处理结果进行对比发现，本发明方法仍能够很好地估计出海底浅地层的反射强度，并且估计的更深沉积层反射界面的反射强度失真较小，而AR预测滤波方法虽然获得了大多数反射界面的反射强度，但后续反射界面反射强度的失真较大，且处理性能不稳定。

本发明方法主要由两次FFT运算组成，计算量小可以实时实现，而AR预测滤波方法需经过多次的FFT及迭代运算才能达到相应的处理效果，计算量大，不易实时实现，且受噪声影响较大。经过多次实验表明，本发明方法不仅稳定性优于AR预测滤波方法，而且在低信噪比下仍能得到很好的处理效果，具有更好的鲁棒性和抗噪声性能。

Claims (1)

1.一种海底沉积层反射强度及时延估计方法，其特征在于：

(1)水听器输出信号首先做包括滤波、放大和采集自动增益控制的信号调理；

(2)经调理的信号分别与信号下边频的正、余弦信号相乘并通过带通滤波器滤除直流偏置；

(3)根据确知信号参数对滤波输出信号做最优分数阶Fourier变换，同时获取时间量纲化过程引入的时延偏移量；

(4)对最优分数阶Fourier变换输出信号做时间量纲化变换；

(5)由量纲化变换结果得到沉积层反射强度序列，同时通过峰值检测获取各沉积层的时延估计信息；

所述步骤(3)中的时延偏移量为：

${\mathrm{\&zeta;}}_{offset}=|\&lsqb;{\mathrm{\&Delta;f}}_L-\frac{\mathrm{\&Delta;}T}{2}\&CenterDot;\frac{1-sin\left({\mathrm{\&alpha;}}_c\right)}{cos\left({\mathrm{\&alpha;}}_c\right)}\&rsqb;\&CenterDot;\frac{S}{\mathrm{\&Delta;}k}|$

式中f_L表示LFM信号的下边频，T为截取待处理信号的长度，α_c＝arctan(k)为LFM信号时频分布与时间轴的夹角，k为LFM信号的频率变化率，Δ表示对应参量的量纲归一化值，S＝sqrt(T/f_s)为量纲归一化时间因子；

所述步骤(4)中的时间量纲化处理过程中，无量纲的u域样本点序列U＝1，2，…N变换为以时间为量纲的时域样本点序列Num_j(j＝N_offset，N_offset+1，…N_offset+N)时，进行变换：

${\mathrm{Num}}_j=\{\&lsqb;{\mathrm{\&Delta;f}}_L-\frac{(U_i-N/2)}{{\mathrm{Sf}}_s\cos \left({\mathrm{\&alpha;}}_c\right)}+\frac{\mathrm{\&Delta;}T}{2}tan\left({\mathrm{\&alpha;}}_c\right)\&rsqb;/\mathrm{\&Delta;}k+{\mathrm{\&Delta;\&zeta;}}_{offset}\}S\&CenterDot;f_s,i=1,2...N$

其中，N为处理信号样本的点数，N_offset＝ζ_offset·f_s为时间量纲化过程引入偏移量所对应的点数，f_s为采样频率。