一种双线性相参的立方相位信号参数估计方法与流程

本发明涉及信号参数估计技术领域，尤其涉及一种双线性相参的立方相位信号参数估计方法。

背景技术：

立方相位信号是一种瞬时频率随着时间而变化的非平稳信号，其普遍存在于自然界和人工应用中，在雷达、声纳、无线通信、地质勘探和医学成像等工程和科学领域有着广泛的应用。立方相位信号参数包含着重要信息，例如：在雷达信号处理中，中心频率、调频率和二次调频率分别对应目标的初始速度、加速度和加速度变化，因此，立方相位信号参数估计的准确性在工程应用中至关重要。

目前，立方相位信号参数估计方法大致可以分为三类：线性方法、多重线性方法和双线性方法。典型的线性方法包括最大似然估计和离散调频率傅里叶变换，其可以获得较高的抗噪声性能，且不存在交叉项，然而较低的分辨率和较高的计算复杂度降低了其实用性。多重线性方法通过多阶相关函数降低立方相位信号阶数，从而实现分辨率的提高和计算复杂度的降低，典型的多重线性方法包括CIGCPF(coherently integrated generalized cubic phase function)、SCFT(scaled Fourier transform)-based algorithm、CRQCRD(chirp rate-quadratic chirp rate distribution)和GDT(generalized decoupling technique)-based algorithm。然而，多阶相关函数严重影响交叉项抑制性能和抗噪性能。双线性方法特点介于线性方法和多重线性方法之间，典型的双线性方法包括CPF(cubic phase function)和(LPWD)local polynomial Wigner distribution，然而立方相位信号中二次调频率的存在导致现有的双线性方法仅能够完成一维的能量积累，因此交叉项干扰和低抗噪性能问题无法解决。

技术实现要素：

针对上述已有技术的不足，本发明的目的在于提供一种双线性相参的立方相位信号参数估计方法，克服了已有技术无法兼顾交叉项抑制性能、抗噪性能、分辨率和计算复杂度的缺陷，能够提高系统的估计精度和处理效率。

为达到上述目的，本发明的实施例采用如下技术方案予以实现。

一种双线性相参的立方相位信号参数估计方法，所述方法包括如下步骤：

步骤1，采集一段含有噪声的信号，用离散多分量立方相位信号z(n)表示，并计算所述离散多分量立方相位信号z(n)的总能量；n表示离散采样时间变量，N为离散信号长度；所述离散多分量立方相位信号z(n)包含P个信号分量；

步骤2，计算所述离散多分量立方相位信号z(n)的参数化双线性瞬时自相关函数R_z(n，m)；m表示离散延时变量；

步骤3，构造参数化双线性瞬时自相关函数R_z(n，m)的时间-调频率分布函数G_z(n，f_m)，并对所述时间-调频率分布函数G_z(n，f_m)进行复数求模操作，并得到模值的整数次幂函数K_z(n，f_m)；f_m表示离散频域变量；

步骤4，对所述模值的整数次幂函数K_z(n，f_m)沿f_m轴做逆傅里叶变换，得到逆傅里叶变换后的函数Q_z(n，b)，再对逆傅里叶变换后的函数Q_z(n，b)采用变尺度傅里叶变换完成沿离散采样时间轴的能量相参积累，得到能量相参积累函数Γ_z(f_[nb]，b)；b表示离散频域变量f_m对应的离散时域变量，f_[nb]表示离散采样时间变量n对应的变尺度频域变量；

步骤5，沿能量相参积累函数Γ_z(f_[nb]，b)的时间轴b做傅里叶变换，得到参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)，f_b表示离散时域变量b对应的频域变量；

步骤6，确定所述参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)的最大峰值，得到最大峰值对应的第j个信号分量的调频率估计值和二次调频率估计值，j＝1，…，P；

步骤7，根据第j个信号分量的调频率估计值和二次调频率估计值对离散多分量立方相位信号z(n)进行多普勒扩散补偿，并确定最大峰值对应的第j个信号分量的中心频率估计值和幅度估计值；

步骤8，根据第j个信号分量的调频率估计值、二次调频率估计值、中心频率估计值以及幅度估计值重构第j个信号分量；

步骤9，从离散多分量立方相位信号z(n)中滤除第j个信号分量，得到剩余离散信号，并计算剩余离散信号的能量；

步骤10，若剩余离散信号的能量与离散信号的总能量的比值大于或者等于预设门限，则将离散信号替换为剩余离散信号，并依次重复执行步骤2至9，直到剩余离散信号的能量与离散信号的总能量的比值小于预设门限时，所述立方相位信号的参数估计结束。

本发明的有益效果为：本发明属于双线性方法，然而相比现有的双线性方法，本发明实现了二维能量积累，在交叉项抑制和抗噪性能方面有较大的改善。由于本发明利用复数取模值操作消除相位中的高次项，因此二维能量积累的实现并不需要较高的计算复杂度。与线性方法相比，本发明分辨率明显改善且计算复杂度降低，与多重线性方法相比，不存在严重的交叉项干扰和抗噪性能降低。因此，本发明克服了已有技术无法兼顾分辨率、计算复杂度、抗噪性能和交叉项抑制性能的缺陷，益于提高系统的估计精度和处理效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种双线性相参的立方相位信号参数估计方法的流程示意图；

图2为仿真实验1中多分量立方相位信号z_s1(n)在二次调频率-调频率-信号幅度值三维空间的分布示意图；

图3为仿真实验2中多分量立方相位信号z_s2(n)在二次调频率-调频率-信号幅度值三维空间的分布示意图；

图4为仿真实验3中立方相位信号z_s3(n)调频率的估计均方误差曲线示意图；

图5为仿真实验3中立方相位信号z_s3(n)二次调频率的估计均方误差曲线示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种双线性相参的立方相位信号参数估计方法，参照图1，所述方法包括如下步骤：

步骤1，采集一段含有噪声的信号，用离散多分量立方相位信号z(n)表示，，并计算所述离散多分量立方相位信号z(_n)的总能量；n表示离散采样时间变量，N为离散信号长度；所述离散多分量立方相位信号z(n)包含P个信号分量。

通过信号采集系统获取一段含有实际噪声的雷达信号、脑电波信号或地震波信号，步骤1具体为：

(1a)采用离散多分量立方相位信号将采集信号近似表示为离散信号z(n)，其表示式为：

其中，n表示离散采样时间变量，N为离散信号长度，T_z表示采样时间间隔，p表示第p个信号分量，P表示离散信号中的信号分量总数，σ_p表示第p个信号分量的幅度、表示第p个信号分量的中心频率、表示第p个信号分量的调频率、表示第p个信号分量的二次调频率；w(n)表示离散信号所含噪声；

(1b)计算所述离散多分量立方相位信号z(n)的总能量其中，|·|表示复数取模值操作。

步骤2，计算所述离散多分量立方相位信号z(n)的参数化双线性瞬时自相关函数R_z(n，m)；m表示离散延时变量。

步骤2具体为：

构造离散信号z(n)的时延函数z(n+m+c)和z(n-m-c)，并计算离散信号z(n)的参数化双线性瞬时自相关函数R_z(n，m)：

其中，m表示离散延时变量，M为离散延时变量长度，c表示离散延时常量，R_cros，n，z(n，m)表示交叉项和噪声。

步骤3，构造参数化双线性瞬时自相关函数R_z(n，m)的时间-调频率分布函数G_z(n，f_m)，并对所述时间-调频率分布函数G_z(n，f_m)进行复数求模操作，并得到模值的整数次幂函数K_z(n，f_m)；f_m表示离散频域变量。

步骤3具体为：

(3a)沿参数化双线性瞬时自相关函数R_z(n，m)的m轴做傅里叶变换，构造参数化双线性瞬时自相关函数R_z(n，m)的时间-调频率分布分布函数G_z(n，f_m)：

其中，f_m表示离散频域变量，F_m为变量f_m离散间隔，δ(·)表示冲激响应函数，G_cros，n，z(n，f_m)表示交叉项和噪声；

(3b)对时间-调频率分布函数G_z(n，f_m)进行复数求模操作，并得到模值的整数次幂函数K_z(n，f_m)：

其中，|·|表示复数取模值操作，x表示幂指数，K_crosn，z(n，f_m)表示交叉项和噪声。

步骤4，对所述模值的整数次幂函数K_z(n，f_m)沿频率轴做逆傅里叶变换，得到逆傅里叶变换后的函数Q_z(n，b)，再对逆傅里叶变换后的函数Q_z(n，b)采用变尺度傅里叶变换完成沿离散采样时间轴的能量相参积累，得到能量相参积累函数Γ_z(f_[nb]，b)；b表示离散频域变量f_m对应的离散时域变量，f_[nb]表示离散采样时间变量n对应的变尺度频域变量。

步骤4具体为：

(4a)沿模值的整数次幂函数K_z(n，f_m)的频率轴做逆傅里叶变换，得到逆傅里叶变换后的函数Q_z(n，b)：

其中，b表示f_m对应的离散时域变量，Q_cros，n，z(n，b)表示交叉项和噪声；

(4b)采用变尺度傅里叶变换完成沿Q_z(n，b)离散采样时间轴的能量相参积累，得到能量相参积累函数Γ_z(f_[nb]，b)：

其中，f_[nb]表示离散采样时间变量n对应的变尺度频域变量，ξ表示常量缩放因子，r_cros，n，z(f_[nb]，b)表示交叉项和噪声。

步骤5，沿能量相参积累函数Γ_z(f_[nb]，b)的时间轴做傅里叶变换，得到参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)，f_b表示离散时域变量b对应的频域变量。

步骤5具体为：

沿能量相参积累函数Γ_z(f_[nb]，b)的时间轴b做傅里叶变换，得到参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)：

其中，f_b表示离散变量b对应的频域变量，Ψ_cros，n，z(f_[nb]，f_b)表示交叉项和噪声。

步骤6，确定所述参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)的最大峰值，得到最大峰值对应的第j个信号分量的调频率估计值和二次调频率估计值，j＝1，…，P。

步骤6具体为：

利用峰值检测技术检测参数化调频率-二次调频率分布函数Ψ_z(f_[nb]，f_b)的最大峰值，获取最大峰值对应的第j个信号分量的调频率估计值和二次调频率估计值

其中，表示第j个信号分量的调频率估计值，表示第j个信号分量的二次调频率估计值，f′_b和f′_[nb]表示第j个信号分量峰值所对应的位置坐标，arg max(·)表示峰值检测操作。

步骤7，根据最大峰值对应的第j个信号分量的调频率估计值和二次调频率估计值对离散信号z(n)进行多普勒扩散补偿，并确定最大峰值对应的第j个信号分量的中心频率估计值和幅度估计值。

步骤7具体为：

(7a)利用第j个信号分量的调频率估计值和二次调频率估计值对离散多分量立方相位信号z(n)进行多普勒扩散补偿，得到多普勒扩散补偿后的信号z_com(n)：

其中，z_com，rest(n)表示噪声和多普勒扩散未被补偿的信号；

(7b)针对多普勒扩散补偿后的信号z_com(n)进行傅里叶变换操作，得到

其中，f_n表示离散变量n对应的频域变量，z_{com，dft，rest}(f_n)表示噪声和由于多普勒扩散能量未被积累的信号；

(7c)利用峰值检测技术确定最大峰值对应的第j个信号分量的中心频率估计值和幅度估计值σ′_j：

其中，σ′_j第j个信号分量的幅度估计值，表示第j个信号分量的中心频率估计值，A′表示峰值幅度，f′_n表示峰值所对应的位置坐标。

步骤8，根据第j个信号分量的调频率估计值、二次调频率估计值、中心频率估计值以及幅度估计值重构第j个信号分量。

步骤8具体为：

根据第j个信号分量的调频率估计值二次调频率估计值中心频率估计值以及幅度估计值σ′_j重构第j个信号分量z_j(n)为：

其中，表示第j个信号分量的调频率估计值，表示第j个信号分量的二次调频率估计值，σ′_j表示第j个信号分量的幅度估计值，表示第j个信号分量的中心频率估计值。

步骤9，从离散多分量立方相位信号z(n)中滤除第j个信号分量，得到剩余离散信号，并计算剩余离散信号的能量。

步骤9具体为：

(9a)利用窄带带阻滤波器作用于z_com，dft(f_n)，将第j个信号分量的中心频率估计值附近的频谱滤除，得到频谱滤除后的信号z_eli(f_n)：

z_eli(f_n)＝z_com，dft(f_n)win_j(f_n)

其中，d的值根据左右的频谱宽度预先设定；

(9b)对频谱滤除后的信号z_eli(f_n)做逆傅里叶变换，并乘以从而获得滤除第j个信号分量的剩余离散信号z_lc(n)：

(9c)计算剩余离散信号的能量

步骤10，若剩余离散信号的能量与离散信号的总能量的比值大于或者等于预设门限，则将离散信号替换为剩余离散信号，并依次重复执行步骤2至9，直到剩余离散信号的能量与离散信号的总能量的比值小于预设门时，所述立方相位信号的参数估计结束。

下面结合仿真实验对本发明的效果做进一步说明。

仿真实验1：为验证本发明在多分量立方相位信号环境下有较好的外部交叉项抑制性能和高分辨率，通过采集系统采集获得无噪声(或噪声可以忽略不计)的多分量立方相位信号，其离散形式为：

其中，n＝-256，-255，-244，…244，255，σ₁＝1，σ₂＝1，σ₃＝1，在仿真实验1中，我们设定ξ＝1和c＝128，采用本发明对多分量立方相位信号z(n)进行参数估计。

参照图2，为仿真实验1中多分量立方相位信号z_s1(n)在二次调频率-调频率-信号幅度值三维空间的分布示意图。图2中，水平两个轴分别表示二次调频率(单位为Hz/s²)和调频率(单位为Hz/s)，竖直轴表示信号幅度值。通过图2可以看出，本发明很好的抑制了交叉项的干扰，三个分量能量得到很好的聚焦(高聚焦对应高分辨率)。通过峰值点位置检测既可以估计出原始信号二次调频率和调频率。表1为利用本发明得出的仿真实验1中立方相位信号z_s1(n)各分量对应的参数估计表。表1中，代表z_s1(n)各分量的幅度估计值，代表z_s1(n)各分量的中心频率估计值，代表z_s1(n)各分量的调频率估计值，代表z_s1(n)各分量的调频率估计值。由表1中的参数估计值可以看出，本发明在多分量立方相位信号环境下具有较高的参数估计精度。

表1三分量立方相位信号z_s1(n)参数估计值

仿真实验2：为验证本发明在噪声环境下对多分量立方相位信号参数有较高的估计精度，通过采集系统采集获得含有白噪声的多分量立方相位信号，信噪比为-10dB，其离散形式为：

其中，n＝-256，-255，-244，…244，255，σ₁＝1，σ₂＝1，σ₃＝1，在仿真实验2中，我们设定ξ＝1和c＝128，采用本发明对多分量立方相位信号z_s2(n)进行参数估计。

参照图3，为仿真实验2中多分量立方相位信号z_s2(n)在二次调频率-调频率-信号幅度值三维空间的分布示意图。图3中，水平两个轴分别表示二次调频率(单位为Hz/s²)和调频率(单位为Hz/s)，竖直轴表示信号幅度值。通过图2可以看出，在输入信噪比-10dB环境下，本发明同样能够很好的分离多分量立方相位信号，通过峰值点位置检测既可以估计出原始信号二次调频率和调频率。表2为利用本发明得出的仿真实验2中z_s2(n)各分量对应的参数估计表。表2中，代表z_s2(n)各分量的幅度估计值，代表z_s2(n)各分量的中心频率估计值，代表z_s2(n)各分量的调频率估计值，代表z_s2(n)各分量的调频率估计值。通过表2我们可以看出，本发明在噪声环境下对多分量立方相位信号参数有较高的估计精度。

表2-10dB信噪比环境下三分量立方相位信号z_s2(n)参数估计值

仿真实验3：为进一步验证本发明的高抗噪性能，通过采集系统获得含有白噪声的单分量立方相位信号，信噪比范围为[-15dB∶1dB∶0dB]，其离散形式为：

其中，n＝-256，-255，-244，…244，255，σ₁＝1，在仿真实验3中，我们设定ξ＝1和c＝128，利用蒙特卡洛方法获取本发明的抗噪声门限。

参照图4和图5，为仿真实验3中z_s3(n)调频率和二次调频率估计的均方误差。图3中，水平轴表示输入信噪比，竖直轴表示z_s3(n)调频率估计的均方误差。图4中，水平轴表示输入信噪比，竖直轴表示z_s3(n)二次调频率估计的均方误差。通过图3和图4可以看出，当信噪比大于等于-11dB时，z_s3(n)调频率和二次调频率估计的均方误差接近克拉美罗界。显然，相比已有立方相位信号参数估计方法，本发明有明显的抗噪性能优势。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。