一种基于功率谱的线性调频信号参数估值方法

摘要

本发明公开了一种基于功率谱的线性调频信号参数估值方法，该方法包括：检测待估值的线性调频信号功率谱P(x)；提取所述线性调频信号功率谱P(x)的形状特征CO(P(x))；计算所述形状特征上升沿的斜率β，并根据调频斜率k和所述上升沿斜率β的关系计算所述线性调频信号的调频斜率k；根据所述形状特征上升沿的峰值频率f

说明书

技术领域

本发明涉及一种基于功率谱形态学运算的线性调频信号参数估值方 法。

背景技术

线性调频信号是一种脉冲压缩信号，具有大时宽带宽积特性，抗干 扰性能良好，是目前最常用的脉冲压缩信号。在雷达系统中，线性调频 信号参数的准确估值是获取雷达的用途、掌握相关信息的前提条件。

线性调频信号是指持续时间内频率随时间连续线性变化的信号，线 性调频信号的时域表达式为：

其中：f_s为起始频率，B为带宽，τ为脉宽，k为调频斜率，B＝kτ。 线性调频脉冲信号的频率在脉冲持续时间内呈线性变化。

线性调频信号功率谱可所示为：

$P_{\mathrm{LFM}}\left(\omega \right)=\frac{A^2\mathrm{\pi \tau }}{B}\{{[C\left(v_1\right)+C\left(v_2\right)]}^2+{[S\left(v_1\right)+S\left(v_2\right)]}^2\}$

其中：

$v_1=\sqrt{D}\frac{1+(\omega -{\omega }_0)/\mathrm{\pi B}}{\sqrt{2}},$$v_2=\sqrt{D}\frac{1-(\omega -{\omega }_0)/\mathrm{\pi B}}{\sqrt{2}}$

$C\left(v\right)=\underset{0}{\overset{v}{\int }}\cos \left(\frac{\pi }{2}{x}^2\right)\mathrm{dx},$$S\left(v\right)=\underset{0}{\overset{v}{\int }}\sin \left(\frac{\pi }{2}{x}^2\right)\mathrm{dx}$

线性调频信号功率谱特点如下：

①功率谱关于中心频率对称；

②频带覆盖范围较宽，具有大的时宽带宽积特性；

③时宽带宽积越大，线性调频信号越接近矩形。

由于线性调频信号是一种非平稳信号，在短时间内的初始频率和终 止频率差距很大，所以它的参数估计问题较一般的信号复杂，因此成为 研究的难点。线性调频信号常用的参数如下所示：

①起始频率f_s、终止频率f_e；

②带宽B，其中B＝|f_s-f_e；

③脉宽τ，也就是时宽；

④调频斜率k。

现有的线性调频信号参数估值方法大都基于时域的信号分析，常用的 参数估计方法有最大似然估计、Radon-Wigner变换(RWT)、 Radon-Ambiguity变换法(RAT)和分数阶傅里叶变换(FRFT)等。2010年9 月15日公开的专利CN101833035A中提出了一种线性调频信号参数估计 方法及其实施装置，该方法利用能量重心插值估计和Radon模糊变换来 获得调频斜率的估计。2011年，罗蓬等人发表在《系统与电子工程》的 文献“高精度LFM信号参数估计的谱校正法”。2011年，宋军发表在《系 统与电子工程》的文献“LFM信号参数估计的插值FRFT方法”。2012 年6月20日公开的专利CN102508206A中提出了一种基于小波包去噪和 功率谱熵的线性调频信号参数估计方法，利用该方法可获得调频信号的 调频斜率。基于时频分析的线性调频信号参数估计的计算量很大，无法 解决估计精度与运算量之间的矛盾，而基于小波包去噪和功率谱熵的线 性调频信号参数估计方法需要计算接收信号的信噪比，增加了参数估计 的复杂度。

发明内容

本发明目的在于：实现未知先验信息条件下的线性调频信号参数的 估值，克服现有技术中低信噪比环境下线性调频信号估值结果受噪声影 响大而且计算量大的缺陷。

本发明的技术方案是：本发明提供了一种基于功率谱的线性调频信 号参数估值方法，其特征在于，该方法包括：

步骤1、检测待估值的线性调频信号功率谱P(x)；

步骤2、提取所述线性调频信号功率谱P(x)的形状特征CO(P(x))；

其中，通过以下过程获得所述形状特征：

a、选取x轴上的扁平结构元素g₁，长度为L₁，结构元素定义域

$D\left[g_1\right]=[x-\frac{L_1}{2},x+\frac{L_1}{2}],$

b、首先对P(x)进行形态学膨胀运算得到P₁(x)，表达式如下所示：

$P_1\left(x\right)=(P\oplus g_1)\left(x\right)=\max \{P\left(z\right):z\in D[g_1\left]\right\}$

其含义为以x点为中心，宽度等于扁平结构元素长度L₁的区域内 P(x)幅值的最大值；

再对P₁(x)进行形态学腐蚀运算得到P₂(x)，表达式如下所示：

P₂(x)＝(P₁Θg₁)(x)＝min{P₁(z):z∈D[g₁]}

其含义为以x点为中心，宽度等于扁平结构元素长度L₁的区域内 P₁(x)幅值的最小值；

最后对P₂(x)重复进行一次腐蚀运算和一次膨胀运算，表达式如下所 示：

$\mathrm{CO}\left(P\left(x\right)\right)=(\left(P_2\Theta g_1\right)\oplus g_1)\left(x\right);$

步骤3、计算所述形状特征上升沿的斜率β，并根据调频斜率k和所 述上升沿斜率β的关系计算所述线性调频信号的调频斜率k，其具体包 括：

a、设定线性调频信号调频斜率k变化区间为[50，2000]，步长为50；

b、计算不同k值对应的所述线性调频信号功率谱形状特征，搜索所 述线性调频信号功率谱形状特征上升沿峰值对应的频率f₁和半峰值对应 的频率f_0.5，根据定义计算出每个所述线性调频信号功率谱形 状特征对应的上升沿的斜率β，则获得一一对应的k-β数据点；

c、选择指数函数k(β)＝a(β^b)+c对所述k-β数据点进行数据拟合；

d、求得系数为：a＝1.4176，b＝-1.959，c＝15.4941；得到调频斜率k 与上升沿斜率β的关系如下式：

k(β)＝1.4176(β^-1.959)+15.4941

根据调频斜率k和所述上升沿斜率β的关系计算所述线性调频信号 的调频斜率k。

步骤4、根据所述形状特征上升沿的峰值频率f₁与起始频率f_s的频 率差Δf和所述调频斜率k的关系计算所述频率差Δf，其具体包括：

a、设定线性调频信号调频斜率变化区间为[50，2000]，步长为50；

b、计算不同k值对应的线性调频信号功率谱形状特征和每个所述线 性调频信号功率谱形状特征对应的上升沿频率差Δf，则获得一一对应的 Δf-k数据点；

c、选择指数函数Δf(k)＝a(k^b)+c对所述Δf-k数据点进行数据拟合；

d、求得系数为：a＝1.007,b＝0.4803，c＝-1.0466；得到调频斜率k与 频率差Δf的关系如下式：

Δf(k)＝1.007(k^0.4803)-1.0466

根据所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率差Δf 和调频斜率k的关系计算所述频率差Δf；

步骤5、计算所述线性调频信号的起始频率f_s和终止频率f_e，并计 算带宽B，其具体包括：

a、搜索所述形状特征上升沿和下降沿的峰值处频率f₁和f₁′；

b、根据所述起始频率f_s和频率差Δf的关系Δf＝f₁-f_s计算起始频率：

f_s＝f₁-Δf

c、根据所述终止频率f_e和频率差Δf的关系Δf＝f_e-f₁′计算终止频率：

f_e＝f₁′+Δf

d、根据公式B＝|f_e-f_s|，则所述线性调频信号的带宽：

B＝f_e-f_s；

步骤6、根据脉宽τ和所述带宽B的关系，计算所述脉宽τ。

其中，在步骤6中：

利用公式B＝kτ确定所述线性调频信号的脉宽：

$\tau =\frac{B}{k}.$

本发明的有益效果是：本发明公开了一种基于功率谱的线性调频信号 参数估值方法，该方法包括：检测待估值的线性调频信号功率谱P(x)； 提取所述线性调频信号功率谱P(x)的形状特征CO(P(x))；计算所述形状 特征上升沿的斜率β，并根据调频斜率k和所述上升沿斜率β的关系计算 所述线性调频信号的调频斜率k；根据所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率差Δf和所述调频斜率k的关系计算所述频率差Δf； 计算所述线性调频信号的起始频率f_s和终止频率f_e，并计算带宽B；根 据脉宽τ和所述带宽B的关系，计算所述脉宽τ。本发明方法的计算量主 要为一次FFT变换的复数乘法量，克服了目前常用的时频分析方法计算 量大的缺陷，在低信噪比环境下，实现未知先验信息条件的线性调频信 号参数的准确估值。

附图说明

图1本发明公开的一种基于功率谱的线性调频信号参数估值方法流 程图。

图2本发明公开的确定线性调频信号调频斜率k和形状特征上升沿 斜率β关系流程图。

图3线性调频信号调频斜率k与形状特征上升沿斜率β的数据点和 关系曲线。

图4不同调频斜率k下的线性调频信号功率谱形状特征。

其中，111为k＝100MHz/μs时的功率谱形状特征；112为 k＝200MHz/μs时的功率谱形状特征；113为k＝1000MHz/μs时的功 率谱形状特征；121为功率谱的起始频率f_s；122为功率谱形状特征上升 沿峰值对应的频率f₁；123为功率谱形状特征上升沿半峰值对应的频率 f_0.5；124为所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率差 Δf；125为所述功率谱形状特征上升沿的斜率β。

图5脉宽τ＝1,10,100μs下的线性调频信号功率谱形状特征。

其中，221为功率谱的起始频率f_s；222为功率谱形状特征上升沿峰 值对应的频率f₁；223为功率谱形状特征上升沿半峰值对应的频率f_0.5； 224为所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率差Δf；225 为所述功率谱形状特征上升沿的斜率β。

图6本发明公开的确定功率谱形状特征上升沿的频率差Δf与调频斜 率k的关系流程图。

图7功率谱形状特征上升沿的频率差Δf与调频斜率k的数据点和关 系曲线。

图8本发明公开的计算线性调频信号功率谱形状特征所使用的扁平 结构元素长度L₁选取图。

其中，*表示L₁＝5时的含噪声线性调频信号功率谱形状特征CO(P₁)与 零噪声线性调频信号功率谱P₀的相似度，○表示L₁＝11时的所述CO(P₁)与 所述P₀的相似度，△表示L₁＝15时的所述CO(P₁)与所述P₀的相似度。

图9本发明公开的线性调频信号功率谱及其形状特征。

其中，图（上）表示理想情况下线性调频信号功率谱，图（中）表 示信噪比为0dB时线性调频信号功率谱，图（下）表示0dB时线性调频 信号功率谱形状特征。

图10本发明公开的含噪声线性调频信号功率谱P₁和含噪声线性调 频信号功率谱形状特征CO(P₁)分别与零噪声线性调频信号功率谱P₀的相 似度分布图。

其中，○表示含噪声线性调频信号功率谱P₁与零噪声线性调频信号功 率谱P₀的相似度分布，+表示含噪声线性调频信号功率谱形状特征CO(P₁) 与零噪声线性调频信号功率谱P₀的相似度分布。

图11本发明实施例中调频斜率估计方法的均方根误差图。

其中，*代表脉宽τ＝1μs时，○代表脉宽τ＝100μs时。

图12本发明实施例中起始频率估计方法的均方根误差图。

其中，*代表脉宽τ＝1μs时，○代表脉宽τ＝100μs时。

图13本发明实施例中带宽估计方法的均方根误差图。

其中，*代表脉宽τ＝1μs时，○代表脉宽τ＝100μs时。

图14本发明实施例中脉宽估计方法的均方根误差图。

其中，*代表脉宽τ＝1μs时，○代表脉宽τ＝100μs时。

具体实施方式

以下将参照图1-14对本发明的实施方式进行说明。

如图1所示，本发明实施例进行基于功率谱的线性调频信号参数估 值方法包括下列步骤：

步骤1、检测待估值的线性调频信号功率谱P(x)；

步骤2、提取所述线性调频信号功率谱P(x)的形状特征CO(P(x))；

其中，步骤2中，通过以下过程获得所述形状特征：

a、选取x轴上的扁平结构元素g₁，长度为L₁，结构元素定义域

$D\left[g_1\right]=[x-\frac{L_1}{2},x+\frac{L_1}{2}],$

b、首先对P(x)进行形态学膨胀运算得到P₁(x)，表达式如下所示：

$P_1\left(x\right)=(P\oplus g_1)\left(x\right)=\max \{P\left(z\right):z\in D[g_1\left]\right\}$

其含义为以x点为中心，宽度等于扁平结构元素长度L₁的区域内 P(x)幅值的最大值；

再对P₁(x)进行形态学腐蚀运算得到P₂(x)，表达式如下所示：

P₂(x)＝(P₁Θg₁)(x)＝min{P₁(z):z∈D[g₁]}

其含义为以x点为中心，宽度等于扁平结构元素长度L₁的区域内 P₁(x)幅值的最小值；

最后对P₂(x)重复进行一次腐蚀运算和一次膨胀运算，表达式如下所 示：

$\mathrm{CO}\left(P\left(x\right)\right)=(\left(P_2\Theta g_1\right)\oplus g_1)\left(x\right);$

完成所述功率谱进行形态学闭-开运算。

需要说明的是，扁平结构元素长度L₁选取的原则如下：

扁平结构元素长度L₁的选取仅与采样点数有关，与信号参数无关， 而扁平结构元素长度的选取会影响最终得到的功率谱形状特征与噪声为 零时功率谱间的相似度以及计算量，因此扁平结构元素长度选取原则是 保证较高相似度的前提下选取较短的扁平结构元素长度。

需要说明的是，功率谱形状特征与噪声为零时功率谱间相似度的定 义如下：

设定含噪声功率谱形状特征序列为{CO(P_1n)}，噪声为零时功率谱序 列为{P_0n}，其中n＝1,2…N，和分别是序列{CO(P_1n)}和{P_0n}的 平均值，定义含噪声功率谱形状特征和噪声为零时功率谱的相似度r为：

$r=\left|\frac{\underset{n=1}{\overset{N}{\Sigma }}(\mathrm{CO}\left(P_{1n}\right)-\overline{\mathrm{CO}\left(P_1\right)})(P_{0n}-{\bar{P}}_0)}{\sqrt{\underset{n=1}{\overset{N}{\Sigma }}{(\mathrm{CO}\left(P_{1n}\right)-\overline{\mathrm{CO}\left(P_1\right)})}^2}\sqrt{\underset{n=1}{\overset{N}{\Sigma }}{(P_{0n}-{\bar{P}}_0)}^2}}\right|.$

步骤3、计算所述形状特征上升沿的斜率β，并根据调频斜率k和所 述上升沿斜率β的关系计算所述线性调频信号的调频斜率k；

其中，如图2所示，步骤3中，所述形状特征CO(P(x))上升沿的斜 率β和调频斜率k的关系的方法包括下列步骤：

a、设定线性调频信号调频斜率k变化区间为[50，2000]，步长为50；

b、计算不同k值对应的所述线性调频信号功率谱形状特征，搜索所 述线性调频信号功率谱形状特征上升沿峰值对应的频率f₁和半峰值对应 的频率f_0.5，根据定义计算出每个所述线性调频信号功率谱形 状特征对应的上升沿的斜率β，则获得一一对应的k-β数据点；

c、选择指数函数k(β)＝a(β^b)+c对所述k-β数据点进行数据拟合， 拟合数据点和拟合曲线如图3所示；

d、求得系数为：a＝1.4176，b＝-1.959,c＝15.4941；得到调频斜率k 与上升沿斜率β的关系如下式：

k（β)＝1.4176（β^-1.959)+15.4941

根据调频斜率k和所述上升沿斜率β的关系计算所述线性调频信号 的调频斜率k。

需要说明的是：

本发明选择指数函数k(β)＝a(β^b)+c对所述k-β数据点进行数据拟 合，拟合度为99%。

需要说明的是：

如图4和图5所示，当脉宽τ在[1,100]us的范围内变化时，所述线 性调频信号功率谱形状特征的上升沿几乎重合，因此可近似认为所述上 升沿斜率β不随所述脉宽τ的变化而变化；而所述上升沿斜率β会随着所 述调频斜率k的变化而变化，所以本发明近似认为上升沿斜率β的大小仅 与调频斜率k有关，与其他参数无关。

步骤4、根据所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率 差Δf和所述调频斜率k的关系计算所述频率差Δf；

其中，如图6所示，步骤4中，计算所述形状特征上升沿的峰值频 率f₁与起始频率f_s的频率差Δf与和所述线性调频信号调频斜率k的关系 的方法包括如下步骤：

a、设定线性调频信号调频斜率变化区间为[50，2000]，步长为50；

b、计算不同k值对应的线性调频信号功率谱形状特征和每个所述线 性调频信号功率谱形状特征对应的上升沿频率差Δf，则获得一一对应的 Δf-k数据点；

c、选择指数函数Δf(k)＝a(k^b)+c对所述Δf-k数据点进行数据拟合， 拟合数据点和拟合曲线如图7所示；

d、求得系数为：a＝1.007,b＝0.4803，c＝-1.0466；得到调频斜率k与 频率差Δf的关系如下式：

Δf(k)＝1.007(k^0.4803)-1.0466

根据所述形状特征上升沿的峰值频率f₁与起始频率f_s的频率差Δf 和所述调频斜率k的关系计算所述频率差Δf。

需要说明的是：

本发明选择指数函数Δf(k)＝a(k^b)+c对所述Δf-k数据点进行数据 拟合，拟合度为99%。

需要说明的是：

如图4和图5所示，当脉宽τ在[1,100]us的范围内变化时，所述线 性调频信号功率谱形状特征的上升沿几乎重合，因此可近似认为所述频 率差Δf不随所述脉宽τ的变化而变化；而所述频率差Δf会随着所述调频 斜率k的变化而变化，所以本发明近似认为频率差Δf的大小仅与调频斜 率k有关，与其他参数无关。

步骤5、计算所述线性调频信号的起始频率f_s和终止频率f_e，并计 算带宽B；

其中，步骤5中，计算所述起始频率f_s、终止频率f_e和带宽B的具 体方法包括以下步骤：

a、搜索所述形状特征上升沿和下降沿的峰值处频率f₁和f₁′；

b、根据所述起始频率f_s和频率差Δf的关系Δf＝f₁-f_s计算起始频率：

f_s＝f₁-Δf

c、根据所述终止频率f_e和频率差Δf的关系Δf＝f_e-f₁′计算终止频率：

f_e＝f₁′+Δf

d、根据公式B＝|f_e-f_s|，则所述线性调频信号的带宽：

B＝f_e-f_s

步骤6、根据脉宽τ和所述带宽B的关系，计算所述脉宽τ。

其中，步骤6中，计算所述脉宽τ的具体方法如下：

利用公式B＝kτ确定所述线性调频信号的脉宽：

$\tau =\frac{B}{k}$

本发明以线性调频信号为实施例对本发明进行详细说明；本发明具 体实施例用缩写LFM表示线性调频信号，用缩写SNR表示信噪比。仿 真参数设置：脉宽τ＝1μs，起始频率为500MHz，调频斜率k＝200MHz/μs。

本发明实施例可以应用到军用电子对抗以及民用频谱监督领域，在 低信噪比环境下准确地估计出信号参数，可以反映雷达的功能和用途， 提高声纳跟踪系统的定位精度与改善跟踪效果具有重要意义。

一、检测待估值的LFM信号功率谱P(x）；

二、提取所述线性调频信号功率谱P(x)的形状特征CO(P(x))；

选取合适的扁平结构元素长度L₁，对所述线性调频信号功率谱进行 形态学的闭-开运算并进行幅值归一化，得到所述线性调频信号功率谱 形状特征CO(P(x))。

1、扁平结构元素g₁长度L₁的选取

设定扁平结构元素长度L₁＝[5,11,15]，分别计算L₁＝[5,11,15]时所述含噪 声LFM信号功率谱形状特征CO(P₁)与噪声为零时功率谱P₀的相似度的分 布，如图8所示，当信噪比SNR＞0时，在L₁取不同值情况下相似度差异 很小；当信噪比SNR＜0时，在L₁取不同值情况下相似度有明显差异，但 随L₁值增大差异减小。因此，按照选取原则选取L₁，则确定L₁＝11。

2、计算LFM信号功率谱形状特征

利用形态学的闭-开运算计算所述LFM信号功率谱形状特征，信号 功率谱及功率谱形状特征如图9所示，图9（上）为噪声为零时信号功率 谱，图9中为信噪比SNR＝0dB时的信号功率谱，图9（下）为图9（中） 的形状特征。由图9可知，信噪比SNR＝0dB时LFM信号功率谱形状特 征与噪声为零时信号功率谱的形状基本不变。因此，本发明专利将LFM 信号功率谱形状特征近似为信号功率谱进行参数估计。

3、所述功率谱形状特征的抗噪性能分析

在信号参数不变的情况下，信噪比变化范围为-10dB~10dB时，利用 所述相似度公式分别计算含噪声LFM信号功率谱P₁与噪声为零时功率谱 P₀的相似度，以及含噪声LFM信号功率谱形状特征CO(P₁)与噪声为零时 功率谱P₀的相似度，分布如图10所示。在整个信噪比变化范围内，CO(P₁) 与P₀间的相似度均大于P₁与P₀间的相似度，且随着信噪比的降低相似度的 差异增大。当信噪比SNR＝0dB时，P₁与P₀间的相似度为0.9273，CO(P₁)与 P₀间的相似度为0.9843，相似度提高了6.15%；当信噪比SNR＝-10dB时， P₁与P₀间的相似度为0.379，CO(P₁)与P₀间的相似度为0.7177，相似度提 高了89.37%，证明在低信噪比情况下LFM信号功率谱形状特征的抗噪 性能更好。

三、计算所述LFM信号的调频斜率k

根据拟合公式k(β)＝1.4176(β^-1.959)+15.4941计算调频斜率k，在SNR 为[-5dB,15dB]时进行100次Monte Carlo实验，由图11可知，计算得到 的调频斜率的估计误差不受脉宽τ变化的影响，而且在SNR=-5dB时调频 斜率估计误差均小于k＝1.2MHz/μs。

四、根据拟合公式Δf(k)＝1.007(k^0.4803)-1.0466计算所述形状特征上升 沿的频率差Δf

五、计算起始频率f_s、终止频率f_e和带宽B

对SNR在区间为[-5dB,15dB]时的LFM信号进行100次Monte Carlo 试验，计算所述LFM信号起始频率和带宽的均方根误差（RMSE），结果 如图12和图13所示。

由图12可知，本发明方法在SNR为[-5dB,15dB]时的起始频率的估 计误差不受脉宽τ变化的影响，而且在SNR=-5dB时起始频率估计误差均 小于3.1MHz。

由图13可知，本发明方法在SNR为[-5dB,15dB]时的带宽的估计误 差不受脉宽τ变化的影响，而且在SNR=-5dB时带宽估计误差均小于 2.3MHz。

六、计算所述LFM信号的脉宽τ

对SNR在区间为[-5dB,15dB]时的LFM信号进行100次Monte Carlo 试验，计算脉宽估计的均方根误差（RMSE），结果如图14所示；

由图14可知，本发明方法在SNR为[-5dB,15dB]时的脉宽的估计误 差不受脉宽τ变化的影响，而且在SNR=-5dB时脉宽估计误差均小于 0.025us。