基于合成孔径雷达干涉测量(InSAR)生成DEM实验研究

Experimental Study of Synthetic Aperture Radar Interferometry to Generate DEM

摘要: InSAR作为一种可全天候、全时段、大范围、快速地获取DEM的新兴技术手段，具备影响深远的应用前景；本文针对InSAR技术生成DEM过程中各种类型的误差，采取“多视处理 + 外部DEM辅助”的方法，实现在低空间分辨率的DEM辅助下生成复杂地形区域的数字高程模型，通过评价干涉获取的DEM精度，分析得出InSAR DEM在平坦地区的精度较高，在地形起伏较大的区域因雷达侧视成像机理导致使得相干程度降低，影响了DEM的精度。最后研究了高程精度与地形因素的关系，分析坡度坡向对DEM精度的影响。结果表明，精度随坡度的增加而减小，而坡向的改变对精度影响较小。

Abstract: In this paper, DEM data from the Shuttle Radar Topography Mission (SRTM) was used to assist the generation of InSAR DEM, which is used to remove some interference errors in the process and im-prove the accuracy of DEM generation. Taking the SRTM DEM as the standard, the height difference between two DEM was obtained, and the accuracy of DEM was evaluated. Through analysis, it is con-cluded that the precision of the InSAR DEM in flat area is higher. On the contrary, the precision of DEM in the region of the terrain fluction was proved lower because of decreasing relevance leaded by radar side-view mechanism. In the end, the relationship between elevation accuracy and terrain factors is studied, and the influence of slope and aspect on DEM accuracy is analyzed. The results show that accuracy decreases with the increase of slope, and slope direction changes have little ef-fect on the accuracy.

1. 引言

数字高程模型(Digital Elevation Model, DEM) [1] 作为地面高程信息的载体，描述了地表起伏的形态特征，是人们研究地形地貌的一种重要方法。在20世纪50年代首次提出就受到了科学界极为广泛的关注，被广泛应用于地质、水文、环境监测、灾害评估等多个领域。

合成孔径雷达干涉测量技术(Synthetic Aperture Radar Interferometry, InSAR) [2] 是一种先进的空间对地观测技术，以星载或机载合成孔径雷达(SAR)系统获取的复影像对为数据源，利用其中的相位信息可获取地面三维信息和探测地表形变信息。近年来，随着卫星雷达技术的不断创新，InSAR技术的迅猛发展，利用InSAR提取大规模、高精度的DEM [3] 已经成为雷达遥感领域一个新的研究方向 [4]。

本论文使用Sarmap公司研发的ENVI SARscape雷达图像处理软件将伊朗Bam地区的Envisat-1卫星ASAR雷达影像进行图像配准 [5]、基线估算等前期处理，由影像对生成干涉图，再去除平地相位、滤除相位噪声，对缠绕相位进行相位解缠，相位转化到高程最后得到DEM，然后以SRTM-3 DEM作为标准评价干涉处理得到的InSAR DEM高程精度及误差影响因素，并以坡度和坡向为例分析地形因子对其精度的影响。

2. InSAR技术原理

干涉测量技术来源于光学成像原理 [6]，通过分析和处理干涉条纹来获取目标的相关信息。InSAR技术将干涉测量原理与合成孔径雷达传感器相结合，通过两次获取的相位差求得传感器到地面的距离差，以此获得地面的地形信息。

2.1. 合成孔径雷达

合成孔径雷达(SAR)是一种主动式的微波传感器，拥有先进的数字处理能力，能够自动识别目标，远距离高分辨率成像。工作原理如图1所示，是利用传感器通过向地面发射微波，然后接收其包含有地物信息的散射回波信号，通过发射信号和接收回波之间的时间差测定距离。

其过程可以用雷达方程式定量表示，即：

$P_r=\frac{P_t{G}_t^2{\lambda }^2\sigma }{\left(4{\text{π}}^3\right)R^4}$ (1)

其中 $\lambda$ 是波长， $P_r$ 雷达接收功率， $P_t$ 为发射功率， $\sigma$ 为雷达接受方向伤目标的有效面积。

Figure 1. Principles of radar image

图1. 雷达成像原理

2.2. InSAR基本原理

在SAR的基础上发展起来的合成孔径雷达干涉测量(InSAR)是SAR成像技术与干涉测量技术相结合的成果，干涉测量的本质是由于两次回波信号的高相干性，同一区域的两幅SAR图像中的相位差可以得到两次传播的路径差。由此可从空间直接获取大范围、高精度的地表形变信息和高程信息，对地表进行长时间的监测 [7]。

InSAR技术的基本原理是通过机载或星载SAR系统对同一地区两次成像，得到两幅SAR影像对，再对这两幅SAR复影像数据进行相干处理，获取地形高程数据及形变信息，可用于DEM的生成和地表形变监测。

根据成像时间分类，InSAR可以分为单次轨道(single-pass)和重复轨道(repeat-pass)两种模式 [8]。单次轨道干涉是指在同一机载或星载平台上的两副天线同时对地观测，一起接收地面回波信号，其中一副发射信号；重复轨道干涉是指一副天线在相邻的重复轨道对地面同一点进行两次近平行的观测，获取地面点的复影像对。无论采取哪种模式，其基本干涉原理是一样的。

根据基线距与飞行方向之间的关系，InSAR又可以分为沿轨(along-track)和横跨轨道(across-track)两种 [9]。沿轨是指基线与飞行方向平行，主要出现在机载双天线模式中，应用于海洋动态监测与分析，能测量运动目标沿距离方向的速度。横跨轨道是指基线与卫星飞行方向垂直，是比较常用的模式，在机载和星载平台中都可以应用。

目前国际上常用的传感器都是利用重复、横跨轨道干涉处理进行测量的，本文以单天线重复轨道模式为例，阐述InSAR DEM建立的基本原理 [10]。

2.2.1. 几何公式

如图2所示，为InSAR基本原理示意图，A₁和A₂分别表示两天线的位置，它们之间的距离B为基线距。基线B和水平方向的夹角为α，H为站台高度。地面点P到天线A₁的斜距为ρ，到天线A₂的斜距为 $\rho +{\delta }_{\rho }$。θ是A₁的基准角，P点高程为h。

两斜距的差 ${\delta }_{\rho }$ 可以通过相应像素的绝对相位差 $\phi$ 估算：

$\phi =-\frac{4\text{π}}{\lambda }{\delta }_{\rho }$ (2)

其中λ为雷达信号波长。

Figure 2. Principles of InSAR technique

图2. InSAR技术基本原理

由图2的几何模型并利用余弦定理可得：

${\left(\rho +{\delta }_{\rho }\right)}^2={\rho }^2+B^{\text{2}}+\text{2}\rho B\sin \left(\alpha -\beta \right)$ (3)

则：

$\sin \left(\theta -\alpha \right)=\frac{{\left(\rho +{\delta }_{\rho }\right)}^2-{\rho }^2-B^2}{2\rho B}$ (4)

因为 $\rho \gg B$ 且 $\rho \gg {\delta }_{\rho }$，将式(2)代入(4)得：

$\sin \left(\theta -\alpha \right)=\frac{{\delta }_{\rho }}{B}=-\frac{\lambda \phi }{4\text{π}B}$ (5)

$\theta =\alpha -\arcsin \left(\frac{\lambda \phi }{4\text{π}B}\right)$ (6)

所以，只要确定天线高度H、相位差 $\phi$ 、基线距B、基线与水平方向夹角α和天线入射角θ，就可以计算出P点的高程值h：

$h=H-{\rho }_1\cos \theta$ (7)

式(6)和式(7)揭示了干涉相位差与高程h之间的数学关系，若已知几何关系和轨道参数等，就可以从相位差计算出地面点的高程值。

需要注意φ是解缠以后的相位，而干涉测量得到的相位值是真实干涉相位在 $[-\text{π},\text{π})$ 内的主值，即缠绕相位，必须要经过相位解缠 [11] 操作才能获得与地表形变相关联的真实干涉相位以及目标的高程信息。

另外，空间基线在垂直和平行于主影像视线方向可分为平行基线和垂直基线 $B_{\parallel }$ 和垂直基线 $B_{\perp }$，由图2可知：

$B_{\parallel }=B\sin \left(\theta -\alpha \right)$ (8)

$B_{\perp }=B\cos \left(\theta -\alpha \right)$ (9)

由于 $\rho \gg B$，可以用 ${\delta }_{\rho }$ 近似替代 $B_{\parallel }$ 从而计算得出 。这两个基线分量在相位转化为高程的过程中十分重要， $B_{\perp }$ 的长度还决定了干涉测量系统中相位对高程的敏感程度， $B_{\perp }$ 越长，敏感程度越高，反之敏感程度越低。

2.2.2. 基本流程

InSAR采用单视复数影像 [12] (SLC)生成DEM的基本流程可分为三个部分：

一是数据导入。包括读取不同格式的卫星数据、精密轨道数据，前置滤波等前期处理工作。

二是生成干涉图。包括SAR影像对的精确配准，基线估算，干涉图和相干系数图的生成，后置滤波和相位解缠等。

三是生成DEM。包括相位到高程的转换，地理编码和生成数字高程模型。

从两幅SLC到生成DEM的数据处理流程可以用图3表示：

Figure 3. InSAR flowchart

图3. InSAR处理流程图

3. 实验数据及其处理

本文使用Sarmap公司研发的ENVI SARscape软件 [13] 进行InSAR处理。SARscape是基于ENVI遥感图像处理软件的雷达图像处理软件，提供图形化的操作界面，具有专业的雷达图像处理和分析功能，能自动检测、选择参考图像及匹配图像，将处理过程简洁化。

用SARscape读取Envisat ASAR文件，获得单视复影像(SLC)文件。对SLC像对干涉处理，产生干涉图，再经相位解缠等步骤可生成高精度的DEM。

3.1. 实验数据

巴姆城(Bam)位于伊朗东南部平原地带，是伊朗最古老的城市之一，被称为“沙漠里的翡翠”。该地区气候干燥，地表植被匮乏，空间失相干很小，因此成为了InSAR实验研究的典型目标区。

本实验选取欧空局Envisat-1卫星分别于2003年6月11日和12月3日获取的2幅伊朗Bam地区SAR图像作为源数据，实现干涉测量生成DEM的过程，所用数据为伊朗BAM地区的ASAR数据。

其基本参数见表1：

Table 1. Basic parameter of experiment data

表1. 实验数据基本参数

3.2. 数据处理

通过基线估算可以看出干涉测量的精度随相干性的增加而提高。具体结果见表2：

Table 2. Result of baseline estimation

表2. 基线估算结果

基线估算结果显示，基线长度在临界基线范围内，多普勒中心差在临界范围内。本次实验所用的干涉像对符合InSAR技术质量要求。

在SARscape中的Interferogram generation模块，输入两景SLC数据，以Bam地区SRTM-3 DEM数据作为参考DEM提供干涉处理所需的地面投影参数和参考椭球信息，经过配准和多视输出干涉图、条纹图和强度图。数据处理变化过程如图4所示：

(a) 干涉条纹图 (b) 去平后的干涉图

Figure 4. Result of Multi-look processing

图4. 多视处理结果

经过干涉处理生成的干涉图还带有较大的相位误差，滤波能去除相位噪声，生成干涉的相干系数图和滤波后的主图像强图如图5所示：

(a) 滤波后干涉图 (b) 相关系数图

Figure 5. Result of adaptive filtering and coherence calculation

图5. 自适应滤波及相干性计算结果

设置相干性阈值为0.18。为了限制区域增长过程中的相位突变，相干性阈值一般设置为设置在0.15~0.2之间，避免相位解缠结果中不连续的区域产生“相位孤岛”。解缠后的相位图如图6所示：

Figure 6. Phase entanglement diagram

图6. 相位纠缠图

即使轨道非常精确，都是需要轨道控制点的，用来在重去平的一步进行相位偏移修正(如恒定相位去除)。

首先需要选取地面控制点生成控制点文件，选取的控制点必须要分布在整幅图像上，避免选择地形残差条纹区域。选好的地面控制点如图7中的红点，本文共选取控制点25个。

Figure 7. Reference point distribution diagram

图7. 控制点分布图

重去平后的干涉图和解缠图结果如图8所示：

(a) 重去平后的干涉图 (b) 重去平后的解缠图

Figure 8. Result of reflattening

图8. 重去平结果

轨道精炼参数如下所示，其中第一行RMS小于15，符合精度要求：

A-priori achievable RMS (m) = 2.6471041704

Orbital Refinement fitting:

Orbital shift in X direction (m) = −1.2700766459

Orbital shift in Y direction (m) = 0.0897130228

Orbital shift in Z direction (m) = −0.6478210078

Dependency of the shift in X direction, from the azimuth position (m/pixel) = −0.0001297462

Dependency of the shift in Y direction, from the azimuth position (m/pixel) = −0.0001526803

Dependency of the shift in Z direction, from the azimuth position (m/pixel) = 0.0004257573

Absolute phase offset (rad) = −64.9712445135

Root Mean Square error (m) = 4.5133342975

Mean difference POINTS HEIGHT after Orbital refinement (m) = 1.4948147195

Standard Deviation POINTS HEIGHT after Orbital refinement (m) = 7.8454153612

3.3. 实验结果及分析

干涉图处理结束后，开始DEM的计算。

由干涉图获取经过轨道矫正和相位解缠的非缠绕相位值，加上根据每个像素点的相干系数和轨道参数，可以计算出对应地面点的高程H。对图像上该像素点的像素坐标(x, y)进行地理编码，转换为大地坐标(L, B)，将SAR数据从斜距或地距投影转换为地理坐标投影，从而生成干涉DEM。

本文中投影信息为UTM [14]，椭球参数来自WGS84坐标系统 [15]，WGS坐标分区号为40，产品相干性阈值为0.2。对相位缺少地区进行内插处理，最后生成的DEM分辨率为25 m。生成的InSAR DEM文件、相干图像、精度图像和分辨率图像如图9所示：

(a) InSAR DEM (b) 相干图 (c) InSAR DEM分辨率图 (d) InSAR DEM精度图

Figure 9. Experiment result

图9. 实验结果

3.3.1. 以SRTM-3 DEM为标准的精度评估

航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission, SRTM)是美国国家地理空间情报局(NGA——前身为国家影像制图局，NIMA)投资约两亿美元与国家航空航天局(NASA)合作开发的一项国际航天测绘项目。

以巴姆地区SRTM-3 DEM，作为参考生成InSAR DEM，同时作为评价标准分析InSAR DEM的精度。相同坐标范围的干涉DEM与参考DEM逐像元对应相减得到高程差异值，生成高程差异图。

图10为实验所得InSAR DEM与剪裁后的SRTM-3 DEM对比图。

(a) InSAR DEM (b) SRTM-3 DEM

Figure 10. Comparison chart of experiment

图10. 实验对比图

通过对比可以看出InSAR DEM的分辨率与SRTM-3 DEM的分辨率相差不大，能清楚地显示出地形起伏的纹理特征。

将两幅DEM逐像素相减，获得高程差异。如图11所示，可以看出高差变化与高程变化相似，说明平坦地区高程精度较高，而地形起伏较大的地区高程精度较差。也就是说，高程差与地形特征有关，即高程精度受地形因子影响。

Figure 11. Elevation differences diagram

图11. 高程差异图

如图12所示，可以看出误差基本遵循正态分布，离散度较小，分布曲线陡峭。统计分析得出高差标准差为12.5 m。

Figure 12. Distribution diagram of geodetic difference

图12. 高差分布图

实验得到的InSAR DEM的高程与SRTM-3 DEM存在误差主要原因可能有：

1) 实验所用的星载SAR影像对的成像时间相隔将近六个月，导致了时间去相干，因此产生较大的干涉相位噪声，使得提取的DEM精度受到影响。

2) SRTM-3 DEM的生成时间是2002年，实验区卫星数据的获取时间是2003年，在这段时间间隔了可能存在地形的变化，产生基线去相干。

3) 采用精密轨道参数并不能完全消除卫星轨道误差，由此产生的附加相位，会影响干涉成果的精度。

由于SAR本身的特性，在山区或地表起伏较大的地区，DEM精度受到的影响更大。因此在起伏较大的区域，InSAR DEM的精度较低。所以在相对平坦的地区，利用InSAR技术生成DEM是比较好的选择；而在地形相对陡峭的山区，两幅雷达影像的相干性降低，干涉处理误差增大，从而影响了InSAR DEM的精度。

要进一步改善InSAR技术提取DEM的精度，可采取的措施有：

1) 提高SAR影像的质量；

2) 利用高精度的地面控制点文件对基线参数进行纠正；

3) 利用成像时的大气数据对大气延迟进行修正；

4) 进行有效的图像失真校正；

5) 对雷达影像对进行多视处理消除热噪声。

3.3.2. 地形分析

基于栅格的DEM计算方法是对栅格图像上的每一个像元进行计算。本文以InSAR DEM为基础，通过Envi软件中的地形分析提取DEM研究区域各个点的坡度坡向，并以此为例分析地形因子对DEM精度的影响。

如图13所示，栅格图形中每一个点的像素单元值为该点的坡度值。

Figure 13. Gradient diagram

图13. 坡度图

将栅格图形分类汇总，得出0˚~68˚的坡度统计图和0˚~360˚范围内的坡向统计图，如图14和图15所示。

Figure 14. Distribution diagram of slope

图14. 坡度分布图

Figure 15. Distribution diagram of aspect

图15. 坡向分布图

图16为不同坡度与DEM的精度关系，图17表明了坡向对DEM的精度的影响。

图16表明，随着坡度增加，DEM标准差逐渐变大，精度逐渐降低。坡度在80˚左右标准差达到最大值。由图17可以看出，DEM的精度与坡向变化没有明显关系。

由此可以得出如下结论：

1) InSAR DEM的精度可以满足一般要求，即InSAR技术提取DEM是一个有效的方法。

2) InSAR DEM在平坦地区精度较高，在地形起伏较大的区域精度比SRTM-3 DEM低，本次实验中两幅DEM的高差标准差为12.5米。

3) 随坡度的增加，DEM精度逐步减小，即坡度越大精度越小。

4) InSAR DEM的精度与坡向没有明显关系。

Figure 16. Relation graph of slope and standard deviation

图16. 坡度与标准差关系图

Figure 17. Relation graph of aspect and standard deviation

图17. 坡向与标准差关系图

4. 结束语

本文对星载InSAR技术生成DEM的原理进行了深入分析，实现了从ENVISAT卫星数据到生成数字高程模型的全过程，研究了从干涉图的生成到由干涉图反演高程的过程中，限制DEM精度的各类因素。

本文利用InSAR技术，以外部DEM数据辅助生成干涉图，在相位转换前去除干涉图中的系统误差，以伊朗Bam地区的ENVISAT卫星数据进行实验，提取出DEM，并与SRTM-3 DEM的精度进行了比较。研究结果表明，利用星载SAR提取DEM的分辨率与SRTM-3 DEM的分辨率相近。