基于高分辨率合成孔径雷达影像的角反射器绝对定位

摘要

详细介绍了合成孔径雷达（SAR）影像大地测量的定位公式和解算流程，利用加权质心法实现了SAR影像中点目标精细定位。针对同轨SAR影像因空间基线过短造成的定位精度和稳定性差问题，提出了一种添加高程限制条件的解决方法。采用数景高分辨率TerraSAR-X和Sentinel-1 SAR影像，对布设于武汉和上海的6个角反射器分别进行了绝对定位实验。结果表明，利用TerraSAR-X影像获得的角反射器内符合定位精度可达到厘米级，利用Sentinel-1影像获得的角反射器内符合定位精度可达到分米级。

关键词

合成孔径雷达; 角反射器; 绝对定位; TerraSAR-X; Sentinel-1A/B

星载合成孔径雷达（synthetic aperture radar，SAR）是20世纪50年代发展起来的一种主动式的微波成像系统，具有全天候、全天时、全球覆盖和高分辨率等独特优势。与光学遥感影像不同，SAR影像中同时包含相位和后向散射强度两类信息。时序干涉SAR测量技术（multi-temporal SAR interferometry，MTInSAR）利用同一地区多景重复轨道SAR影像间的干涉相位时间序列提取高精度、大面积的连续地表高程和形变信息，已广泛应用于各类地表形变监测中^［1-3］。尽管MTInSAR技术能够获取毫米级的视线向形变信息，但在SAR影像中相干点的定位精度只有米级或更差，严重限制了该技术在城市大型基础设施结构形变监测中的应用。为了正确解译MTInSAR技术获取的地面形变信息，基于SAR影像强度信息的像素点绝对定位问题一直备受关注。Curlander^［4］最早根据SAR卫星成像几何提出了利用距离‒多普勒（Range-Doppler，R-D）模型和地球椭球模型的SAR影像地理编码方法，并采用SEASAT SAR影像进行实验，但定位精度只有200 m。Schwabisch^［5］采用数字高程模型（digital elevation model，DEM）优化了地球椭球模型，有效改善了R-D方法的定位精度。其后，周金萍^［6］、杨杰^［7］、袁修孝^［8］等结合RadarSAT、ERS1、ENVISAT ASAR等卫星影像验证了R-D方法的定位效果，但由于早期SAR卫星轨道精度较差和系统测距精度不高等原因，定位精度只有米级。

随着星载SAR技术的发展，TerraSAR-X、COSMO-SkyMed和Sentinel-1A/B等卫星的轨道精度、测距精度和影像分辨率有了较大幅度的提升，为SAR影像中点目标的高精度绝对定位提供了数据支撑。2011年，Eineder等^［9］率先提出了SAR影像大地测量（SAR imaging geodesy）新技术，并利用多景SAR影像基于R-D模型成功确定了SAR影像中目标点的精确二、三维坐标。Blass^［10］、Gisinger^［11］、Dheenathayalan^［12］等分别利用X波段高分TerraSAR-X和C波段EnviSAT、Sentinel-1A/B卫星的多景SAR影像，对人工布设的角反射器（corner reflector，CR）和自然分布的永久散射体（permanent scatterers，PS）分别进行了三维定位实验，结果表明SAR影像目标点的定位精度与影像数量、目标在影像中的信杂比、空间基线大小、斜距大气改正、地表运动改正（固体潮、板块运动等）等因素有关，在理想情况下，利用多景异轨TerraSAR-X卫星影像可以获得目标点厘米级的定位精度，利用Sentinel-1A/B卫星影像可以获得目标点分米级的定位精度。我国学者在SAR影像点目标三维高精度绝对定位方面的研究较少，相关研究主要针对SAR影像的地理编码，定位精度大多在几米甚至几十米^［13-15］。为此，本文基于布设在我国境内的地面角反射器，开展SAR影像点目标高精度三维绝对定位研究。此外，采用同轨SAR影像开展目标点三维定位时，较小的卫星空间距离交会角使得定位精度受观测误差影响严重，降低了定位解算的稳定性。鉴于此，本文给出了一种添加高程限制条件的观测方程，成功解决了目标点定位解算不稳定问题。下面首先给出SAR影像点目标绝对定位原理和相关数据处理流程，然后给出基于TerraSAR-X聚束模式SAR影像和Sentinel-1A/B干涉宽幅模式SAR影像的角反射器目标三维绝对定位实验结果，最后对定位结果进行分析和总结。

1 SAR影像点目标绝对定位

1.1　基本原理

SAR影像点目标绝对定位的基本原理是采用多景SAR影像的斜距观测值和卫星精密轨道信息，通过空间距离‒多普勒交会的方法实现目标精确定位。根据R-D模型，目标定位的观测方程如下：

$$\begin{array}{l}\left|X_{\mathrm{S}}(t_{\mathrm{A}})-X_{\mathrm{T}}\right|-c\cdot {\tau }_{\mathrm{\rho }}/\mathrm{2}=\mathrm{0}\\ \frac{V_{\mathrm{S}}(t_{\mathrm{A}})(X_{\mathrm{T}}-X_{\mathrm{S}}(t_{\mathrm{A}}))}{\left|V_{\mathrm{S}}(t_{\mathrm{A}})\right|\left|X_{\mathrm{T}}-X_{\mathrm{S}}(t_{\mathrm{A}})\right|}=\mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}$$

（1）

式中：$X_{\mathrm{S}}(t_{\mathrm{A}})$表示$t_{\mathrm{A}}$时刻卫星所在的位置；$X_{\mathrm{T}}$为待求的目标点三维坐标；$c$表示光速；${\tau }_{\mathrm{\rho }}$表示雷达发射电磁波到目标点经历的双程时间，即SAR成像的距离向时间；$X_{\mathrm{S}}$和$V_{\mathrm{S}}$分别表示卫星的位置和速度，它们的值可利用已知的卫星精密轨道参数由卫星拍摄目标点的时刻，即方位向时间$t_{\mathrm{A}}$来计算；$\alpha$为卫星与目标点连线与卫星轨道垂向的夹角，对于TerraSAR-X和Sentinel-1A/B系统，它们在成像时已完成了零多普勒几何转换，此时$\alpha =\mathrm{0}$。图1是SAR影像目标三维定位的示意图。图中，S1~S4表示卫星。显然这种定位方式的精度和稳定性与空中SAR卫星之间距离（即空间基线的长度）有关。

图1 SAR影像目标绝对定位几何示意图

Fig. 1 Typical geometry of SAR absolute geolocation

因此，理论上只要采用两景或两景以上的SAR影像来联立R-D方程，线性化后通过最小二乘法即可获取目标点的三维坐标。假设输入$n$景SAR影像，则对于待定位的目标像元点，其观测值向量可表示为：$l=[{\tau }_{\mathrm{\rho }}^{\mathrm{1}},t_{\mathrm{A}}^{\mathrm{1}},\dots ,{\tau }_{\mathrm{\rho }}^n,t_{\mathrm{A}}^n]$。对式（1）进行线性化，误差方程可简化表示为

$$\mathit{B}\mathit{v}+\mathit{A}\mathit{x}+\mathit{w}=\mathrm{0}$$

（2）

式中：$\mathit{x}$为待估参数，即目标点的坐标改正值：$\mathit{x}=[\mathrm{d}{x}_{\mathrm{T}},\mathrm{d}{y}_{\mathrm{T}},\mathrm{d}{z}_{\mathrm{T}}]$。由于观测方程是相互独立的，设计矩阵$\mathit{B}$可逆，式（2）可转换为

$$\mathit{v}=\widehat{\mathit{B}}\mathit{x}-\mathit{L}$$

（3）

式中：$\widehat{\mathit{B}}=-{\mathit{B}}^{-\mathrm{1}}\mathit{A}$，$\mathit{L}={\mathit{B}}^{-\mathrm{1}}\mathit{w}$。上述问题转换为经典的参数平差问题，利用最小二乘法解算点目标的最优估值：

$$\widehat{\mathit{x}}={\mathit{N}}^{-\mathrm{1}}{\widehat{\mathit{B}}}^{\mathrm{T}}\mathit{P}\mathit{L}$$

（4a）

$$\mathit{N}={\widehat{\mathit{B}}}^{\mathrm{T}}\mathit{P}\widehat{\mathit{B}}$$

（4b）

方位向时间和距离向时间分别与卫星飞行速度和光速相联系，它们在定位解算中的权重不同。由于无法精确已知两类观测值的先验中误差，故不能给出较为合理的权比，因此采用方差分量估计来确定其权比，步骤如下：

（1） 首先根据SAR卫星的斜距精度和影像分辨率，给定距离向时间和方位向时间观测值较为合理的先验中误差${\sigma }_{\mathrm{\rho }}$和${\sigma }_{\mathrm{A}}$，进行第一次定权，即：

$$p_{\mathrm{\rho }}^{(\mathrm{1})}=\frac{{\sigma }_{\mathrm{\rho }}^{\mathrm{2}}}{{\sigma }_{\mathrm{\rho }}^{\mathrm{2}}}=\mathrm{1},\hspace{1em}{p}_{\mathrm{A}}^{(\mathrm{1})}=\frac{{\sigma }_{\mathrm{\rho }}^{\mathrm{2}}}{{\sigma }_{\mathrm{A}}^{\mathrm{2}}}$$

（5）

（2） 第一次解算后，利用残差估计出两类观测值的单位权方差${\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}$和${\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{A}}}^{\mathrm{2}}$：

$$\begin{array}{l}{\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}=\frac{{\mathit{v}}_{\mathrm{\rho }}^{\mathrm{T}}{\mathit{P}}_{\mathrm{\rho }}{\mathit{v}}_{\mathrm{\rho }}}{n_{\mathrm{\rho }}-\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}({\mathit{N}}^{-\mathrm{1}}{\widehat{\mathit{B}}}_{\mathrm{\rho }}^{\mathrm{T}}{\mathit{P}}_{\mathrm{\rho }}{\widehat{\mathit{B}}}_{\mathrm{\rho }})}\\ {\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{A}}}^{\mathrm{2}}=\frac{{\mathbf{v}}_{\mathrm{A}}^T{\mathbf{P}}_{\mathrm{A}}{\mathbf{v}}_{\mathrm{A}}}{n_{\mathrm{A}}-\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}({\mathit{N}}^{-\mathrm{1}}{\widehat{\mathit{B}}}_{\mathrm{A}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{A}}{\widehat{\mathit{B}}}_{\mathrm{A}})}\end{array}$$

（6）

（3） 重新定权：

$$\begin{array}{l}{p}_{\mathrm{\rho }}^{(k)}=\frac{{\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}}{{\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}/p_{\mathrm{\rho }}^{(k-\mathrm{1})}}=\mathrm{1}\\ p_{\mathrm{A}}^{(k)}=\frac{{\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}}{{\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{A}}}^{\mathrm{2}}/p_{\mathrm{A}}^{(k-\mathrm{1})}},k=\mathrm{2,3},\mathrm{4},\dots \end{array}$$

（7）

（4） 反复进行步骤（2）和步骤（3），直至${\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{\rho }}}^{\mathrm{2}}$、${\widehat{\sigma }}_{{\mathrm{0}}_{\mathrm{A}}}^{\mathrm{2}}$两者之差小于预定阈值为止。

为了保证空间距离交会的几何稳定性，定位计算应尽可能采用具有较大空间基线的SAR影像集，通常优先考虑采用升降轨、不同平台的SAR影像。然而，传统的三角形角反射器只能用于选定的SAR卫星轨道，难以在升、降轨影像中均被成功观测到。当采用同轨SAR影像时，为了保证干涉测量较好的空间相干性，其轨道空间基线通常只有几百米，较小的空间基线导致空间距离交会角度过小，定位解算在高程方向的精度较差，偏差可达到几米乃至几十米，同时也影响了距离向和方位向上的定位精度。为此，本文提出通过将目标点的高程固定为一个较为合理的值，提高定位解算的稳定性和精度。实际应用中，可采用LiDAR点云数据获取观测区域的DSM，利用空间分布特征将目标点与DSM进行匹配，即可获得目标点较为合理的高程信息，此时在定位模型中加入如下对应高程不变限制条件的观测方程：

$$\mathrm{\Delta }H=(\mathrm{c}\mathrm{o}\mathrm{s}B\mathrm{c}\mathrm{o}\mathrm{s}L,\mathrm{c}\mathrm{o}\mathrm{s}B\mathrm{s}\mathrm{i}\mathrm{n}L,\mathrm{s}\mathrm{i}\mathrm{n}B)\cdot \mathit{x}=\mathrm{0}$$

（8）

式中：$B$和$L$表示目标点初始纬度和经度，一般用SAR影像地理编码获得的经纬度作为近似值代入，在迭代求解中不断更新。联立此方程后，目标点的高程在定位迭代解算中保持不变，提高了定位解算的稳定性。

1.2　角反射器精细点目标分析

SAR影像绝对定位的方位向时间和距离向时间观测值是通过目标点的像素坐标间接计算得到的：

$$\begin{array}{l}{t}_{\mathrm{A}}=t_{\mathrm{0}}+(x_{\mathrm{p}\mathrm{t}}-\mathrm{1})\cdot \mathrm{\Delta }t\\ {\tau }_{\mathrm{\rho }}={\tau }_{\mathrm{0}}+(y_{\mathrm{p}\mathrm{t}}-\mathrm{1})\cdot \mathrm{\Delta }\tau \end{array}$$

（9）

式中：$(x_{\mathrm{p}\mathrm{t}},y_{\mathrm{p}\mathrm{t}})$为目标点在SAR影像中的方位向和距离向像素坐标，需要通过点目标精细分析获得；${\tau }_{\mathrm{0}}$、$t_{\mathrm{0}}$分别表示距离向起始像元的往返距离时间和方位向起始像元的UTC（universal time coordinated）时刻；$\mathrm{\Delta }\tau$、$\mathrm{\Delta }t$分别为像素点间的距离向和方位向采样时间间隔，这些参数可从SAR影像的产品描述文件中直接读出。因此要实现SAR影像目标精确三维定位，首先要解决两个问题：一是不同SAR影像中同名点识别，二是要确定同名点在各自影像中的精确像素位置。安装在地面上的角反射器后向散射信号强度高，在SAR强度影像中表现为一个具有较高信噪比的亮点目标，可较为容易地被识别出来。另一方面，在SAR影像上，角反射器通常呈现十字形亮点，十字中心的后向散射能量最强，即为角反射器的相位中心所在像素位置。但由于SAR影像的一个像素点代表其对应分辨率下所有地物的后向散射总贡献，因此要选取角反射器产生的信号窗口进行精细点目标分析，获取能够满足其高精度绝对定位的理论上的亚像素坐标位置。

首先利用角反射器的大致地理坐标位置根据R-D模型反算出其对应的影像坐标，寻找附近区域的强度最大值位置，通过人工目视判读，确定该位置是否是角反射器所在的像素位置；然后根据强度分布特性，在SAR强度影像中截取包含角反射器散射信号的$n\times n$窗口的子块；接下来采用双线性插值函数或SINC函数对截取子块的距离向和方位向同时进行$m$倍过采样；然后计算过采样影像每一行、每一列的质心坐标和平均强度值^［16］：

$$\begin{array}{cc}{\overline{x}}_i=\frac{{\displaystyle \sum _{k=\mathrm{1}}^n}\left(\frac{k}{m}\cdot I_k\right)+{\displaystyle \sum _{j=\mathrm{1}}^{m\times n}}\left[j\cdot f(i,j)\right]}{{\displaystyle \sum _{k=\mathrm{1}}^n}{I}_k+{\displaystyle \sum _{j=\mathrm{1}}^{m\times n}}f(i,j)}& \overline{f_i}=\frac{{\displaystyle \sum _{j=\mathrm{1}}^{m\times n}}f(i,j)}{m}\\ {\overline{y}}_j=\frac{{\displaystyle \sum _{k=\mathrm{1}}^n}\left(\frac{k}{m}\cdot I_k\right)+{\displaystyle \sum _{i=\mathrm{1}}^{m\times n}}\left[i\cdot f(i,j)\right]}{{\displaystyle \sum _{k=\mathrm{1}}^n}{I}_k+{\displaystyle \sum _{i=\mathrm{1}}^{m\times n}}f(i,j)}& \overline{f_j}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{m\times n}}f(i,j)}{m}\end{array}$$

（10）

式中：$I_k$为被插值像素在原始影像中所处行、列的强度值；$f(i,j)$为插值后像素的强度值；$({\overline{x}}_i,{\overline{y}}_i)$和$(\overline{f_i},\overline{f_j})$分别为过采样影像第$i$行和第$j$列的质心坐标以及对应的平均强度值。最后计算过采样影像的质心坐标$(x_{\mathrm{c}},y_{\mathrm{c}})$：

$$\begin{array}{l}{x}_{\mathrm{c}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^m}{\overline{x}}_i\overline{f_i}}{{\displaystyle \sum _{i=\mathrm{1}}^m}\overline{f_i}}\\ y_{\mathrm{c}}=\frac{{\displaystyle \sum _{j=\mathrm{1}}^m}{\overline{y}}_j\overline{f_j}}{{\displaystyle \sum _{j=\mathrm{1}}^m}\overline{f_j}}\end{array}$$

（11）

根据截取子块起始像素在原始SAR影像中的像素坐标，可最终得到角反射器精确的像素坐标，并代入式（9），即可得到方位向时间和距离向时间观测值。

1.3　观测值改正

为了获取高精度的绝对定位结果，需要对方位向时间和距离向时间观测值进行误差改正。

（1） 地表运动改正

板块构造运动、固体潮、海潮等引起的大范围地表形变会对SAR影像高精度定位产生干扰。固体潮是由于日月引力引起的固体地球表面变形，在一天内引起的周期性形变可达到几个分米。对于固体潮改正，首先采用Milbert提供的程序（https：//geodesyworld.github.io/SOFTS/solid.htm）计算其在卫星拍摄日期固体潮引起的北、东、上方向的位移，时间采样为60 s，通过插值可得SAR影像拍摄待定位目标时刻的位移${\delta }_{\mathrm{N}}$、${\delta }_{\mathrm{E}}$和${\delta }_{\mathrm{U}}$，并投影到方位向和斜距方向：

$$\begin{array}{l}{\delta }_{\mathrm{A}\mathrm{z}}=-{\delta }_{\mathrm{E}}\mathrm{s}\mathrm{i}\mathrm{n}\beta -{\delta }_{\mathrm{N}}\mathrm{c}\mathrm{o}\mathrm{s}\beta \\ {\delta }_{\mathrm{R}\mathrm{g}}=-{\delta }_{\mathrm{U}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{l}\mathrm{o}\mathrm{c}}+{\delta }_{\mathrm{E}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{l}\mathrm{o}\mathrm{c}}\mathrm{c}\mathrm{o}\mathrm{s}\beta -{\delta }_{\mathrm{N}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{l}\mathrm{o}\mathrm{c}}\mathrm{s}\mathrm{i}\mathrm{n}\beta \end{array}$$

（12）

式中：${\theta }_{\mathrm{l}\mathrm{o}\mathrm{c}}$为SAR成像的局部入射角；$\beta$表示北方向顺时针旋转到卫星飞行方向的角度。板块构造运动、海潮、极潮、大气负荷等对陆地产生的影响相对较小，量级一般在几个厘米以内，本文在定位解算中忽略了这些因素的干扰。

（2） 大气延迟改正

SAR卫星到目标点的斜距是通过距离向时间乘以真空光速计算得到的，信号传播路径的大气延迟会严重影响定位结果的精度。大气延迟分为对流层延迟和电离层延迟，其中对流层延迟与地形高程和入射角有关，引起的斜距误差高达数米，而电离层延迟的影响相对较小，一般给斜距观测值引入厘米级至分米级的误差。对流层延迟可采用连续GNSS观测站获取的对流层天顶延迟（zenith path delay，ZPD）进行改正：

$${\delta }_{\mathrm{t}\mathrm{r}\mathrm{o}}=\frac{I_{\mathrm{Z}\mathrm{P}\mathrm{D}}}{\mathrm{c}\mathrm{o}\mathrm{s}\theta }{\mathrm{e}}^{-(\frac{h_{\mathrm{S}\mathrm{A}\mathrm{R}}-h_{\mathrm{G}\mathrm{P}\mathrm{S}}}{h_{\mathrm{0}}})}$$

（13）

式中：$h_{\mathrm{0}}=\mathrm{6}\mathrm{000}\mathrm{m}$；$h_{\mathrm{S}\mathrm{A}\mathrm{R}}$和$h_{\mathrm{G}\mathrm{P}\mathrm{S}}$分别为SAR影像地面目标点高程和GPS测站高程。该方法一般要求GNSS（全球导航卫星系统）观测站要位于目标点附近50 km以内的范围^［11］。如果没有合适的GNSS观测数据，可根据气压、温度、湿度数据采用大气模型进行该项改正。电离层延迟改正则采用全球竖向总电子含量（ertical total electron content，VTEC）格网数据进行估计：

$${\delta }_{i\mathrm{o}\mathrm{n}}=K\cdot \frac{I_{\mathrm{V}\mathrm{T}\mathrm{E}\mathrm{C}}}{f^{\mathrm{2}}\cdot \mathrm{c}\mathrm{o}\mathrm{s}\theta }$$

（14）

式中：$f$表示载波频率；$K=\mathrm{40.28}{\mathrm{m}}^{\mathrm{3}}\cdot {\mathrm{s}}^{-\mathrm{2}}$。由于电离层大致分布在大气层100 ~ 1 000 km的范围内，而SAR卫星的轨道高度约为500 km，因此由式（14）计算得到的电离层延迟需要一定程度的减小，已有学者通过回归分析得出实际电离层延迟值约为按式（14）计算值的75 %^［17］。

（3） SAR成像改正

TerraSAR-X和Sentinel-1系统分别使用各自的TMSP（TerraSAR-X multimode SAR processor）和IPF（instrument processing facility）程序对卫星采集的雷达回波信息进行成像处理，将原始数据组织为二维（方位向和距离向）的单视复数（single look complex，SLC）影像，并同时将每个像元的成像几何转换为零多普勒几何。然而，最终的SLC影像中仍然存在由于成像处理带来的方位向时间和距离向时间偏差。其中方位向时间偏差主要是由SAR影像成像中广泛使用的“停‒走‒停”假设引起的，即假设卫星对地面发射雷达脉冲后停留在原地，等接收到地物回波信息后再继续向前飞行。但实际上，SAR卫星在不断地向前飞行并同时发射和接收雷达信号。TerraSAR-X卫星的TMSP程序未采用“停‒走‒停”假设，但是Sentinel-1卫星的IPF程序采用了该假设，因此需要对Sentinel-1 SAR影像的方位向时间进行零多普勒偏移改正：

$$\mathrm{\Delta }{t}_{\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{o}-\mathrm{D}\mathrm{o}\mathrm{p}}={\tau }_{\mathrm{m}\mathrm{i}\mathrm{d}}/\mathrm{2}+{\tau }_{\mathrm{\rho }}/\mathrm{2}-{\tau }_{\mathrm{0}}$$

（15）

式中：${\tau }_{\mathrm{m}\mathrm{i}\mathrm{d}}$是一个常数，它等于Sentinel-1 TOPS SAR影像第二个子条带中心像元的距离向时间，可从Sentinel-1影像的产品描述文件中读出。

SAR卫星运动同时会导致成像处理在对距离向脉冲聚焦时产生多普勒频移，因此需要对距离向时间进行多普勒频移改正：

$$\mathrm{\Delta }{\tau }_{\mathrm{D}\mathrm{o}\mathrm{p}-\mathrm{s}\mathrm{h}\mathrm{i}}=f_{\mathrm{D}\mathrm{C}}/K_{\mathrm{r}}$$

（16）

式中：$K_{\mathrm{r}}$表示距离向信号的调制速率，可从影像产品描述文件中获取；$f_{\mathrm{D}\mathrm{C}}$表示多普勒中心频率。此外，对于TOPS模式SAR影像，模拟的方位向多普勒调频斜率与真实值之间的差异会给方位向时间引入与地形相关的偏差。TerraSAR-X在计算方位向调频斜率时不断更新地形高程，而Sentinel-1只是采用平均高程进行计算，因此要对Sentinel-1 SAR影像的方位向时间进行相应的改正：

$$\mathrm{\Delta }{t}_{\mathrm{F}\mathrm{M}-\mathrm{m}\mathrm{i}\mathrm{s}}=f_{\mathrm{D}\mathrm{C}}\cdot (\frac{\mathrm{1}}{K_{\mathrm{g}\mathrm{e}\mathrm{o}}}-\frac{\mathrm{1}}{K_{\mathrm{a}}})$$

（17）

式中：$K_{\mathrm{a}}$表示多普勒调频斜率，具体的计算方法可参考ESA发布的Sentinel-1 IPF技术文档（https：//sentinel.esa.int/web/sentinel/user-guides/document-library），$K_{\mathrm{g}\mathrm{e}\mathrm{o}}$表示真实的方位向多普勒调频斜率：

$$K_{\mathrm{g}\mathrm{e}\mathrm{o}}=\frac{\mathrm{2}\cdot [({\mathit{X}}_{\mathrm{S}}-{\mathit{X}}_{\mathrm{T}})\cdot {\mathit{A}}_{\mathrm{S}}^{\mathrm{T}}+{\mathit{V}}_{\mathrm{S}}\cdot {\mathit{V}}_{\mathrm{S}}^{\mathrm{T}}]}{\lambda \cdot \left|{\mathit{X}}_{\mathrm{S}}-{\mathit{X}}_{\mathrm{T}}\right|}$$

（18）

式中：$A_{\mathrm{S}}$为卫星的加速度；$\lambda$为卫星载波波长。

2 实验分析

SAR影像目标绝对定位实验采用了6个地面三角形角反射器，其中，有3个角反射器布设于武汉某特高压试验基地，直角边长为0.8 m，分别编号WH1、WH2和WH3，如图2a所示。另外3个角反射器布设于上海市，直角边长为1.2 m，其中1个布设于上海市崇明岛东滩湿地公园内，编号SH1，如图2b所示，2个布设于上海长江大桥上，分别编号SH2和SH3，如图2c所示。考虑到野外防风和防积水，在上海角反射器上安装了电磁波可穿透的聚乙烯塑料材质盖板。为了达到对SAR卫星发射微波脉冲最佳的反射效果，角反射器安装的朝向垂直于卫星航向，并通过调整倾角使得角反射器的中心指向线对准雷达微波的入射方向。

图2 布设于武汉/上海的角反射器

Fig. 2 CRs located at Wuhan/Shanghai

实验使用了3景覆盖武汉角反射器的TerraSAR-X高分辨率聚束（high-resolution spotlight，HS）模式SAR影像，15景覆盖上海角反射器的Sentinel-1A/B干涉宽幅（interferometric wide swath，IW）模式SAR影像，基本参数信息如表1所示。同时，为了验证定位结果的可靠性，采用差分GPS（全球定位系统）技术测量了6个角反射器相位中心的大地坐标。

表1 实验使用的SAR影像基本信息

Tab. 1 Basic parameters of SAR images used in this study

测站	SAR卫星	拍摄日期	拍摄模式	轨道方向	距离向（方位向像素间隔） /m
WH1	TerraSAR-X	2008/09/26~2008/12/23	HS	降轨	0.45（0.85）
WH2	TerraSAR-X	2008/09/26~2008/12/23	HS	降轨	0.45（0.85）
WH3	TerraSAR-X	2008/09/26~2008/12/23	HS	降轨	0.45（0.85）
SH1	Sentinel-1A/B	2020/07/29~2020/10/27	IW	升轨	2.33（13.95）
SH2	Sentinel-1A/B	2020/09/09~2020/10/27	IW	升轨	2.33（13.95）
SH3	Sentinel-1A/B	2020/09/09~2020/10/27	IW	升轨	2.33（13.95）

在每幅SAR影像的强度图上能较为容易地识别到6个角反射器的所在像素位置，如图3所示。图4则分别展示了6个角反射器在平均强度影像中的局部放大图，图中每个像元的数字代表该像素的后向散射强度，单位为dB。可以看到每个角反射器的回波信号极强，且对其周边邻近的多个像素产生影响。

图3 TerraSAR-X/Sentinel-1 SAR影像强度图

Fig. 3 TerraSAR-X/Sentinel-1 SAR intensity map with six CRs highlighted in red box

图4 角反射器邻近区域后向散射能量分布

Fig. 4 Backscatter energy distribution around the adjacent area of each CR

采用1.2节所述的角反射器精细点目标分析方法，计算了6个角反射器在SAR影像中的质心坐标，并由式（9）得到距离向时间和方位向时间观测值。接下来，对观测值进行固体潮、大气延迟和卫星成像偏移改正。其中武汉和上海角反射器的电离层延迟改正都采用了NASA CDDIS提供的时间间隔为15 min的全球VTEC格网数据，经纬度间隔分别为2.5°和5.0°。采用最邻近的4个点进行二次插值，再对时间插值，即可获取角反射器所处位置在卫星拍摄时刻的VTEC值，并根据式（14）计算得到雷达视线向的电离层延迟。武汉角反射器的对流层改正则采用距离角反射器约17 km的IGS站（WUH200CHN）发布的ZPD数据进行改正，上海角反射器则采用距离角反射器约26 km的位于同济大学的连续GPS测站计算的ZPD数据完成对流层改正。通过给距离向时间和方位向时间的误差改正值分别乘以光速和卫星轨迹速度，将其单位换算为米，表2则给出了上述各项改正的最大、最小和平均值。可以看出，对流层延迟是主要的误差，可达到几米以上，而电离层延迟可达到几个分米。固体潮改正在距离向有几个厘米乃至分米，而在方位向上影响较小，只有几个毫米。Sentinel-1影像成像偏移改正中，方位向的零多普勒偏移改正可达到米级，距离向多普勒频移改正可达到几个分米，而由于上海地区地形较为平坦，方位向多普勒调频斜率差异引起的方位向改正量相对较小。此外，TerraSAR-X影像的产品描述文件提供了方位向时间由于相对多普勒效应和仪器计时误差引起的系统偏移值，为-1.331 m。

表2 距离向时间和方位向时间观测值改正量

Tab. 2 Range time and azimuth time observation corrections

SAR卫星	改正量	维度	最小值 /m	最大值 /m	平均值 /m
TerraSAR-X	固体潮	方位向	0.002 4	0.012 5	0.007 4
		距离向	0.049 2	0.113 3	0.082 4
	对流层	距离向	2.772 3	3.035 7	2.869 3
	电离层	距离向	0.025 1	0.026 4	0.025 8
Sentinel-1	固体潮	方位向	0.004 4	0.042 7	0.024 0
		距离向	0.073 5	0.028 6	0.056 0
	对流层	距离向	3.224 9	3.472 8	3.409 9
	电离层	距离向	0.177 8	0.229 0	0.197 7
	零多普勒偏移	方位向	1.327 2	1.710 6	1.507 8
	多普勒频移	距离向	0.676 1	0.786 4	0.730 1
	调频斜率差异	方位向	-0.021 2	-0.126 6	-0.073 4

为了评估点目标分析提取的角反射器质心坐标的精度，首先根据GPS观测得到的角反射器相位中心坐标和SAR影像提供的轨道状态向量，反算出角反射器理论上不受任何误差干扰时的SAR影像像素位置，将其作为角反射器在SAR影像中的精确参考位置。然后采用点目标分析获得角反射器的质心坐标，按式（9）计算对应的方位向和距离向时间，并且对其添加大气延迟、固体潮和卫星成像偏移改正，用改正后的值得到一组新的角反射器像素位置，并且计算相对其精确参考位置的偏移量。图5给出了每个角反射器在TerraSAR-X和Sentinel-1A/B SAR影像中的偏移量分布。可以看到，TerraSAR-X影像中的三个角反射器像素位置在距离向和方位向上均有分米级的偏移量，但它们的分布较为集中。Sentinel-1A/B影像中的角反射器像素位置在距离向上同样有分米级的偏移量，但在方位向上的偏移量可达到米级。此外，相比TerraSAR-X角反射器，Sentinel-1A/B角反射器像素位置偏移量的离散程度可达到分米级。

图5 SAR影像中角反射器的距离向和方位向像素位置偏移量

Fig. 5 Pixel position offsets in range and azimuth direction of CRs in SAR images

最后采用式（1）所示的R-D定位模型进行迭代解算实现每个角反射器的三维绝对定位。由于采用的是同轨SAR影像，TerraSAR-X影像的最大空间基线约600 m，Sentinel-1A/B影像的最大空间基线约400 m，而卫星到角反射器的斜距超过500 km，空间交会角度极小，一个很小的观测值误差即会引起很大的定位偏差，故在定位模型中加入式（8）所示的高程限制条件。表3则展示了6个角反射器最终的定位结果。其中$({\sigma }_{\mathrm{R}},{\sigma }_{\mathrm{A}})$是距离向时间和方位向时间观测值验后精度，$({\sigma }_{\mathrm{X}},{\sigma }_{\mathrm{Y}},{\sigma }_{\mathrm{Z}})$是点位内符合精度，$(N,E,U)$是解算坐标与GPS坐标在北、东、上方向的差值。由表3可以看出，武汉角反射器的观测值验后精度和角反射器的点位内符合精度均可达到厘米级，而上海角反射器可达到分米级。但是6个角反射器的定位结果与GPS定位结果有一些偏差，这主要是由于GPS是基于点目标观测的，而SAR影像分辨率限定了其最小地物分辨单元。可以看到Sentinel-1影像角反射器定位结果与GPS测量的差值更大，且Sentinel-1 TOPS SAR影像较低的方位向分辨率也体现在其北方向上具有较大的偏差。此外，武汉角反射器和上海布设在东滩湿地公园内的角反射器SH1，具有良好的观测环境，故采用GNSS静态连续测量，获取了角反射器顶点（即理论上的相位中心）厘米级精度的点位坐标。而对于大桥上的两个角反射器SH2和SH3，由于观测条件所限，采用Trimble RTK小碟获取了角反射器邻近位置的点位坐标，未获取角反射器顶点的坐标。同时角反射器在SAR影像中的实际相位中心会受到背景地物散射信号的干扰而产生偏移，武汉的3个角反射器所在的外界环境相近，而对于上海的3个角反射器，安装在野外环境的角反射器SH1相比安装在大桥上的角反射器SH2和SH3受到的背景干扰较小。基于以上两个因素，由表3可看出，相比SH2和SH3，SH1的定位结果与GPS测量结果的偏差较小。

表3 角反射器的绝对定位结果

Tab. 3 Absolute geolocation results of CRs

测站	${\sigma }_{\mathrm{R}}$ /m	${\sigma }_{\mathrm{A}}$ /m	${\sigma }_X$ /m	${\sigma }_Y$ /m	${\sigma }_Z$ /m	$N$ /m	$E$ /m	$U$ /m
WH1	0.035	0.054	0.033	0.018	0.022	0.145	0.718	0.165
WH2	0.038	0.047	0.035	0.017	0.020	0.232	0.773	0.130
WH3	0.046	0.048	0.043	0.019	0.021	0.421	1.169	0.134
SH1	0.147	0.866	0.098	0.079	0.174	-1.150	-0.008	0.006
SH2	0.182	1.003	0.148	0.118	0.257	-1.317	0.459	0.348
SH3	0.173	0.361	0.084	0.065	0.093	-1.113	0.496	0.170

3 结语

合成孔径雷达观测技术具有空间分辨率高、全球覆盖、观测周期短等独特优势，采用SAR影像点目标绝对定位，只要布设地面角反射器，甚至直接利用天然的角反射点，无需充电维护，即可获得厘米级定位精度，是对传统卫星大地测量技术的有效补充，也是大地测量学科发展的一个新的技术增长点。本文给出了SAR影像目标三维绝对定位的原理公式和数据处理流程，采用TerraSAR-X和Sentinel-1A/B SAR影像对布设在地面的6个角反射器进行了定位实验，结果表明TerraSAR-X卫星可以获得角反射器厘米级的内符合定位精度，Sentinel-1A/B卫星可以获得角反射器分米级的内符合定位精度，展示了TerraSAR-X和Sentinel-1A/B SAR影像可以作为一种新的大地测量数据源，实现角反射器的精确定位测量，非常适合长期野外形变监测，在山体滑坡、地震、火山、活动断层形变监测等方面具有很强的应用潜力。

作者贡献声明

作者贡献声明：

宋瑞庆：完成实验设计、数据处理、实验结果分析及初稿写作。

伍吉仓：相关文献资料的收集与分析，研究的构思，指导实验设计、数据分析、论文写作与修改。

马明雷：完成角反射器的设计、加工及实地安装，野外数据采集，定位结果分析解译。

张香：完成角反射器的设计、加工及实地安装，野外数据采集，定位结果分析解译。

李陶：参与论文中点目标分析方法研究。