张朝星,赵兴刚,韩晓军
(1.解放军32137部队,河北 张家口 075000;2.解放军66136部队,北京 100042;3.解放军31675部队,河北 张家口 075000)
时频分析[1]通过将一维时域信号转换为二维时频域图像,从而提取有效的图像特征来判断有无目标,是实现海杂波背景下目标检测[2]的一种有效方法。短时傅里叶变换(STFT)作为一种常用的时频分布[3],其基本思想是用一个确定宽度的窗函数沿时间轴逐段对信号进行傅里叶分析,得到信号的谱图。STFT具有许多优点,然而也有其缺点和不足,根据不确定性原理[4],其时间分辨率与频率分辨率的乘积恒定,不能同时提高,并且窗的宽度不能根据不同频率实现自适应调整,时频分辨力较低,导致使用STFT方法对海杂波背景下的弱小目标进行检测时效果较差。
曲线拟合选择一条合适的曲线来最佳地拟合观测数据,在欧式空间中,常见的曲线类型有二维直线、抛物线、椭圆等。在功率谱的信息几何理论[5]中,所有功率谱的集合可以看做一个黎曼流形[6],则对回波数据进行STFT处理后所得到的不同时间点的功率谱就对应了流形上的若干点。类比于欧式空间中的曲线拟合,本文提出了一种使用流形上的距离和测地线[7]对不同时间点的功率谱进行曲线拟合以提高其时频分辨力的方法。由于STFT是一种线性变换,我们选择欧式空间中的直线在流形上的推广(即测地线)对不同点的功率谱进行拟合。在具体的实现过程中,欧式空间使得给定数据与所选择曲线对应点的欧式距离总的平方和最小(最小二乘法)[8]。转换到流形上以后,实现方法不变,但欧式距离要替换为两功率谱间的距离。我们选择M-K运输问题[9]中的运输距离作为两功率谱之间的差异度量,得到了相应的黎曼度规,对STFT后的谱图实现了测地线拟合。拟合后的谱图对于单一点的功率谱进行了重构,使其更加平滑和准确。其次拟合所用测地线对不同时间点的功率谱实现了时域上的平滑,时间分辨率也得到提高,进而整体提高了谱图的时频分辨力,STFT的检测性能得到改善。
在上节中已提到,具体实现流形上的拟合时,要将最小二乘法中的欧式距离替换为2个不同时间点功率谱之间的距离,且该距离除了满足对称性和三角不等式等距离的一般性质外,还必须要满足弱连续性[10],以使得功率谱能够收敛,即随着2个功率谱的“靠近”或“远离”,该距离能相应地减小或增大。因此,我们选用了M-K运输问题中的运输距离来度量两功率谱之间的差异,对于2个归一化的密度函数f0和f1,其运输距离定义如下[9]:
(1)
式中:X=[-π,π];ψ(x)为一个f0和f1之间的保测度映射,即:|detψ(x)|f1(ψ(x))=f0(x)。
在一般情况下,ψ(x)可通过下式计算得到[9]:
F0(θ)=F1(ψ(θ))
(2)
(3)
当f1是f0的一个小的增量时,可由运输距离得到一个黎曼度规,并称之为运输度规[11]:
(4)
(5)
定义了流形上的度规,就可得到f0和f1之间由梯度流决定的测地线[12]:
Fτ((1-τ)θ+τψ(θ))=F0(θ)
(6)
式中:ψ(θ)可以通过式(2)计算;Fτ为fτ的积累函数。
有了2个功率谱之间的距离和测地线,我们就可以在流形上对功率谱进行曲线拟合。假设给定一个功率谱序列:G={gτi(θ):θ∈[-π,π],i=0,1,…,n},其中τi是一个指代时间的归一化序列,且τ0=0,τn=1。这些功率谱通常可以通过对观测数据进行短时傅里叶变换(STFT)得到,τi(i=0,1,…,n)可看做STFT中对应时间窗的中点。令d代表两功率谱之间的距离,则曲线拟合的过程就是寻找一条测地线fτ,τ∈[0,1],使得下式取得最小值:
(7)
利用运输距离及其相应的测地线来实现该过程,则整个测地线拟合的目标就是要确定一条测地线fτ,τ∈[0,1],使得下式取得最小值:
(8)
任意一条测地线都可以完全由两“点”确定,在这里就是指f0和f1。另外,它也可以根据式(2)和式(6)中的Ψ确定。f0,f1和Ψ是根据STFT处理后得到的功率谱G来确定的,根据式(3),有:
(9)
根据式(6),由ψ(θ)可得f0和f1之间测地线fτ的具体表达式为:
F0(θ)=F1(ψ(θ))=Fτi((1-τi)θ+τiψ(θ))
(10)
(11)
将式(11)代入式(9),则目标函数可写为:
(12)
(13)
(14)
所以Jg(f0,f1)可以通过以下的求和式进行近似:
(15)
所以,功率谱的测地线拟合就转化为了以式(15)为目标函数的最优化问题,其线性约束为:
(16)
图1 功率谱测地线拟合示意图
本节采用加拿大McMaster大学提供的IPIX雷达#320号数据[14]对本文所提方法的检测性能进行验证,该数据文件包含14个距离单元,每个单元由131 072个采样样本构成,其中第7个单元为主目标单元,第6、8、9单元为次目标单元,其余为纯杂波单元。
首先将第1个距离单元的采样数据作为杂波背景,通过STFT分析海杂波本身的谱特性。取观测数据样本长度为256,窗函数宽度为32,时间窗重叠区域宽度为16,脉冲重复频率fr=1 000 Hz,可得第1个距离单元的海杂波二维谱图,如图2所示。
图2 第1个距离单元的海杂波二维谱图
从图2中可以看出,海杂波谱是非对称的,由于海浪的移动,使得谱的主瓣偏离零频,杂波能量主要集中在[-200 Hz,0 Hz]这个区间内,下面将仿真目标信号加到第1个单元的纯海杂波数据中,设雷达在该距离单元观测的复包络信号为:
x(n)=s(n)+v(n),n=1,…,N
(17)
式中:s(n)=aej(2πfdn/fr),a为信号幅度,fd为目标多普勒频率。
令信噪比分别取0 dB和4 dB,目标多普勒频率分别取-100 Hz(主杂波谱区)和100 Hz(非主杂波谱区),分别给出相应的STFT变换后的谱图,如图3所示。
图3 海杂波加仿真信号后STFT三维谱图
从图3(a)和图3(c)中可以看出,当目标频率处于非主杂波区(fd=100 Hz)时,可以非常轻易地从STFT变换后的时频图像中判断目标的存在。但当目标频率处于主杂波区(fd=-100 Hz)时,如图3(b)和图3(d)所示,信噪比较高时,还可以从谱图中对目标进行辨别;但在低信噪比(SNR=0 dB)时,目标被杂波淹没,受谱图分辨力的限制,很难判断目标的存在。
然后就可以用Matlab中的最优化工具箱对该问题进行求解。对应图2中的4种情况,拟合得到的4条测地线的2个端点如图4所示。
图4 拟合得到的测地线端点图
图5 拟合处理后的谱图
从图4中单个时间点的功率谱已经可以初步判断各个情况下目标的存在,但当目标处于主杂波区时,还不够明显。根据式(11),有了2个端点就可以确定一条测地线,进而就得到了测地线上各个时刻的功率谱,如图5所示。
从图5中可以看出,当目标处于主杂波区(fd=-100 Hz)并且信噪比较低时(SNR=0 dB),与图3(d)相比,图5(d)中可以明显地看到一条凸出的“山脊”,其所在频率位置与目标多普勒频率相近,进而就可以较容易地判断目标的存在。
以上结果都是将仿真目标信号加到纯杂波数据上进行处理的,下面直接对实际有目标的第7个距离单元数据进行STFT变换,并进一步进行测地线拟合,给出处理结果如图6所示。
图6 第7个单元(主目标单元)STFT+测地线拟合处理后的谱图
IPIX实测数据是由雷达对漂浮于海面直径约1 m的金属球持续照射得到的,该金属球会随着海浪的波动进行漂移,但其运动速度与海浪运动速度相比相对较小,因而体现在功率谱上即其能量主要集中在零频附近,且其运动是无规律的,目标RCS和速度都会随着时间推移而发生变化,因而其功率谱也会随时间变化。从图5中可以非常明显地看到,在零频附近有一条凸出的脊线,而且这条脊线并不像加仿真目标信号时形状比较规则,而是在时域上有一定起伏和变化,这与实际目标状况是相符的,因而就可以判断目标的存在。
本文根据功率谱的信息几何理论,提出利用流形上的距离和测地线对STFT变换后的不同时间点的功率谱实现曲线拟合,以提高其时频分辨力;并通过推导,将该拟合过程转换为可用数值方法解决的具有线性约束的最优化问题,使其容易求解;最后在仿真实验中,分别利用IPIX实测纯海杂波数据加仿真目标信号和本身带目标的海杂波数据2种方式对所提方法的性能进行了测试,验证了本文方法的有效性。