基于傅里叶梅林变换的宽带非线性调频信号时差/尺度差估计方法与流程

1.本发明属于信号处理领域，适用于解决双星无源定位系统中的时差/尺度差估计难题。


背景技术：

2.电子战是信息化战争的重要组成部分，对目标的精确定位是电子战的一项重要内容。目标定位系统通常可分为有源定位系统和无源定位系统，典型的有源定位系统为雷达，由于需要对外辐射电磁信号以获取目标回波，战场上容易被敌方侦察而遭受反辐射武器打击。相较而言，无源定位系统自身不辐射电磁信号，隐蔽性好，典型体制为双星无源定位。宽带条件下，双星无源定位系统根据两路接收信号估计时差(tdoa)和尺度差(sdoa)，再利用卫星、目标和地球的几何关系对辐射源进行精确定位。为估计时差/尺度差，传统方法是计算两路接收信号宽带互模糊函数(wbcaf)，该方法能够同时估计2个参数，且适用多种源信号形式，存在的问题是计算量大。wbcaf计算过程包含两步，首先计算一路信号的尺度版本，然后将其与另一路信号进行相关处理，由于第2步可利用fft、ifft在频域快速实现，因此制约wbcaf快速实现的难点为第1步，针对该问题，文献[王钊,王鹏毅,苏卫民,等.宽带无源时差尺度差联合估计算法研究[j].无线电工程,2016,46(11):34
‑
37]提出多速率采样重构滤波方法，但当尺度因子为多位小数时，该方法无法获得精确的尺度版本。部分研究考虑仅在变换域计算wbcaf，以规避尺度版本计算问题，其中，文献[niu x x,ching p c.wavelet based approach for joint time delay and doppler stretch measurements[j].ieee transaction on aerospace and electronic systems,1999,35(3):1111
‑
1119]利用小波变换内积定理，提出基于交互小波变换的wbcaf实现方法，但母小波的选取对计算结果影响较大。文献[王军,李亚安.基于梅林变换的宽带水声信号处理研究[j].兵工学报,2007,28(1):87
‑
90]提出基于梅林变换的wbcaf实现方法，但梅林变换的数值计算方法有多种，所用方法计算精度不高。另外，如果源信号形式已知，时差/尺度差估计问题会变得相对容易，如针对线性调频源信号，文献[郭付阳,张子敬,杨林森.基于尺度变换的宽带线性调频信号时差/尺度差估计算法[j].电波科学学报,2017,32(4):441
‑
448]和文献[黄清顺,邓兵.基于分数阶傅里叶四阶中心矩的尺度差/时差估计[j].信号处理,2021,37(1):111
‑
119]利用分数阶傅里叶变换(frft)测量两路接收信号调频斜率，根据两路信号调频斜率比例关系估计尺度差，再对第1路信号“拉伸”或者“压缩”，与第2路信号相关处理估计时差。该改进方法通过两步级联处理，降低了运算整体复杂度，低信噪比条件下时差/尺度差估计均方误差逼近克拉美罗下界，遗憾的是当源信号为非线性调频信号时，该方法无法适用。
[0003]
针对上述问题，本发明利用两路信号频谱包络是彼此尺度版本的特点，提出基于傅里叶梅林变换的宽带非线性调频信号时差/尺度差估计方法。核心思想是利用傅里叶变换“时移不变性”、梅林尺度变换“尺度不变性”对尺度差/时差进行级联估计。仿真结果表明，当snr大于5db，本发明能够准确估计二次调频(qfm)、指数调频(efm)、双曲调频(hfm)、
正弦调频(sfm)4种宽带非线性调频信号时差和尺度差，时差估计平均相对误差小于1％，尺度差估计平均相对误差小于5％。


技术实现要素：

[0004]
本发明的目的是提出基于傅里叶梅林变换的宽带非线性调频信号时差/尺度差估计方法，以解决双星无源定位系统中的时差/尺度差估计难题。为阐述本发明，首先，简要介绍传统时差/尺度差估计方法与改进方法；其次，引入2个基本概念，时间尺度(ts)和尺度估计(se)；再次，对本发明的基本原理进行介绍；然后，对比传统基于宽带互模糊函数的时差/尺度差估计方法，分析本发明计算复杂度；最后，介绍本发明方法步骤。
[0005]
(1)传统时差/尺度差估计方法与改进方法
[0006]
宽带条件下，利用双星无源定位系统对地面一个固定辐射源进行定位，以第1路信号为参考基准，则2个独立接收机接收到的信号可表示为：
[0007][0008]
式中：x1(t)、x2(t)分别为第1路和第2路接收信号，t≥0为时间自变量，a1、a2为源信号s(t)到达2个接收机的路径增益，δt为两路信号时差(tdoa)，δα为尺度差(sdoa)，与卫星相对辐射源径向速度和光速有关，n1(t)、n2(t)为接收机基底复高斯白噪声。传统时差/尺度差估计方法是利用宽带互模糊函数(wbcaf)实现。连续信号z(t)与x(t)的wbcaf定义为：
[0009][0010]
式中：γ为wbcaf尺度因子，τ为时移，*为共轭符号。令z(t)＝x2(t)、x(t)＝x1(t)，易知当γ＝1 δα、τ＝δt时，|x
zx
(γ,τ)|取得最大值，峰值搜索可得到时差和尺度差(需要减1)。wbcaf可表示成卷积形式频域快速实现，即：
[0011][0012]
式中：为卷积符号，f[
·
]、f
‑1[
·
]为傅里叶变换与逆变换表示符号。为计算两路接收信号wbcaf，需设定尺度因子搜索区间，然后根据不同的尺度因子对第1路信号进行“拉伸”或“压缩”，搜索区间的扩大会导致运算量成倍增加。当源信号形式已知时，时差/尺度差估计问题会变得相对容易，如源信号为单频信号，则尺度上的差异反应为接收信号频率上的不同，如果源信号为lfm信号则反应为调频斜率上的不同，通过估计两路信号频率或者调频斜率可得sdoa，然后根据估计的sdoa对x1(t)进行“拉伸”或者“压缩”，再与x2(t)时域相关处理可得tdoa，该改进方法虽然能够显著降低运算量，但当源信号为nlfm信号时无法适用。雷达、通信系统中常用的nlfm信号包括二次调频(qfm)、指数调频(efm)、双曲调频(hfm)和正弦调频(sfm)等，复数形式依次为：
[0013][0014]
[0015][0016][0017]
式中：t为时宽，f1、f2分别为4种信号的起始频率和截止频率，t0为hfm延时时间。相位求导可得qfm、efm、hfm、sfm瞬时频率分别为：
[0018][0019][0020][0021][0022]
(2)时间尺度与尺度估计
[0023]
所谓时间尺度(ts)，即连续信号x(t)到y(t)的映射过程，可表示为：
[0024][0025]
式中：ts
α
[
·
]为ts表示符号，α∈r

为尺度因子自变量，是为保持ts前后信号能量相同，即：
[0026][0027]
当0＜α＜1时，ts后信号时域扩张，当α＞1时，ts后信号时域收缩。根据傅里叶变换尺度性质可知：
[0028][0029]
式中：y(f)、x(f)分别为y(t)、x(t)的傅里叶变换，f为频率自变量。可以看出，对信号进行ts还会导致带宽成比例的减小(0＜α＜1)或增大(α＞1)。设连续信号α1为z(t)尺度因子。所谓尺度估计(se)，即z(t)、x(t)已知时对α1的估计，可表示为：
[0030][0031]
式中：为α1的估计值，φ
zx
(α)为z(t)与x(t)的尺度互相关函数(scf)，即：
[0032][0033]
式中：
◇
为尺度互相关符号。当τ＝0时，wbcaf与scf近似，区别在于不同的积分区间。尺度变换(st)是梅林变换(mt)复自变量实部取0.5的特例，表示式为：
[0034][0035]
式中：s{
·
}为st表示符号，d
f
(c)为信号f(t)的st，c为尺度自变量。逆尺度变换(ist)为：
[0036][0037]
式中：s
‑1{
·
}为ist表示符号。令t＝e
u
，得到：
[0038][0039]
可以看出，st实际为的傅里叶变换，可对f(t)指数采样，再利用fft快速实现，即快速尺度变换(fst)，一般称其为指数采样型fst。文献[sena d a,rocchesso d.a fast mellin and scale transform[j].eurasip journal on advances in signal processing,2007:089170]通过分析各实现环节，得出指数采样型fst计算复杂度为o[(n ln n)log2(n ln n)]，n为原始信号等间隔采样点数，总共需要的复乘次数为n ln n 0.5n ln nlog2(n ln n)。易知：
[0040][0041]
可以看出，ist同样可利用ifft快速实现，得到快速逆尺度变换(ifst)，计算过程与fst相反，先计算d
f
(c)的逆傅里叶变换，再进行对数采样。为解决ts快速实现问题，引用st尺度不变性，即：
[0042][0043]
式中：d
x
(c)、d
z
(c)分别为x(t)、z(t)的st，两者包络相同、相位不同。计算上式ist，容易得到：
[0044][0045]
利用上式可得x(t)在某一尺度因子下的ts后信号。se同样利用st快速实现，具体为：
[0046][0047]
(3)本发明基本原理
[0048]
为介绍本发明，暂不考虑接收机噪声，容易得到：
[0049][0050]
式中：x1(f)、x2(f)分别为x1(t)、x2(t)的傅里叶变换，f≥0为频率自变量。易知：
[0051][0052]
可以看出，|x2(f)|实际为|x1(f)|的尺度版本，尺度因子为1/(1 δα)，与源信号形式无关。当|x1(f)|、|x2(f)|已知时，对1/(1 δα)的估计属尺度估计概念，计算两路信号频谱包络尺度互相关函数，可估计1/(1 δα)，进而得到尺度差δα，然后对x1(t)进行时间尺度(尺度因子取1 δα)，得到时间尺度后的第1路信号计算x2(t)与的时域互相关函数τ为时延自变量，取包络得到峰值搜索得到对应的时延τ0即为估计的时差δt。时域互相关函数计算公式为：
[0053][0054]
其中，ψ
zx
(τ)为连续信号z(t)与x(t)的时域互相关函数，
□
为时域互相关符号。ψ
zx
(τ)利用fft和ifft实现，即：
[0055]
ψ
zx
(τ)＝f
‑1{f[z(t)]f
*
[x(γt)]；τ}
[0056]
(4)计算复杂度分析
[0057]
首先，分析时间尺度(ts)计算复杂度。ts实现过程包含“fst 相位修正 ifst”三个环节，计算复杂度为o[(n ln n)log2(n ln n)]，总共需要进行3n ln n n ln nlog2(n ln n)次复数乘法。其次，分析尺度估计(se)计算复杂度。se实现的关键在于计算尺度互相关函数，具体需要进行2次fst、1次复数乘积和1次ifst，计算复杂度为o[(n ln n)log2(n ln n)]，总共需要进行4n ln n 1.5n ln nlog2(n ln n)次复数乘法。再次，分析传统基于宽带互模糊函数的时差/尺度差估计方法计算复杂度。对于1个尺度因子，计算wbcaf需要进行1次ts、2次fft、1次ifft和1次点乘，总共需要进行mn 1.5mnlog2n 3mn ln n mnln nlog2(n ln n)次复乘运算，计算复杂度为o[(mn ln n)log2(n ln n)]，m为尺度因子取值个数。最后，分析本发明计算复杂度。为估计尺度差需要进行2次fft、1次se，而估计时差需要进行1次ts和1次时域相关处理(2次fft、1次ifft和1次点乘)，总共需要进行n 2.5nlog2n 7n ln n 2.5n
·
lnnlog2(n ln n)次复乘运算，计算复杂度为o[(n ln n)log2(n ln n)]。由于传统方法需要进行1次二维搜索(本发明需要进行2次一维搜索)，实际运算量会更大。
[0058]
(5)本发明方法步骤
[0059]
本发明主要包括以下步骤：
[0060]
步骤(一)计算两路信号傅里叶变换，取包络得到两路信号频谱；
[0061]
步骤(二)计算第2路信号频谱与第1路信号频谱的尺度互相关函数，峰值搜索估计尺度差；
[0062]
步骤(三)根据估计的尺度差对第1路信号进行时间尺度；
[0063]
步骤(四)计算第2路信号与时间尺度后的第1路信号时域互相关函数，峰值搜索估计时差。
[0064]
上述步骤具体为：
[0065]
步骤(一)计算第1路信号x1(t)和第2路信号x2(t)的傅里叶变换，取包络分别得到|x1(f)|、|x2(f)|，其中，t≥0为时间自变量，f≥0为频率自变量，x1(f)、x2(f)分别为x1(t)、x2(t)的傅里叶变换；
[0066]
步骤(二)计算|x2(f)|与|x1(f)|的尺度互相关函数α∈r

为尺度因子自变量，取包络得到峰值搜索得到对应的尺度因子α0，估计得到尺度差δα为1/α0‑
1；
[0067]
步骤(三)令α＝1 δα，对第1路信号x1(t)进行时间尺度，得到经时间尺度操作后的第1路信号为ts
α
[
·
]为时间尺度表示符号；
[0068]
步骤(四)计算x2(t)与的时域互相关函数τ为时延自变量，取包络得到峰值搜索得到对应的时延τ0即为估计的时差δt。
[0069]
本发明的有益效果说明：
[0070]
相比传统基于wbcaf的时差/尺度差估计方法，本发明计算复杂度更低，相比基于frft的时差/尺度差估计的改进方法，本发明对宽带lfm、nlfm等信号形式均适用。
附图说明
[0071]
附图1是本发明的方法步骤流程图；
[0072]
附图2是复乘次数随采样点数变化曲线；
[0073]
附图3是nlfm信号时域波形；
[0074]
附图4是nlfm信号时频分布；
[0075]
附图5是接收信号时域波形和频谱；
[0076]
附图6是接收信号频谱的尺度互相关函数；
[0077]
附图7是时间尺度后的接收信号时域波形和频谱；
[0078]
附图8是第2路信号与时间尺度后的第1路信号时域互相关函数；
[0079]
附图9是尺度差估计平均相对误差随信噪比变化曲线；
[0080]
附图10是时差估计平均相对误差随信噪比变化曲线。
[0081]
具体实施方法
[0082]
下面结合附图对本发明进行详细描述。附图1为本发明的方法步骤流程图。本发明具体步骤如下所示：
[0083]
步骤(1)计算两路信号傅里叶变换，取包络得到两路信号频谱；
[0084]
步骤(2)计算第2路信号频谱与第1路信号频谱的尺度互相关函数，峰值搜索估计尺度差；
[0085]
步骤(3)根据估计的尺度差对第1路信号进行时间尺度；
[0086]
步骤(4)计算第2路信号与时间尺度后的第1路信号时域互相关函数，峰值搜索估计时差。
[0087]
设传统方法wbcaf尺度因子取值个数分别为10、100、1000，附图2为本发明与传统方法复乘次数随采样点数变化曲线，传统方法计算量随尺度因子取值个数成倍增加，而本发明计算量明显低于传统方法。
[0088]
实施条件：在以下参数条件下进行仿真实验：
[0089]
表1辐射源和接收信号
[0090][0091]
设辐射源发射qfm、efm、hfm、sfm 4种nlfm信号，时宽2.5μs，带宽250mhz，起始频率和截止频率分别为10mhz和260mhz，hfm延时时间为1个时域采样间隔，双星无源定位系统采样频率为1ghz，接收的两路信号时差/尺度差分别为3.75μs和0.1。4种nlfm信号时域波形如图3所示，利用同步压缩变换得到nlfm信号时频分布如图4所示。设源信号为sfm信号，接收信号时长为10μs，按照图1流程，对本发明可行性进行验证。两路接收信号时域波形和频谱
如图5所示。计算两路信号频谱的尺度互相关函数估计尺度差，结果如图6所示。根据估计的尺度差，对第1路信号进行时间尺度，时间尺度后的接收信号时域波形和频谱如图7所示。第2路信号与时间尺度后的第1路信号时域互相关函数如图8所示。分别选取lfm信号(起始频率、截至频率同上)、4种nlfm信号和零均值高斯白噪声作为源信号，以平均相对误差(mre)为指标，评估传统方法和本发明抗噪效能。snr取值
‑
10～10db，间隔1db，运行蒙特卡洛仿真500次，得到尺度差估计平均相对误差随信噪比变化曲线如图9所示，图10为时差估计平均相对误差随信噪比变化曲线。通过附图5可以看出，由于尺度差的存在，第2路信号中的sfm信号相比第1路信号中的sfm信号存在时域收缩和频谱扩张现象；通过附图6可以看出，当对频谱取包络，估计得到的尺度因子为0.9092，尺度差为0.0999，与仿真所设尺度差0.1一致，如果不对频谱取包络，估计结果完全错误；通过附图7可以看出，两路接收信号中的sfm信号时宽、带宽基本一致，不再存在收缩、扩展现象；通过附图8可以看出，时域互相关函数近似为辛格函数，峰值搜索估计得到时差为3.75μs，与仿真所设参数完全一致；通过附图9、附图10可以看出，传统方法与本发明对不同信号形式时差/尺度差估计效能不尽相同，总体上看，传统方法更适用低snr条件，当snr大于
‑
4db时，对lfm、qfm、efm、sfm 4种信号尺度差估计mre小于5％，时差估计mre小于1％，而本发明则要求snr大于5db，需要注意的是，高snr条件下本发明估计精度要高于传统方法，原因是因为计算wbcaf需人为设定尺度因子取值范围和间隔(仿真中设尺度因子取值范围为0.5～1.5，间隔0.01)，当取值范围一定时，降低间隔能够提高估计精度，但同时也需要更大的计算量。传统方法对hfm尺度差估计效果极不理想，不同snr条件下尺度差估计mre大于13％，估计误差明显高于本发明，原因是因为hfm具备“多普勒不变性”，其宽带模糊函数脊线近似为线性，脊线幅度较为接近，峰值搜索无法准确估计时差/尺度差，而本发明是对信号频谱包络进行尺度估计，不存在该问题。同时可以看出，传统方法与本发明对噪声信号尺度差估计mre小于5％的临界snr分别为2db和3db，时差估计mre小于1％的临界snr分别为1db和2db，估计效能相差不大。