基于GPS-R的双基地SAR移动目标成像方法研究

何振宇， 杨扬， 陈武， 翁多杰

1. 香港理工大学 土地测量及地理资讯学系, 香港 999077;

2. 香港理工大学 深圳研究院, 深圳 518057;

3. 河海大学 计算机与信息学院, 南京 211100

遥感技术和导航定位技术是当前最具应用价值和发展潜力的时空信息采集技术。两种技术的相互融合是遥感和导航综合应用的必然发展趋势（万玮 等，2016；
陈锐志 等，2019）。全球导航卫星系统反射GNSS-R（Global Navigation Satellite System Reflectometry）就是一种典型的融合应用。国内外的学者已经开展了一系列关于GNSS-R技术的理论与应用研究，例如DDM 仿真（夏俊明 等，2019；
Beltramonte等，2020）、海面测高（Xie等，2018）、风速反演（王笑蕾 等，2021）、雪深反演（王泽民 等，2018；
Yu 等，2019）、土壤水分反演（彭学峰 等，2017）、土壤盐度反演（Wu等，2019）以及目标探测（Ma等，2018；
Pastina等，2018）等。

在目标探测领域，GNSS-R 可看作是一种无源雷达系统。由于没有专门的发射机设备，GNSS-R用于目标探测具有低成本、低功耗、隐蔽性强等诸多优点（He 等，2020）。近年来，利用GNSS-R做双基地合成孔径雷达GNSS-SAR（GNSS-based Bistatic Synthetic Aperture Radar）引起了学者们的广泛关注（Antoniou 等，2007；
Liu 等，2013；
Ma 等，2015；
Zeng等，2016；
Liu等，2019）。当前GNSSSAR 成像的主要对象是地表静止目标，而关于动目标成像的研究非常之少。为了扩大GNSS-SAR的应用范围，本文致力于动目标成像方法研究。王博等（2019）对GNSS-SAR 动目标成像做了技术探讨并提出了一种可行的成像方法，其考虑的是运动目标具有径向速度的情况；
此时，目标回波的多普勒频率恒定，可直接通过频域滤波分离动、静目标并根据多普勒频率的频移量估计目标速度，方便之后的成像。但是，如果目标具有切向速度，则目标回波的多普勒频率不再恒定而是随时间变化；
此时，频域滤波法不再适用。此外，该方法成像使用的时域后向投影算法适用于各种双基地几何构型并且没有近似运算使得成像精度高，但是，需要对成像场景中每个点进行计算，导致运算量极大不利于GNSS-SAR 系统的实时性应用（Zeng等，2016）。

本文以全球定位系统GPS（Global Positioning System）卫星作为机会照射源，针对具有切向速度的目标，采用频域GNSS-SAR 算法成像动目标。频域算法虽然有一定的近似运算，但其最大的优点在于高效的信号处理效率，适用于实时处理（Antoniou 和Cherniakov，2013）。前期工作（何振宇 等，2020）已经初步探讨了成像海面移动船只的可行性；
本文将继续扩展之。频域GNSS-SAR成像算法主要包含两个步骤：距离单元徙动校正和方位向压缩。然而，对于未知运动速度的目标，传统的GNSS-SAR 算法难以完成上述两个步骤，SAR 图像上的动目标将会发生散焦、位置偏移等问题。因此，动目标成像需要解决：未知的距离单元徙动校正、目标运动速度估计以及方位向压缩。为了解决以上问题，本文首先构建了成像场景的双基地几何模型如图1所示。其中，接收机固定于海岸，同时接收GPS 卫星信号以及目标回波信号；
在观测时间内GPS 卫星处于准静态状态；
移动目标的运动轨迹垂直于接收机天线视线方向以便形成合成孔径。基于该双基地几何模型，推导了一种近似的双基地距离历程，方便描述目标回波的方位向相位变化。然后，采用Keystone变换（Perry等，1999；
王娟和赵永波，2011）补偿未知的距离单元徙动以及组合短时傅里叶变换STFT（Short Time Fourier Transform）（Auger等，2012）和随机抽样一致性RANSAC（Random Sample Consensus）（Fischler 和Bolles，1981）算法估计目标运动速度。最后，推导了方位向匹配滤波器完成方位向压缩从而实现动目标成像。现场实验将用于验证本文提出的动目标成像方法。

图1 动目标成像的双基地几何示意图Fig. 1 Bistatic acquisition geometry for moving target imaging

本节首先介绍动目标SAR 成像的双基地几何从而推导近似的双基地距离历程；
然后，给出GPS直达波和目标回波的雷达信号模型。

2.1 双基地几何模型

如图1 所示，GPS 接收机被固定在海岸边。它包含两个通道分别命名为参考通道和监视通道；
参考通道天线负责接收GPS卫星发射的直达波作为参考信号；
监视通道天线则接收移动目标反射的回波。为了方便描述，图2 对图1 的双基地几何做了抽象简化。在图2（a）中，移动目标Tg 和接收机Rx处在相同的高度（即地表）并且假设在整个观测时间内卫星Tx处于准静态状态（例如，观测时长小于120 s）以及移动目标Tg 始终处在监视通道天线3 dB 波束角照射范围内。卫星Tx 到目标Tg 距离表示为Rt，卫星Tx 到接收机Rx 距离表示为Rb，目标Tg 到接收机Rx 距离表示为Rr。定义双基地距离Rbi为直达波与目标回波的传播距离之差，表示为

式中，u∈[1，T]称为方位时间或者慢时间；
T表示整个观测时长。式（1）可以精确地描述观测时间内双基地距离历程；
但是，它不便于下一节的动目标成像处理。因此，本节推导一种近似的双基地距离历程。

考虑到GPS卫星到地表距离有20000 —24000 km（刘立东 等，2004），远远大于目标到接收机距离（例如小于10 km）。所以，图2（a）中卫星到目标距离以及卫星到接收机距离近似于平行关系。根据三角函数关系可知，式（1）可近似为

式中，θ表示双基地角。杨振起等（1998）指出双基地角θ具有如下角度关系：

式中，α表示卫星高度角；
β是图2（a）中由接收机、目标和星下点构成的角度。从图2（b）可以看出，β角有如下角度关系：

式中，φ= cos-1(RsRr)表示监视通道天线相位中心对目标的斜视角；
Rs为目标到接收机垂直距离；
az定义为本地卫星方位角（将在4.1 节结合具体实验场景进行讨论）；
ε是图2（b）中由接收机、星下点和目标构成的角度。然而，ε几乎可以忽略不计，因为星下点到接收机距离和星下点到目标距离依然远远大于目标到接收机距离。由此，结合式（3）和（4），式（2）可表示为

式中，u∈[1，2L v]；
L表示合成孔径长度的一半；
v表示目标移动速度。然而，如果图2（b）中的目标自左向右移动，则双基地距离历程变为

式（5）和（6）表明，当目标的运动方向不同时，对应的双基地距离历程表示式也不相同。需要注意的是，图2（b）仅显示了星下点在雷达视线方向右侧的例子，如果星下点位于左侧并且目标依然自右向左运动，则双基地距离历程表示式变为式（6）。综上可知，在目标穿过雷达视线方向之前，如果目标和星下点位于视线方向的同侧，则双基地距离历程表示式为式（5），反之，则为式（6）。该特点可用于判断目标移动方向，将在4.2 节进行讨论。

图2 不同视角的双基地几何示意图Fig. 2 Different views of bistatic acquisition geometry

利用仿真实验来评估本节推导的近似双基地距离适用范围，即，计算式（1）定义的双基地距离与式（5）或式（6）之间的最大误差。表1列出了仿真所需参数；
图3 给出了在不同的卫星高度角、本地卫星方位角条件下，两种双基地距离的最大误差与目标到接收机垂直距离之间的关系。从中可以看出，最大误差随着高度角、方位角以及垂直距离的增大而增大。然而，以图3（b）为例，最大误差在大约3000 m 范围内没有超过GPS L1信号载波波长（约0.19 m）。因此，在3000 m范围以内，推导的近似双基地距离历程可安全地用于接下来的动目标成像处理。

图3 不同卫星高度角、本地卫星方位角和目标到接收机垂直距离下两种双基地距离之间的最大差值Fig. 3 Maximum errors between two bistatic ranges versus the target-to-receiver vertical ranges in the case of different satellite elevation angles and local satellite azimuth angles

表1 仿真参数Table 1 Simulation parameters

2.2 信号模型

本文选用的GPS L1信号SL1表示为

式中，t表示观测时刻；
A表示信号幅度；
C表示测距码；
D表示导航电文；
fc为载波中心频率；
φ0为信号初相。GPS-R 做目标探测需要直达波作为参考信号，然而，参考通道接收的直达波能量很低，小于-30 dB（Antoniou 和Cherniakov，2013）。这样的直达波不能直接用于回波信号处理，需要做信号同步。信号同步通过捕获和跟踪直达波来提取信号中的码相位、载波相位、多普勒频率、导航电文等观测量，再利用这些观测量生成接近于无噪声的直达波作为参考信号（Antoniou 等，2015）。将提纯后的直达波Sdi表示为雷达信号格式且忽略信号的幅值和初相：

式中，t∈[1，PRI]称为快时间；
PRI表示脉冲重复间隔（即测距码持续时间为1 ms）。τdi、fdi和φdi分别表示码延迟、载波频率和载波相位。他们的取值依赖于卫星与接收机之间的相对位置和运动。此外，除了传播距离产生的码延迟和载波相位延迟，τdi和φdi还包含了由大气延迟和接收机钟差所引起的码延迟误差和载波相位误差。同理，考虑单个点目标，监视通道内的目标回波Sre表示为

式中，τre、fre和φre分别表示回波信号的码延迟、载波频率和载波相位。

图4 为本文提出的GPS-R 双基地SAR 动目标成像方法流程图。整个过程可划分为3个部分，分别为预处理、速度估计和方位向压缩。下面将分别描述这3个部分。

图 4 提出的动目标成像方法流程图Fig. 4 Flowchart of the proposed moving target image formation algorithm

3.1 预处理

2.2 节已经给出了GPS 直达波和目标回波的信号模型。可以看出，GPS信号是连续波，不能反映SAR 成像所需要的目标回波在二维（距离，方位时间）域上的情况。此外，导航电文会引起信号相位跳变，破坏载波多普勒频率的连续性以及GPS信号内的各种误差项有待消除（王博 等，2019）。为了解决上述问题，需要做距离向压缩处理。距离向压缩即为式（8）的参考信号和式（9）的目标回波在快时间域上做互相关运算，得到二维距离压缩信号rc，表示为

式中，τ∈[0，PRI- 1]表示测距码偏移量，*表示取共轭操作，Δf表示快时间域两路信号的载波频率差，Δφ表示慢时间域两路信号的相位差。需要指出的是：（1）D(t-τdi(u))×D(t-τre(u))乘积为1，即导航电文被消除，因为在6000 km 范围内参考信号和目标回波的导航电文相同 （Zeng，2013）；
（2）Δf→0，因为相比于卫星速度，目标速度在快时间域内引起的载波多普勒频率变化几乎忽略不计（Pastina等，2018）；
（3）Δφ消除了大气延迟和接收机钟差所产生的误差相位，仅剩由双基地距离产生的传播延迟相位（He 等，2020），因为接收机和目标上方的大气延迟误差几乎相同（Ma等，2018）以及接收机双通道共用时钟使得钟差相同。综上，积分后，式（10）表示为

式中，CF表示互相关函数包络；
c为光速；
λ为载波波长。

从式（11）可以看出，在方位时间域上双基地距离历程会造成互相关函数包络平移即距离单元徙动以及复载波相位变化即多普勒频率徙动。多普勒频率徙动与目标的运动速度有关，将在3.2 节用于速度估计。执行距离单元徙动校正需要已知目标的运动速度。但是，对于非合作的运动目标其速度信息是未知的。Keystone 变换可以在目标速度未知的情况下消除距离单元徙动。将2.1 节推导的近似双基地距离历程式（5）或（6）代入式（11），结合王娟和赵永波（2011）给出的Keystone变换快速实现算法即可完成距离单元徙动校正。因此，这里不再赘述具体的实施步骤。图5展示了距离单元徙动校正前和校正后的二维距离压缩信号。从中可以看到，图5（a）中的互相关函数包络发生了明显的包络偏移，这将不利于之后的方位向信号处理；
而经过Keystone 变换校正后，图5（b）中的互相关函数包络已经对齐在相同的距离单元。以目标在合成孔径中心位置作为参考位置，距离单元徙动校正后，式（11）可重写为

图5 Keystone变换校正前后的二维距离压缩信号Fig. 5 2D Range-compressed signal before and after the Keystone transform correction

式中，Rc=Rs(1 + cosαcosaz)表示目标在合成孔径中心位置时的双基地距离。由此，目标到接收机垂直距离Rs可表示为

垂直距离Rs将会用于接下来的速度估计和方位向匹配滤波器生成。

3.2 速度估计

目标运动在方位向上给回波信号带来相位调制。在设计方位向匹配滤波器前，需要已知目标的运动速度（Kulpa 等，2013）。本文假设目标以恒定的速度移动，例如以巡航速度航行的船只。因此，本节处理目标速度估计问题。

将式（5）或（6）代入式（12）并取二阶泰勒级数展开，目标回波的方位向相位Φ表示为

目标回波的多普勒频率fd可由式（14）取一阶导数得到：

式中，γ表示调频率；
f0表示多普勒质心。式（15）表明目标回波在方位向上表现为线性调频信号并且目标速度的平方正比于调频率绝对值。因此，可以把速度估计问题转换为对调频率的估计。时频分析技术是一种有效的估计线性调频信号参数的方法，例如，Radon-Wigner 变换（Wood 和Barry，1994）、chirp-Fourier 变换（Xia，2000）、分数阶Fourier 变换（Ozaktas 等，1996）等。这些时频分析方法的计算复杂度阶数为Ο(N2logN)。在GNSS-SAR 中，由于目标回波能量较低，长孔径时间非常有必要（Antoniou 等，2015）。一方面可以提升方位向分辨率，另一方面可以增强目标回波的信噪比。然而，长孔径时间意味着大时间带宽积，导致上述时频分析方法非常耗时（保铮 等，2005）。因此，这里提出一种组合STFT 和RANSAC 的调频率估计方法，其计算复杂度阶数为Ο(NlogN+M)，原因在于STFT 采用快速傅里叶变换实现，其计算复杂度阶数为Ο(NlogN)，而RANSAC 需要若干次的迭代运算，其计算复杂度为Ο(M)。

利用STFT 将式（12）沿方位时间投影到二维（时间，频率）域。目标回波的STFT表示为

式中，ω∈[ - PRF 2，PRF 2]表示多普勒频率；
PRF 表示脉冲重复频率（即2.2 节PRI的倒数）；
η∈[ -h2，h2]；
winh表示长度为h的短时窗函数。假设有5 个散射点依次分布在船只长度方向上，图6（a）给出了散射点回波经STFT 投影到时频域的仿真结果。从中可以看出，由多个散射点回波组成的线性调频信号在时频域上表现为一条有斜率的宽直线，其斜率即为调频率，可以通过拟合直线的方法来估计斜率。因此，首先需要将直线上的时频点提取出来。这里通过设定阈值Γ做二值化，用阈值Γ将时频点分为大于Γ的集合和小于Γ的集合，而图6（b）则为二值化后的结果。

图6 5个散射点回波的时频分布和对应的二值图像Fig. 6 Time-frequency distribution of five scattering points and the corresponding binary image

由于散射点回波的信噪比有限，除了直线上的时频点（称为内点），图6（b）中还包含有大量的噪声点（称为离群点）。大量的离群点使得传统的最小二乘法LSM（Least Squares Method）不再有效，这是由于LSM 将所有的时频点参与拟合问题。这里采用计算机视觉常用的RANSAC 算法实现在大量离群点中拟合直线。对传统的RANSAC 算法做了修改，具体步骤如下：

步骤1：将大于阈值Γ的时频点集合命名为S并定义Ns为S内元素总数；
令K为最大迭代次数。

步骤2：从S 中随机选择两个时频点，利用式（15）构建拟合直线。由于船只速度有限，设置斜率范围来约束随机选点过程从而增强算法的鲁棒性。

步骤3：对于任意的时频点(ui，ωi)，计算该点到步骤2拟合直线距离di，表示为

式中，|•|表示取绝对值；
上标k表示第k次迭代。定义thtol为时频点到拟合直线的最大距离（简称为容许阈值），如果di≤thtol，则时频点(ui，ωi)被认为是内点。

步骤4：统计内点的数目并记为Nk。定义thscale为Nk Ns的最小值（简称最小比例），如果Nk Ns≥thscale，则把该拟合直线的斜率放入候选解集合C中。

步骤5：让k=k+ 1。如果k≤K，则返回步骤2。否则进入步骤6。

步骤6：使用中值滤波器从候选解集合C 中确定估计的调频率，这是由于图6（b）中的直线有一定的宽度，会有若干组符合拟合条件的调频率。中值滤波解可以提高估计的准确性和鲁棒性。最后，利用式（15）的调频率公式和式（13）的垂直距离公式计算目标速度。

3.3 方位向压缩

方位向压缩即为在方位频率域中利用匹配滤波器来消除目标回波的调制相位项从而实现脉冲压缩。因此，本节需要推导方位向匹配滤波器。

将式（5）或（6）代入式（12）并在方位时间域做傅里叶变换：

式中，σ表示频谱的幅值；
ψ是复指数相位，表示为：

式中，第一项为调制相位项，后两项分别为线性项和常数项。因此，方位向匹配滤波器H0即为调制相位项的共轭形式：

式（21）根号内的“±”号表明当目标的运动方向不同时，对应的方位向匹配滤波器也不相同；
换句话说，利用两种方位向匹配滤波器以及卫星与接收机之间的方位关系可以判断目标的移动方向，下节实验将进行讨论。由式（21）的组成可知，生成方位向匹配滤波器需要卫星高度角、卫星方位角、监视通道天线视线方位角、目标到接收机垂直距离以及目标速度，其中前三个参数可直接获取，后两个参数分别在3.1 节和3.2 节得到；
最后，将式（19）乘以式（21），再利用逆傅里叶变换将乘积结果转换回方位时间域即可完成方位向压缩。最后，成像结果I表示为：

式中，x=(c×τ)( 1 + cosαcosaz)表示距离单元；
y=v×u表示横向距离单元；
xTg和yTg表示目标在（距离，横向距离）域的坐标。

4.1 实验场景

2019 年5 月16 日，在香港数码港海滨公园进行了海面船只探测实验。图7（a）为双通道GPS接收机前端。关于接收机前端的详细信息可参考Yang等（2019）。图7（b）中的圆天线和方天线分别连接接收机前端的参考通道和监视通道，其中，圆天线为商用现成的右旋圆极化天线用于接收GPS直达波，而方天线为定制的左旋圆极化天线，用于接收海面船只反射的目标回波。图8为现场实验的俯视图示意，以正东方向作为横轴和正北方向作为竖轴建立坐标系，接收机Rx位于坐标系原点。从图8 标记的方位角关系可知，式（4）定义的本地卫星方位角az在该现场实验场景下为：

图7 现场实验硬件设备Fig. 7 Hardware devices at the field test

图8 现场实验俯视图示意Fig. 8 Top view of the field test

式中，azsat代表GPS 卫星方位角（以正北顺时针方向开始）；
azant等于方天线视线方位角减去180°。

4.2 实验结果与分析

图10为两艘船只的时频分布图，其中分别绘制了RANSAC、LSM以及真值拟合直线。从图10（a）可以看出，在大量的离群点下，LSM 已经失效；
其拟合线（虚线）的斜率为正数与式（15）的调频率公式相悖；
然而，RANSAC 方法的拟合线（点虚线）则非常贴近真值拟合线（实线），这意味着估计的目标速度将非常接近于真实速度。在图10（b）中，由于该目标表面的散射率起伏较大，特别是在零多普勒附近时频点较多而远离零多普勒的两端时频点较少，造成RANSAC 拟合线与真值拟合线间存在一个小的夹角；
而对于LSM方法，虽然该组实验的离群点很少，但是由于离群点恰好集中在图10（b）的下半部分，导致LSM的拟合结果偏离了真值。表4列出了RANSAC方法估计的两艘船只的调频率及目标速度，从中可以看到，两组实验的最大速度估计误差不超过0.6 m/s。考虑到目标表面的散射状态起伏，估计值与真值之间的误差是可以接受的。

图10 两艘船只的拟合结果Fig. 10 Fitting results for two ships

图11（a）和11（b）是用表4 的估计速度生成的SAR 图像。图像中的亮斑即为两艘船只表面散射的回波能量。由于GPS L1 信号的距离向分辨率很差，图像中的亮斑在距离轴方向上表现为宽主瓣。两个亮斑的最大峰值对应在距离轴上分别位于1609.0 m 和766.8 m，非常接近于表3 列出的船只离岸垂直距离并且误差远小于GPS L1 信号的最佳距离向分辨率（约150 m），表明了本文方法估计垂直距离的有效性和准确性。长孔径时间使得GNSS-SAR在方位向上有很高的分辨率。因此，通过测量横向距离轴方向上的船长度来进一步验证本文方法。图11（a）和11（b）已经标注了测量结果分别为254 m 和185 m。对比表3 列出的真实船长度可知，估计的船长度误差不超过真实船长的10%。这意味着可以利用估计的船长度来初步识别雷达观测区域内过往船只类型。此外，作为对比，图11（c）和11（d）给出了使用估计错误的速度（如3 m/s）生成的SAR图。可以看出，图像上的船目标已经发生了严重的散焦以及位置偏移。

表3 船舶的AIS信息和选取的GPS卫星Table 3 AIS information for two cargo ships and the selected GPS satellites

图11 两艘船的SAR图像以及散焦的SAR图像Fig. 11 SAR images and defocused SAR images of two ships

表4 估计结果与AIS真值Table 4 Estimated results and AIS ground truth

本文3.3 节的方位向匹配滤波器推导表明：利用两种方位向匹配滤波器以及卫星与接收机之间的方位关系可以判断目标的移动方向。这里以万海506为例来说明。首先，标记式（21）根号中第二个圆括号内取“+”号时为“AMF+”，取“-”号时为“AMF-”；
分别用两个匹配滤波器生成对应的SAR 图像，其中，“AMF-”对应图11（a）而“AMF+”对应图12。然而，图12 中的亮斑“断裂”成了两块，所以，“AMF-”是正确的方位向匹配滤波器，其由式（6）得到。根据2.1 节可知，在目标穿过雷达视线方向之前，式（6）对应目标和星下点位于雷达视线方向的不同侧。紧接着，从表2 和表3 给出的方天线视线方位角和PRN3 卫星方位角可知，星下点在方天线视线方向的下半侧，即图8星下点位置所示。综上可知，目标的移动方向为自西北向东南移动；
从接收机视角来看，目标正自右向左移动，符合图9（a）箭头标记的运动方向。需要说明的是：如果式（21）根号中az 等于0°，即准单站几何配置，此时，两个方位向匹配滤波器表示式相同，则无法判断目标移动方向。

表2 现场实验和信号处理参数Table 2 Field test and signal processing parameters

图9 海面船只照片Fig. 9 Photographs of two ships on the sea surface

图12 使用错误的方位向匹配滤波输出的SAR图像Fig. 12 SAR image with wrong azimuth matched filter

本文以GPS信号作为机会照射源，提出了一种动目标成像方法。该方法利用Keystone变换校正未知的距离单元徙动；
组合STFT 和RANSAC 算法估计未知的目标运动速度以及推导了方位向匹配滤波器实现移动目标成像。最后，采集了两组现场实验数据来验证本文提出的方法。实验结果显示，在大量离群点下，组合STFT 和RANSAC 算法依然能够估计出船只速度，表明了该估计算法的有效性和鲁棒性；
SAR 图像上的目标特征能够反映船只到接收机垂直距离以及船只长度，表明了该成像算法的有效性。此外，通过对比两种方位向匹配滤波器的成像图，还可以判断船只移动方向。本文方法要求目标运动轨迹垂直于雷达天线视线方向，该要求在河道场景下易于实现。考虑到GNSS-SAR 使用的L 波段可在雨雪、大风、雾霾等恶劣环境中工作并且接收机成本较低，因此，本文提出的成像方法可发展应用于河道监测，估计过往船只的速度、长度、离岸垂直距离以及移动方向。

本文没有考虑GPS 卫星星座的多卫星优势，即多基地雷达系统。融合多颗卫星信号生成的SAR 图像不但可以提高距离向分辨率还可以提高SAR 图像质量有助于目标特征的分类与识别。因此，下一步工作将考虑多基地图像融合问题。