南京航空航天大学学报  2017, Vol. 49 Issue (5): 699-706   PDF    
基于多普勒频率反投影的NCW-SAR动目标成像方法
汪玲, 岳怡然     
南京航空航天大学电子信息工程学院,南京,211106
摘要: 传统合成孔径雷达(Synthetic aperture radar, SAR)成像依赖大带宽发射信号,研究证明SAR也可利用窄带连续波形(Narrowband continuous wave, NCW)具有的高多普勒分辨率实现场景的高分辨成像。NCW-SAR成像的优势在于占有很少的频带资源、降低系统成本和规模,更适用于通常带宽较窄的利用外辐射源波形的无源SAR。文中首先指出NCW-SAR的应用前景,简要回顾窄带动目标成像的发展,并分析NCW-SAR成像技术的现状。然后给出通用的SAR回波模型,提出NCW-SAR动目标成像方法,实现静止场景中动目标的聚焦,获得目标的速度估计,并分析速度分辨率。实验验证了NCW-SAR动目标成像方法的有效性。
关键词: 信号与信息处理     雷达     成像     合成孔径雷达     连续波     动目标    
Doppler-Projection-Based NCW-SAR Moving Target Imaging Method
WANG Ling, YUE Yiran     
College of Electronic and Information Engineering, Nanjing University of Aeronautics & Astronautics, Nanjing, 211106, China
Abstract: Traditional synthetic aperture radar (SAR) imaging relies on wideba nd transformed waveforms. SAR using narrowband continuous waveforms (NCW) has be en proved to be also able to achieve high resolution imaging exploiting high Dop pler resolution. NCW-SAR has the following advantages: it occupies very small band resources, it promotes the development of low-cost and small SAR, it is much applicable to passive SAR using illuminators of opportunity, which usually trans mits waveforms with relatively narrow bandwidth. Firstly, the application pot ential of NCW-SAR is pointed out. The development of narrowband imaging is brie fly reviewed and the state of art NCW-SAR imaging technology is analyzed. Then, ageneral SAR received signal model is presented and afterwards, a moving target imaging method for NCW-SAR is proposed to focus the moving targets on the ground scene. The velocity estimation method as well as the velocity resolution analys is is also presented. The simulation results demonstrate the effectiveness of th e proposed NCW-SAR moving target imaging method.
Key words: signal and information processing     radar     imaging     synthetic aperture radar     continuous wave     moving target    

合成孔径雷达(Synthetic aperture radar, SAR)高分辨成像技术发展至今,已经进入比较成熟的阶段,尤其在工程应用方面,SAR成像理论得到了较深入的探讨。SAR成像依赖发射大带宽信号获得距离/延时高分辨,利用雷达平台和目标之间相对运动在方位维合成大孔径,通过相干积累所有回波获得方位向高分辨。本文探讨一种新的SAR成像体制下的动目标成像技术,该成像体制不需要宽带发射波形,不依赖宽带波形的距离分辨率,而是采样窄带连续波(Narrowband continuous wave, NCW)发射波形,利用连续波波形自身具有的高多普勒分辨率进行高分辨率SAR成像,简称这种SAR为窄带连续波SAR(NCW-SAR)。

NCW-SAR的应用优势和未来应用前景主要体现在以下方面:第一,节约频带资源,缓解频带拥挤。随着电子信息和通信技术的迅速发展,可用频带资源日益减少。NCW-SAR在获得高分辨的同时频带占用率低,可以节约频带资源和提高频谱利用率,在未来雷达发展中具有优势。第二,降低SAR系统复杂度,促进SAR系统小型化和低成本。宽带SAR一般工作在脉冲模式,雷达发射系统较为复杂。NCW-SAR采用连续波(Continuous wave, CW)发射波形,而且波形带宽窄,可简化雷达结构,降低雷达系统成本。第三,促进无源SAR的发展。近年来,广播台、电视台、移动通信基站、导航卫星等民用机会照射源的数目迅速增加,基于外部照射源的无源雷达可缓解频带拥挤问题。目前无源SA R研究依赖发射源波形的带宽,但是外部照射源信号基本都为带宽较窄的连续波形式,NCW -SAR成像为无源SAR提供了新的成像手段。

国外研究人员从20世纪80年代开始窄带微波成像的研究[1],而窄带成像思想最早也用于医学成像领域,目的是为减小医疗设备成本、开发便携式设备和降低诊断对人体的伤害,在本质上与微波领域的窄带成像方法互为借鉴。已有的窄带动目标成像研究可分为:(1)针对远场静止点目标和运动目标的频域图像重建;(2)基于匹配滤波处理(Matched filter processing, MFP)的时域图像重建;(3)基于压缩感知(Compressive sensing, CS)理论的图像重建。

针对远场静止稀疏点目标和运动目标的频域图像重建方法也称之为层析成像。Mensa采用单频连续波对远场稀疏目标进行成像[1],称之为相干多普勒层析成像。Sun等人在2010年对这种依靠单频连续波信号进行转台目标成像的方法进行了实验验证[2]。Mensa等人[1]给出的成像方法与远场近似条件和平面波前假设下的极坐标格式算法相似。Munson研究小组开展的对飞机目标的无源多基成像研究也可看作频域非均匀采样的图像重建算法[3]。基于MFP的时域图像重建方法大多数用于固定分布式孔径的窄带运动目标成像。分布式雷达的运动目标窄带成像,称之为动目标层析成像(Tomography of Moving Targets, TMT)。文献[4]中提出的分布式孔径动目标无源成像,验证了单频CW信号对动目标定位和速度估计的优良性能。文献[5]中针对单频CW SAR成像,给出基于MFP的时域图像重建方法,称之为多普勒成像。基于压缩感知(Compressi ve sensing, CS)理论的图像重建方法的优点在于可以使用较少的观测数据进行成像[6]。窄带雷达带宽较小,采样数据在波数域呈稀疏分布,因此可以采用基于CS理论的图像重建方法进行窄带雷达成像。

国内近些年也开展了窄带雷达成像研究[7-9],包括基于窄带外辐射源单/多基配置下的动目标二维和三维ISAR成像的研究[7, 10],提出了基于PFA的波数域/频域重建方法[11]和基于匹配思想的时域重建方法[7],以及利用CS方法解决窄带信号距离分辨率不足[10]

窄带雷达使用的CW波形具有良好的多普勒分辨率,已有研究人员注意到这个良好的性质,开展了基于多普勒信息的成像。南京航空航天大学对基于多普勒反投影的UCW SAR成像方法进行了跟踪[8-9],研究了单频CW双基SAR成像。对比已有的窄带雷达成像方法,基于多普勒反投影的方法不需要平面波前假设和远场近似,适用于宽孔径和大场景,适用于任意接收机和发射机航迹和单/双/多基配置,因此本文着重研究基于多普勒反投影的UCW-SAR动目标成像。

本文首先给出通用的SAR回波模型,适用于窄带发射波形、静止和运动目标、单基和双基SAR配置,以及任意发射和接收机航迹,然后给出NCW-SAR动目标成像方法和目标速度估计方法。最后通过仿真实验,验证NCW-SAR动目标成像方法的性能。

1 SAR回波模型

本节给出SAR回波模型,适用于任意发射波形,在“停-走-停”假设下,该回波模型可退化为常见的传统宽带脉冲SAR回波模型。

假设γT(t),γR(t)为发射机和接收机轨迹,p(t)为发射波形,V(X)=ρ(x)δ(xψ (x))为目标散射率函数。利用波方程,运动目标的接收信号可表示为[12]

$ \begin{array}{l} r\left( t \right) = {E^{sc}}\left( {{\mathit{\boldsymbol{\gamma }}_R}\left( t \right),t} \right) = \\ \int {\frac{{\delta \left( {t - t' - \left| {{\mathit{\boldsymbol{\gamma }}_R}\left( t \right) - \mathit{\boldsymbol{z}}} \right|/{c_0}} \right)}}{{4{\rm{ \mathsf{ π} }}\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( t \right) - \mathit{\boldsymbol{z}}} \right|}}V\left( {\mathit{\boldsymbol{\alpha }}\left( {\mathit{\boldsymbol{z}},t'} \right)} \right)} \times \\ \int {\frac{{\delta \left( {t' - t'' - \left| {\mathit{\boldsymbol{z}} - {\mathit{\boldsymbol{\gamma }}_T}\left( {t''} \right)} \right|/{c_0}} \right)}}{{4{\rm{ \mathsf{ π} }}\left| {\mathit{\boldsymbol{z}} - {\mathit{\boldsymbol{\gamma }}_T}\left( {t''} \right)} \right|}}\ddot p\left( {t''} \right){\rm{d}}t''{\rm{d}}t'{\rm{d}}\mathit{\boldsymbol{z}}} \end{array} $ (1)

式中:Esc(γr(t),t)表示散射场函数;zΓ(Xt)代表运动目标在时间t的位置, X=(x, ψ(x));Γ(∙,t)的反函数Γ- 1(∙,t)表示为α(∙,t);t″为发射信号时间;t ′代表发射信号入射至目标处的时间;t代表接收时间。

考虑起始时刻为s的一小段时间,式(1)可以表示为

$ \begin{array}{l} r\left( {t + s} \right) = \int {\frac{{\delta \left( {t - t' - \left| {{\mathit{\boldsymbol{\gamma }}_R}\left( {t + s} \right) - \mathit{\boldsymbol{z}}} \right|/{c_0}} \right)}}{{4{\rm{ \mathsf{ π} }}\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( {t + s} \right) - \mathit{\boldsymbol{z}}} \right|}}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;V\left( {\mathit{\boldsymbol{\alpha }}\left( {\mathit{\boldsymbol{z}},t' + s} \right)} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\int {\frac{{\delta \left( {t' - t'' - \left| {\mathit{\boldsymbol{z}} - {\mathit{\boldsymbol{\gamma }}_T}\left( {t'' + s} \right)} \right|/{c_0}} \right)}}{{4{\rm{ \mathsf{ π} }}\left| {\mathit{\boldsymbol{z}} - {\mathit{\boldsymbol{\gamma }}_T}\left( {t'' + s} \right)} \right|}}} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\ddot p\left( {t'' + s} \right){\rm{d}}t''{\rm{d}}t'{\rm{d}}\mathit{\boldsymbol{z}} \end{array} $ (2)

zΓ(zt′+s)=X+Vx(t′+s),其中Vx=(vx, ▽xψ(x)∙vx), 表示目标初始时刻速度。代入式(2),经过变换后可得

$ \begin{array}{l} r\left( {t + s} \right) = \\ \int {\frac{{\ddot p\left( {\alpha t - \tau + s} \right)q\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right)}}{{{{\left( {4{\rm{ \mathsf{ π} }}} \right)}^2}\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( s \right) - \left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right)} \right|\left| {\left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right) - {\mathit{\boldsymbol{\gamma }}_T}\left( s \right)} \right|}}{\rm{d}}\mathit{\boldsymbol{x}}{\rm{d}}\mathit{\boldsymbol{v}}} \end{array} $ (3)

式中:V=(v, ▽xψ(x)∙v), q(xv)定义为考虑目标运动的场景反射率函数

$ q\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right) = \rho \left( \mathit{\boldsymbol{x}} \right)\delta \left( {\mathit{\boldsymbol{v}} - {\mathit{\boldsymbol{v}}_x}} \right) \approx \rho \left( \mathit{\boldsymbol{x}} \right)\delta \left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{v}}_x}} \right) $ (4)

式中:ρ(x)表示地面反射率;α为目标和发射机、接收机运动引入的波形尺度因子

$ \begin{align} &\alpha =\frac{1-\left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot {{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)/{{c}_{0}}}{1+\left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot {{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)/{{c}_{0}}}\cdot \\ &\ \ \ \ \ \ \frac{1+\left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot \mathit{\boldsymbol{V}}/{{c}_{0}}}{1-\left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot \mathit{\boldsymbol{V}}/{{c}_{0}}} \\ \end{align} $ (5)

式中:${\mathit{\boldsymbol{\dot r}}_{T\left( R \right)}}$表示发射机/接收机的速度。式(3)中τ表示回波延迟

$ \begin{align} &\tau =\left[ \left| {{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)-\mathit{\boldsymbol{X}} \right|+\left| {{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)-\mathit{\boldsymbol{X}} \right| \right]/{{c}_{0}}- \\ &\ \ \ \ \ \left[ \left( \left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)+\left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right) \right)\cdot \right. \\ &\ \ \ \ \ \left. \mathit{\boldsymbol{V}}s \right]/{{c}_{0}} \\ \end{align} $ (6)

由式(6)可见,相比静止目标成像,动目标的回波延时中增添了由目标速度和雷达视线方向决定的延时项。

对于窄带信号$p\left( t \right) = {e^{{\rm{j}}{\omega _0t}}}\tilde p\left( t \right)$,其中ω0为载波频率,$\tilde p\left( t \right)$为信号复包络,式(3)可进一步写为

$ \begin{array}{l} r\left( {t + s} \right) \approx - \omega _0^2\\ \int {\frac{{\tilde p\left( {\alpha t - \tau + s} \right){{\rm{e}}^{{\rm{j}}{\omega _0}\left( {\alpha t - \tau + s} \right)}}q\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right)}}{{{{\left( {4{\rm{ \mathsf{ π} }}} \right)}^2}\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( s \right) - \left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right)} \right|\left| {\left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right) - {\mathit{\boldsymbol{\gamma }}_T}\left( s \right)} \right|}}{\rm{d}}\mathit{\boldsymbol{x}}{\rm{d}}\mathit{\boldsymbol{v}}} \end{array} $ (7)

对于窄带信号,可以忽略尺度因子引起的复包络变化,在上式中使用了近似$\tilde p\left( {\alpha t} \right) \approx \tilde p\left( t \right)$

2 NCW-SAR动目标成像

本节推导适用于窄带SAR运动目标的成像方法,并给出目标速度估计方法。

2.1 成像前向模型

基于式(7)给出的SAR回波模型,首先建立如下成像前向模型

$ d\left( {s,\mu } \right) = \int {r\left( {t + s} \right){p^ * }\left( {\mu t} \right)\phi \left( t \right){\rm{d}}t} $ (8)

式中:d表示成像数据;μ为引入的尺度因子;s为平移变量,ϕ(t)为时间窗函数。该成像前向模型通过引入尺度因子μ利用发射波形的多普勒频率信息,变量s类似传统宽带SAR成像的慢时间变量,对应孔径位置采样,ϕ(t)决定了每一孔径采样处的参与成像处理的信号长度。将式(7)代入式(8),对α进行近似,d可被写为如下形式

$ d\left( {s,\mu } \right) = \int {{{\rm{e}}^{ - {\rm{i}}\phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right)}}A\left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right)q\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right){\rm{d}}\mathit{\boldsymbol{x}}{\rm{d}}\mathit{\boldsymbol{v}}{\rm{d}}t} $ (9)

式中:ϕA分别为成像前向模型的相位函数和幅度函数

$ \phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right) = 2{\rm{ \mathsf{ π} }}{f_0}t\left[ {\left( {\mu - 1} \right) + {f_d}\left( {s,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right)/{f_0}} \right] $ (10)
$ \begin{array}{l} A\left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right) = \\ \frac{{\tilde p\left( {t - \tau + s} \right){{\tilde p}^ * }\left( t \right){{\rm{e}}^{i{\omega _0}\left( {s - \tau } \right)}}\omega _0^4}}{{\left( {4{\rm{ \mathsf{ π} }}} \right)2\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( s \right) - \left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right)} \right|\left| {\left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right) - {\mathit{\boldsymbol{\gamma }}_T}\left( s \right)} \right|}} \end{array} $ (11)

式中:fd(sxv)为多普勒频率

$ \begin{align} &{{f}_{d}}\left( s,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}} \right)=\frac{{{f}_{0}}}{{{f}_{d}}}\left[ \left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{T}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot \left( {{{\mathit{\boldsymbol{\dot{ }\!\!\gamma\!\!\rm{ }}}}}_{T}}\left( s \right)-\mathit{\boldsymbol{V}} \right)+ \right. \\ &\ \ \ \ \ \ \ \left. \left( \overset\frown{{{\mathit{\boldsymbol{ }}\!\!\gamma\!\!\rm{ }}_{R}}\left( s \right)-\left( \mathit{\boldsymbol{X}}+\mathit{\boldsymbol{V}}s \right)} \right)\cdot \left( {{{\mathit{\boldsymbol{\dot{ }\!\!\gamma\!\!\rm{ }}}}}_{R}}\left( s \right)-\mathit{\boldsymbol{V}} \right) \right] \\ \end{align} $ (12)

对式(9)应用驻留相位原理可知,成像前向模型的相位主导项位于(x, v)四维空间中多普勒频率相同的点构成的等多普勒流形上,表示为F(s, μ)={(x, v):fd(s, x, v)=f0(1- μ)}。

从上述对成像模型的分析可以看出,式(8)构建的成像数据实际上是场景反射系数q(xv)在等多普勒流形F(s, μ)上的投影。由式(12)多普勒频率定义可知,该等多普勒流形取决于发射机和接收机的飞行轨迹和速度,场景中散射点位置和速度。图 1给出发射机和接收机均沿直线飞行,针对给定速度分布v0下的位置空间等多普勒曲线图,Fv0(s, μ)= {x:fd(s, x, v0)=f0(1-μ)}。图 1(a)对应含有动目标场景,假设在[5, 10, 0] km处有一运动目标,运动速度Vx=[100, 50, 0] m/s,实心圆表示动目标的初始位置。图 1(b)对应纯静止场景,实心圆给出与动目标初始位置对应的静止目标位置。图中虚线与空心三角形表示发射机轨迹与发射机位置,实线与实心三角形表示接收机轨迹与接收机位置。发射机和接收机轨迹方程分别为γT=[3.5, 261t, 6.5]km和γR=[261t, -7+261t, 6.5] km,发射机和接收机速度均为261m/s。比较图 1图 2可以看出,场景中存在动目标时的等多普勒曲线与纯静止场景对应的等多普勒曲线分布不同。因此两种情况下获得数据d也不同。

图 1 SAR直线轨迹下的双基导多普勒曲线示意图 Figure 1 Bistatic iso-doppler contours for SAR traversing linear trajec tory

图 2 SAR直线轨迹下的成像几何关系示意图 Figure 2 Imaging geometry for SAR traversing linear trajectory

2.2 成像

成像前向模型给出场景q(x, v)和成像数据d(s, μ)之间的映射关系,但是由两维数据d无法直接反演出四维散射率分布,即两维速度和两维位置组成的四维空间的目标散射率分布。因此,首先基于某一个速度假设vh,采用滤波反投影方法(Filtered back projection, FBP)进行成像,得到该速度假设下的场景散射率分布qvh(z),称之为vh-场景散射率图像。针对一定的速度范围,获得一组vh -场景散射率图像。若速度vh与场景中动目标实际速度相等,则该动目标将被正确重建。下面给出具体的成像过程。

给定速度假设vh,利用FBP进行成像

$ {q_{{v_h}}}\left( \mathit{\boldsymbol{z}} \right) = \int {{{\rm{e}}^{{\rm{i}}{\phi _{{v_h}}}\left( {t,\mathit{\boldsymbol{z}},s,\mu } \right)}}{Q_{{v_h}}}\left( {\mathit{\boldsymbol{z}},t,s} \right)d\left( {s,\mu } \right){\rm{d}}t{\rm{d}}s{\rm{d}}\mu } $ (13)

式中:ϕvh(tzsμ)=ϕ(tzvhsμ)由式(10)给出。Qvh(zts)表示针对给定假设速度vh的滤波器,可按照不同准则来设计。本文按照使点目标扩展函数为Delta函数来设计,具有以下形式[13]

$ {Q_{{v_h}}}\left( {\mathit{\boldsymbol{z}},t,s} \right) = \frac{{{A^ * }\left( {t,\mathit{\boldsymbol{z}},{\mathit{\boldsymbol{v}}_h},s,\mu } \right)}}{{{{\left| {A\left( {t,\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{v}}_h},s,\mu } \right)} \right|}^2}}}\frac{1}{J} $ (14)

式中:1/J为补偿变量代换引入的雅克比行列式因子。该滤波器在幅度上对A(t, x, v, s, μ)进行补偿,在vh=vx时,使重建的动目标具有正确的幅度。

将式(13)与式(9)对比可见,成像过程是对成像数据d的相位进行匹配,同时幅度进行一定调整。相位匹配决定了目标的能量是否可被相干积累,并重建在正确的位置。从反投影角度,式(13)表明成像实际是将成像数据d向等多普勒线Fvh(s, μ)反投影,并将不同孔径采样时刻的投影结果累加。显然,只有当vh=vx时,数据才能被投影到正确的等多普勒线,结果才可以相干累加,得到良好聚焦的动目标图像。

2.3 动目标速度估计

利用反射率图像的对比度进行速度估计,定义对比度图像如下

$ I\left( {{\mathit{\boldsymbol{v}}_h}} \right) = \frac{{\sqrt {M\left[ {{{\left| {{q_{{v_h}}}\left( \mathit{\boldsymbol{z}} \right) - M\left[ {{q_{{v_h}}}\left( \mathit{\boldsymbol{z}} \right)} \right]} \right|}^2}} \right]} }}{{M\left[ {{q_{{v_h}}}\left( \mathit{\boldsymbol{z}} \right)} \right]}} $ (15)

式中:vh=(vh1vh2)为对比度图像的坐标;M[∙]表示空间平均算子;qvh表示速度为vh时的反射率图像。图像对比度实际为重建图像标准差与均值的比值。当vh为目标真实速度,图像聚焦,图像对比度最大,因此由对比度图像峰值可以确定估计速度的大小。如果场景中存在多个不同速度的运动目标,则对比度图像可能有几个峰值,每个峰值对应于不同运动目标的速度。因此可以通过检测对比度图像中的局部最大值来检测不同运动目标的不同速度。具体做法是:首先确定对比度图像峰值,作为第一个检测出的目标速度,然后屏蔽以此速度为中心的一小块区域,可简单置零,再确定新的对比度图像的第二大峰值,以此重复,直至对比度图像没有明显的峰值。

关于动目标速度估计运算量,速度采样空间为N×N时,算法复杂度相对于静止场景成像提高N2倍。静止目标成像算法中的反投影步骤的时间复杂度在整个成像算法中最高[13],反投影步骤的复杂度决定了成像的运算复杂度。但是无论是反投影成像算法,还是速度估计都可以采用并行处理,并使用快速算法[14-15]和图形处理单专用硬件进行加速,因此实际应用的可以确保高的计算效率。

需要指出,速度估计精度受限于速度空间离散化网格的疏密程度,网格越细,估计精度越高。在实际中,为提高目标速度估计精度和确保成像质量,可以第一步采用较大的速度范围和采样间隔,确定目标真实速度附近的小区域,然后第二步再在较小的速度区域内用更细的网格,以获得更准确的速度估计。

2.4 速度分辨率分析

假设目标运动速度vx0,位于场景x0处,则有

$ q\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}}} \right) = \delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_0}} \right)q\left( {{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right) $ (16)

对式(10)在v=0处进行泰勒展开,可以得到

$ \begin{array}{l} \phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right) \approx \phi \left( {t,\mathit{\boldsymbol{x}},0,s,\mu } \right) + \\ \;\;\;\;\;\;\;{\nabla _v}\phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right)\left| {_{v = 0}} \right. \cdot \mathit{\boldsymbol{v}} \end{array} $ (17)

将式(16,17)代入式(9)可得

$ \begin{array}{l} d\left( {s,\mu } \right) \approx \\ \;\;\;\;\;\;\;\;\int {{{\rm{e}}^{ - {\rm{i}}\phi \left( {t,\mathit{\boldsymbol{x}},0,s,\mu } \right)}}{{\rm{e}}^{ - {\rm{i}}{\nabla _v}\phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right)\left| {_{v = 0}} \right.}}} \times \\ \;\;\;\;\;\;\;\;q\left( {{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right)A\left( {t,{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}},s,\mu } \right){\rm{d}}\mathit{\boldsymbol{v}}{\rm{d}}t \end{array} $ (18)

其中

$ \phi \left( {t,{\mathit{\boldsymbol{x}}_0},0,s,\mu } \right) = 2{\rm{ \mathsf{ π} }}t\left[ {\left( {\mu - 1} \right){f_0} + {f_d}\left( {s,\mathit{\boldsymbol{x}},0} \right)} \right] $ (19)
$ {\nabla _v}\phi \left( {t,\mathit{\boldsymbol{x}},\mathit{\boldsymbol{v}},s,\mu } \right)\left| {_{v = 0}} \right. = 2{\rm{ \mathsf{ π} }}t{\nabla _v}{f_d}\left( {s,{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right)\left| {_{v = 0}} \right. $ (20)

fd(sx,0)表示目标速度为0时的多普勒频率。将式(19,20)代入式(18)得到

$ \begin{array}{l} d\left( {s,\mu } \right) \approx \int {\left( {{{\rm{e}}^{ - {\rm{i2 \mathsf{ π} }}t{\nabla _\mathit{\boldsymbol{v}}}{f_d}\left( {s,{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right)\left| {_{\mathit{\boldsymbol{v}} = 0}} \right. \cdot \mathit{\boldsymbol{v}}}}q\left( {{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right)A\left( {t,{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}},} \right.} \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {s,\mu } \right){\rm{d}}\mathit{\boldsymbol{v}} \times {{\rm{e}}^{ - {\rm{i2 \mathsf{ π} }}t\left[ {\left( {\mu - 1} \right){f_0} + {f_d}\left( {s,\mathit{\boldsymbol{x}},0} \right)} \right]}}{\rm{d}}t \end{array} $ (21)

tvfd(sx0v)|v=0ζζ表示与v相关的“频域”傅里叶矢量,其长度和方向分布决定了动目标的速度分辨率。ζ的模长为

$ \begin{array}{l} \zeta = \left| {t{\nabla _\mathit{\boldsymbol{v}}}{f_d}\left( {s,{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{v}}} \right)\left| {_{\mathit{\boldsymbol{v}} = 0}} \right.} \right| = {L_\phi }\left| {D\psi \cdot } \right.\\ \;\;\;\;\;\;\left[ {2\cos \frac{{{\theta _{TR}}\left( {{\mathit{\boldsymbol{X}}_0},\mathit{\boldsymbol{V}},s} \right)}}{2}{{\hat b}_{TR}}\left( {{\mathit{\boldsymbol{X}}_0},\mathit{\boldsymbol{V}},s} \right) + } \right.\\ \;\;\;\;\;\;\left. {\left. {\frac{{{{\mathit{\boldsymbol{\dot \gamma }}}_T}\left( s \right) - {{\left. \mathit{\boldsymbol{V}} \right)}^ \bot }s}}{{\left| {{\mathit{\boldsymbol{\gamma }}_T}\left( s \right) - \left( {{\mathit{\boldsymbol{X}}_0} + \mathit{\boldsymbol{V}}s} \right)} \right|}} + \frac{{{{\mathit{\boldsymbol{\dot \gamma }}}_R}\left( s \right) - {{\left. \mathit{\boldsymbol{V}} \right)}^ \bot }s}}{{\left| {{\mathit{\boldsymbol{\gamma }}_R}\left( s \right) - \left( {{\mathit{\boldsymbol{X}}_0} + \mathit{\boldsymbol{V}}s} \right)} \right|}}} \right]} \right| \end{array} $ (22)

式中:θTR(X0Vs)表示发射机到目标视线方向与接收机到目标视线方向的夹角,即双基角${{\hat b}_{TR}}({\mathit{\boldsymbol{X}}_0}, \mathit{\boldsymbol{V}}, s)$,表示$[(\overset\frown{{\mathit{\boldsymbol{\mathit{\boldsymbol{\gamma}} }}_T}\left( s \right)-\left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right)})+(\overset\frown{{\mathit{\boldsymbol{\mathit{\boldsymbol{\gamma}}}}_R}\left( s \right)-\left( {\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{V}}s} \right)})]$的单位矢量,即双基SAR视线方向;${({{\mathit{\boldsymbol{\dot r}}}_{TR}}\left( s \right)-{\mathit{\boldsymbol{V}}_x})^ \bot }$表示天线与目标相对速度在垂直视线方向的分量。

由式(22)可以看出,速度分辨率与窗函数长度Lϕ成正比,Lϕ的值越大,速度分辨率越高。除此之外,速度分辨率与目标距发射机和接收机距离、发射机和接收机与目标的相对速度、双基角大小有关。速度分辨率与天线距离成反比,与天线和目标的相对速度成正比,与双基角大小成反比。由于雷达作用距离较远,式(22)中后两项的贡献通常较小,因此双基角大小相对于天线和目标的位置对速度分辨率影响更大。此外,速度分辨率也与孔径采样点数有关,采样点数越多,速度分辨率越高。

3 实验仿真

本节对NCW-SAR动目标成像方法与目标速度估计方法进行验证。发射波形采用单频CW,载波频率为800MHz。仿真中选取场景大小为256×256m2,将其离散为128×128个像素,其中[1, 1]和[128, 128]的像素点分别对应场景中的[0, 0] m和[256, 256] m。

速度估计所需的速度采样空间设置为[-10, 10]×[-10, 10] m/s,并以1m/s为采样间隔将整个速度采样空间离散为21×21个采样点。

仿真验证考虑以下两种SAR配置:(1)SAR沿直线轨迹飞行,平台轨迹方程为γ(t)=[8750+vt, 0, 6500] m,速度为261m/s。(2) SAR沿圆形轨迹飞行,平台轨迹方程为γ(t)=[11+11×cos(vt/R), 11+11*sin(vt/R), 6.5] km,其中速度为261m/s,半径R=11 km。

3.1 SAR平台直线轨迹飞行情况下的运动目标成像

仿真成像几何关系如图 2所示。黑色区域为仿真场景,SAR平台沿X轴正方向匀速直线飞行。运动目标在成像开始时刻位于场景中心,对应像素[65, 65],箭头方向为目标运动方向,运动速度vx=[6,-5] m/s。雷达工作在正侧视模式下,合成孔径长度为5500m,采样点数为2048。场景中心到合成孔径中心的水平距离为11km。

对于每个速度采样值,按照2.2节给出的方法进行成像,计算相应的对比度函数,获得对比度图像,如图 3(a)所示。图中峰值处由白色空心圆标示,对应速度为[6,-5] m/s,与目标实际速度相同。对比度图像峰值所在单元的垂直与水平方向的剖面局部放大图分别如图 3(b, c)所示。

图 3 SAR直线轨迹下单个运动目标速度估计 Figure 3 Velo city estimation of single moving target for SAR traversing straight line ar trajectory

vh=[6, -5] m/s=vx,成像结果如图 4(a)所示。图 4(b)为成像结果沿X方向的剖面图,图 4(c)为沿Y方向的剖面图。可以看出,当速度vh与实际速度vx相等时,峰值点位于(65, 65),目标被重建于正确位置,而且聚焦良好。从而验证了上述UCW-SAR动目标成像方法的有效性。

图 4 SAR直线轨迹下单个运动目标的成像结果(vh=vx=[6, -5] m/s) Figure 4 Imaging results of single moving ta rget for SAR traversing straight linear trajectory (vh=vx=[6, -5] m/s)

如果令速度vh=[5.5,-5] m/s≠v x,即成像所用速度与目标真实速度有偏差,得到目标重建结果如图 5(a)所示,峰值点位于(68, 89)。与图 4结果比较可知,当速度估计出现误差,成像结果峰值点并不是目标真实位置,而且图像散焦,对比度下降。

图 5 SAR直线轨迹下单个运动目标的成像结果(vh=[5.5,-5] m/s≠vx) Figure 5 Imaging result of single moving target for SAR traversing linear trajectory (vh=[5.5, -5] m/s≠vx)

3.2 SAR平台圆轨迹飞行情况下的运动目标成像

图 6所示给出仿真成像几何关系示意图。合成孔径为半径11 km的圆,对应合成孔径长度为2π×11 km,采样点数为4096。运动目标设置与3.1节中配置相同,即运动目标开始时刻位于场景中心,速度vx=[6,-5] m/s。

图 6 SAR圆轨迹下的成像几何关系示意图 Figure 6 Imaging geometry for SAR traversing c ircular traje ctory

按照第2节给出的成像方法获得对比度图像,如图 7(a)所示。图中峰值处由白色空心圆标示,对应速度为[6,-5] m/s,与目标实际速度相同。令vh=[6, -5] m/s=vx,对目标进行成像,成像结果如图 7(b)所示。目标重建于正确位置,而且聚焦效果良好。由于积累了更长的孔径,成像分辨率和成像质量优于3.1节直线轨迹下的成像结果。

图 7 SAR圆轨迹下单个运动目标的成像结果 Figure 7 Imaging results of single moving target for SAR traversing circ ular trajectory

4 结束语

本文首先给出了非“停-走-停”假设下,适用于静止目标和运动目标的SAR回波信号模型,然后给出了基于多普勒反投影的NCW-SAR动目标成像方法,可同时给出速度估计。对NCW-SAR的动目标速度分辨率进行了理论分析。最后通过不同发射机和接收机轨迹下的仿真实验,验证了本文NCW-SAR动目标成像方法的有效性。在后续工作中将对NCW-SAR动目标成像的速度分辨率的理论分析结果进行实验验证,并深入探讨多目标的成像性能,以实际外部照射源信号为发射源的NCW-SAR动目标成像性能。

参考文献
[1] MENSA D, HALEVY S, WADE G. Coherent Doppler tomography for microwave imaging[J]. Proceedings of the IEEE, 1983, 71(2): 254–261. DOI:10.1109/PROC.1983.12563
[2] SUN H, FENG H, LU Y. High resolution radar tomographic imaging using single-tone CW signals[J]. 2010 IEEE Radar Conference, 2010, 29(16): 975–980.
[3] WU Y, MUNSON D C. Multistatic synthetic aperture imaging of aircraft using reflected television signals[C]//Aerospace/Defense Sensing, Simulation, and Controls.[S.l.]: International Society for Optics and Photonics, 2001.
[4] WANG L, YAZICI B. Passive ima ging of moving targets exploiting multiple-scattering using sparse distributed apertures[J]. Inverse Problems, 2012, 28(12): 5009.
[5] GARRY J L, BAKER C J, SMITH G E, et al. Doppler imaging for passive bistatic radar[C]// 2013 I EEE Radar Conference. Ottawa, ON, Canada: IEEE, 2013:1-6.
[6] BARANIUK R G. Compressive sensing[J]. IEEE Signal Processing Magazine, 2007, 24(4): 118–121. DOI:10.1109/MSP.2007.4286571
[7] 刘玉春, 王俊, 杨杰, 等. 基于单频连续波的无源雷达成像研究[J]. 电子与信息学报, 2013, 35(5): 1108–1113.
LIU Yuchun, WANG Jun, YANG Jie, et al. Research of passive radar imaging ba sed on single-frequency continuous wave[J]. Journal of Electronics and Information Technolo gy, 2013, 35(5): 1108–1113.
[8] WANG L, YAZICI B. Bistatic synthetic aperture radar imaging using ultra-narrowband continuous waveforms[J]. IEEE Transactions on Image Processing, 2012, 21(8): 062–067.
[9] 汪玲, 张璇, 韦渊斐. 基于多普勒分辨的连续波合成孔径雷达成像方法[J]. 系统工程与电子技术, 2014, 36(6): 1057–1061. DOI:10.3969/j.issn.1001-506X.2014.06.07
WANG Ling, ZHANG Xuan, WEI Yuanfei. Doppler resolution based im aging method for synthetic aperture radar using continuous waveforms[J]. Systems Engineering and Electronics, 2014, 36(6): 1057–1061. DOI:10.3969/j.issn.1001-506X.2014.06.07
[10] 徐浩, 尹治平, 刘畅畅, 等. 基于压缩感知的稀疏无源雷达成像[J]. 系统工程与电子技术, 2011, 33(12): 2632–2630.
XU Hao, YI Zhiping, LIU Changchang, et al. Space passive radar ima ging based compressive sensing[J]. Systems Engineering and Electronics, 2011, 33(12): 2632–2630.
[11] 张馨文, 王俊. 基于多电视台子孔径综合的无源雷达成像算法[J]. 电子与信息学报, 2007, 29(3): 528–531.
ZHANG Xinwen, WANG Jun. Passive radar imaging algorithm based on sub-apertures synthesis of multiple televisio n stations[J]. Journal of Electronics and Information Technology, 2007, 29(3): 528–531.
[12] WANG L, YAZICI B. Bistatic synthetic aperture radar imaging of moving targets using ultranarrowband continuous waveforms[J]. SIAM Journal on Imaging Sciences, 2014, 7(2): 824–866. DOI:10.1137/130906714
[13] WANG L, YARMAN C E, YAZICI B. Doppler-Hitchhiker: A novel passive synthetic aperture radar using ultranarrowband sources of opportunity[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(10): 3521–3537. DOI:10.1109/TGRS.2011.2142420
[14] ULANDER L M H, HELLSTEN H, STENSTROM G. Synthetic-aperture radar processing using fast factorized back-pro je ction[J]. IEEE Transactions on Aerospace & Electronic Systems, 2003, 39(3): 760–776.
[15] CANDÈS E, DEMANET L, YING L. Fast computation of Fourier integral operators[J]. Siam Journal on Scientific Computing, 2006, 29(6): 2464–2493.