ENVISAT和ALOS卫星SAR影像提取地表沉降的对比试验

Comparative Experiment of Land Subsidence Extraction from ENVISAT and ALOS Satellite SAR Images

摘要: 本文介绍双轨法差分干涉测量中的数学模型及计算公式，并以ENVISAT和ALOS两种卫星得到的影像数据分析得到某一地区的地表沉降值，分析二者影像的不同点以及各自数据的特点。进行差分干涉处理，并进一步提取出地表沉降量，并且展开分析。对比解译和分析ENVISAT ASAR数据和ALOS PALSAR数据提取的差分干涉图和形变值，并最终得出二者卫星在提取地表形变中可行性的结论。

Abstract: This paper introduces the mathematical model and calculation formula in the two-track differential interferometry, and analyzes the land subsidence value of a certain area by using the image data obtained from ENVISAT and ALOS satellites, and analyzes the differences between the two images and the characteristics of their respective data. Carry out differential interference processing, and further extract and analyze the surface settlement. Comparative interpretation and analysis of dif-ferential interferograms and deformation values extracted from ENVISAT ASAR data and ALOS PALSAR data, and ultimately a conclusion that the two satellites are feasible in extracting surface deformation is reached.

1. 引言

地面沉降是一种可由多种因素引起的地面标高缓慢降低的环境地质现象，严重时会成为一种地质灾害。可能为人类社会带来一系列严重危害，如适宜人类生存发展的安全地面高度的丧失，房屋、道路、地铁及桥梁等建筑物的损坏，潮灾、涝灾灾情的加重等。地面沉降不仅给社会发展带来重大的经济损失，也严重威胁着人类的生命和财产安全。地表沉降，主要是由于地下水和其他地下资源(比如石油，天然气，煤炭)的过分开采引起的。出现地面沉降的地区一般范围大，沉降过程缓慢，所以早期一般不易察觉，也不易引起人们的重视。它多发生在大中城市，对人们的生产、生活、交通等影响极大，造成的损失和危害也大，成为一种严重的环境地质问题，影响和制约着当地国民经济的可持续发展。InSAR (合成孔径雷达干涉测量)是近几年迅速发展起来的高新技术，它是使用雷达信号的相位信息提取地球表面三维信息。主要应用于测量地面点的高程及其动态变化的测量 [1] 。InSAR技术能全天候、全天时地获取大面积地面精确三维信息，空间分辨率高，基本不受气候条件的影响。干涉成像雷达遥感与GPS的区别在于它具有比GPS更高的形变观测精度、采样密度高(100 m之内)、空间延续性好和无需建立地面接收站。INSAR与其它测量形变位移仪器的区别在于那些测量仪器成本高，布点少，且难以大范围追踪地壳形变。相关文献做了大量研究并取得了丰硕的成果 [2] 。

本论文以ERS-1/2卫星获取某地区的雷达影像，藉由GAMMA干涉处理程序以合成孔径雷达干涉技术生成垂直位移图。用不同时期产生的精度地形计算地表沉降。并探讨干涉生成过程中，不同模式影像对于精度地形计算地表沉降之影响。

2. 二轨法差分误差特性分析

单视复影像slc1与slc3分别于形变前后获取，以slc1为主图像与slc3进行干涉。图1为图像的成像几何关系， ${\bar{S}}_3={\left[x_3,y_3,z_3\right]}^{\text{T}}$ 、 ${\bar{S}}_3={\left[x_3,y_3,z_3\right]}^{\text{T}}$ 分别为成像中心时刻卫星在地心笛卡儿坐标系下的位置， $B_{13}$ 为基线矢量， ${\alpha }_1$ 为雷达下视角， ${\beta }_1$ 为雷达斜视角， ${\beta }_3$ 与 ${\phi }_3$ 分别为基线 ${\bar{B}}_{13}$ 的方位角与高度角，O点为地心，T为目标点， $\bar{h}$ 为形变矢量 [3] 。

为了计算方便，通过配准和共轭相乘后得到的干涉相位可以表示为

${\varphi }_{13}\approx {\varphi }_{topo}+{\varphi }_{defo}$ (1)

其中，

Figure 1. The geometric relationship between slc1 and slc2 imaging

图1. Slc1与slc2成像几何关系

$\begin{array}{c}{\varphi }_{topo}=-\frac{4\text{π}}{\lambda }{B}_{13pr}=-\frac{4\text{π}}{\lambda }\langle {\bar{B}}_{13},{\bar{r}}_1^u\rangle \\ =-\frac{4\text{π}}{\lambda }{B}_{13}\left[\sin a_1\cos {\phi }_3\cos \left({\beta }_1+{\beta }_3\right)-\sin {\phi }_3\cos {\alpha }_1\right]\end{array}$ (2)

${\varphi }_{defo}=-\frac{4\text{π}}{\lambda }{h}_r=-\frac{4\text{π}}{\lambda }\langle {\bar{r}}_1^u,\bar{h}\rangle$ (3)

利用观测区已知的DEM (数字高程)数据确定出地形相位 ${\varphi }_{topo}$ ，并将其从干涉相位中去除掉，利用剩余的形变相位反演地表形变 [4] 。

当我们得到观测区的DEM数据后，就可以通过R-D定位方法得到T的位置，此时利用下式就可以得到地形相位。

${\varphi }_{topo}=-\frac{4\text{π}}{\lambda }\left(r_1-{\langle \bar{T}-{\bar{S}}_3,\bar{T}-{\bar{S}}_3\rangle }^{\frac{1}{2}}\right)$ (4)

其中 $r_1$ 由回波延迟读出， ${\bar{S}}_3$ 由卫星测量设备得到。

地形相位去除后，就可以从形变相位中反演出实现向形变 $h_r$ ：

$h_r=\frac{\lambda }{4\text{π}}\left({\varphi }_{13}-{\varphi }_{topo}\right)$ (5)

将地形相位公式(4)中的矢量 $\bar{T}-{\bar{S}}_3$ 记为 $r_{TS}$ ，则 $r_{TS}$ 代表目标点T到 $s_3$ 的距离，其数值可由(6)式计算

$r_{TS}=\sqrt{\langle {\bar{r}}_1-{\bar{B}}_{13},{\bar{r}}_1-{\bar{B}}_{13}\rangle }$ (6)

地形相位表示成标量形式为 ${\varphi }_{topo}=-\frac{4\text{π}}{\lambda }\left(r_1-r_{TS}\right)$ 。假设已知基线的三维分量 ${\varphi }_{topo}=-\frac{4\text{π}}{\lambda }\left(r_1-r_{TS}\right)$ ， ${\bar{r}}_1$ 有回波延迟得出，则 ${\bar{r}}_1$ 表示为：

${\bar{r}}_1=r_1\left(\sin {\alpha }_1\cos {\beta }_1,\sin {\alpha }_1\sin {\beta }_1,\cos {\alpha }_1\right)$ (7)

此时，可以将(6)写成如下形式：

$\begin{array}{c}{r}_{TS}=\sqrt{{\left(B_{13x}-r_1\sin {\alpha }_1\cos {\beta }_1\right)}^2+{\left(B_{13y}-r_1\sin {\alpha }_1\sin {\beta }_1\right)}^2+\left(B_{13z}-r_1\cos {\alpha }_1\right)}\\ =\sqrt{B_{13}^2+r_1^2-2r_1{B}_{13pr}}\end{array}$ (8)

其中角度 ${\alpha }_1$ 与T点的地心距H有关，设R为 $S_1$ 到地心O的距离，则有：

$\cos {\alpha }_1=\frac{R^2+r_1^2-H^2}{2R{r}_1}$ (9)

下面计算干涉相位 ${\varphi }_{13}$ 、高程H、基线三维值 ${\left(B_{13x},B_{13y},B_{13z}\right)}^{\text{T}}$ 误差对视线向形变的误差传播系数，首先，将视线向形变反演公式(5)对 ${\varphi }_{13}$ 求偏导得

$\frac{\partial h_r}{\partial {\varphi }_{13}}=\frac{\lambda }{4\text{π}}$ (10)

同理，求 $h_r$ 对其他各变量的偏导为

$\frac{\partial h_r}{\partial B_{13x}}=\frac{r_1\sin {\alpha }_1\cos {\beta }_1-B_{13x}}{r_{TS}}$ (11)

$\frac{\partial h_r}{\partial B_{13y}}=\frac{r_1\sin {\alpha }_1\sin {\beta }_1-B_{13y}}{r_{TS}}$ (12)

$\frac{\partial h_r}{\partial B_{13z}}=\frac{r_1\cos {\alpha }_1-B_{13z}}{r_{TS}}$ (13)

公式(10)、(11)、(12)、(13)为二轨法差分干涉处理中的几个主要的误差传播系数计算公式。

由上面的误差传播系数可以看出，在“二轨法”差分干涉处理中，干涉相位的误差传播系数仅仅与波长有关，当波长较短时干涉相位误差对视线向形变的影响会明显减小。所以，在相同精度要求下，利用ENVISAT ASAR卫星图像进行“二轨法”差分干涉处理比利用PALSAR卫星图像进行“二轨法”差分干涉处理对干涉相位误差的容忍程度要高得多，因为ENVISAT ASAR所搭载的合成孔径雷达工作在C波段，PALSAR所搭载的合成孔径雷达工作在L波段 [5] 。

3. ALOS卫星的数据处理

日本地球观测卫星计划主要包括2个系列：大气和海洋观测系列以及陆地观测系列。先进对地观测卫星ALOS是JERS-1与ADEOS的后继星，采用了先进的陆地观测技术，能够获取全球高分辨率陆地观测数据，主要应用目标为测绘、区域环境观测、灾害监测、资源调查等领域。ALOS卫星载有三个传感器：全色遥感立体测绘仪(PRISM)，主要用于数字高程测绘；先进可见光与近红外辐射计-2 (AVNIR-2)，用于精确陆地观测；相控阵型L波段合成孔径雷达(PALSAR)，用于全天时全天候陆地观测。ALOS卫星采用了高速大容量数据处理技术与卫星精确定位和姿态控制技术，其分辨率可达2.5米。作为一颗科研卫星，ALOS和同类卫星相比具有影像质量高，价格低廉等特点。

3.1. 数据处理

以ALOS卫星获取的某一地区2019年01月10日到2019年02月25日的数据为例。

1) 生成干涉图

生成干涉图如图2所示。

Figure 2. Comparison of generated interferograms

图2. 生成的干涉图对比

2) 去除地形影响

去除地形影响后如图3所示。

Figure 3. Comparison after removing the influence of terrain

图3. 去除地形影响后的对比

3) 滤波

滤波处理之后效果如图4所示。

Figure 4. Comparison after filtering

图4. 滤波处理后的对比

4) 相位解缠

相位解缠之后效果如图5所示。

Figure 5. Comparison after phase unwrapping

图5. 相位解缠后的对比

5) 二轨差分

二轨差分处理之后效果如图6所示。

Figure 6. Comparison after two-track differential processing

图6. 二轨差分处理后的对比

3.2. 提取剖面线的形变信息

由于时间短和SAR技术的测量精度不很高，本实验就提取此地区的某几条重要的线，利用提取这几条线的垂直形变信息，每个象元点的灰度值均表示形变量。依次累加的形变信息图中提取了一条线的形变信息，得到每次测量的形变量累加值，并依次按时间显示。红色的曲线表示从2019年1月10日至2019年2月25日被提取的线的每个象元点形变量；黄色的曲线表示从2019年2月25日至2020年1月13日被提取的线的每个象元点形变量；绿色的曲线表示从2020年1月13日至2020年2月28日被提取的线的每个象元点形变量；纵轴的单位为m。下沉量统计如表1所示，下沉示意图如图7所示。

Table 1. Settlement statistics of profile line

表1. 剖面线下沉量统计

从上面所得的数据以及图表可知剖面线在2019年1月10日到2019年2月25日下沉较小，最大为250 mm，平均月下沉量为50 mm/月。在2019年2月25日到2020年1月13日这段时间下沉较大，最大为500 mm，平均月下沉量为30 mm/月。在2020年1月13日到2020年2月28日下沉最大为90 mm，平均月下沉量为35 mm/月。因为坐标点640-826，622-822，582-826经过采煤区域，所以下沉量比较大，符合实际情况。

Figure 7. Schematic diagram of settlement of profile line

图7. 剖面线的下沉量示意图

4. ENVISAT数据处理

Envisat是由欧空局发射的一颗先进的极轨对地观测卫星，载有10种传感器，其中有先进的合成孔径雷达ASAR (Advanced Synthetic Aperture Radar)。ASAR工作在C波段，具有主动相控天线系统，5种成像模式，7种成像条带及交替极化成像功能。ENVIASAR数据可以很好的区分不同地物，而且正确反映了地物的后向散射特征，可以有效的应用于地物探测，地物识别，土地覆盖分类和农作物生长监测等领域。

4.1. 数据处理

以ENVISAT获取的某一地区2020年1月240日到2020年10月24日的数据为例生成干涉图。

1) 生成干涉图

生成的干涉图效果如图8所示。

Figure 8. Comparison of generated interferograms

图8. 生成的干涉图的对比

2) 去除地形影响

去除地形影响之后的效果如图9所示。

Figure 9. Comparison after removing the influence of terrain

图9. 去除地形影响后的对比

3) 滤波

滤波处理之后效果如图10所示。

Figure 10. Comparison after filtering

图10. 滤波处理后的对比

4) 相位解缠

相位解缠之后效果如图11所示。

Figure 11. Comparison after phase unwrapping

图11. 相位解缠后的对比

4.2. 结果分析

由于时间短和SAR技术的测量精度不很高，本实验就提取此地区的某几条重要的线，如图12所示，利用提取这几条线的垂直形变信息，每个象元点的灰度值均表示形变量。依次累加的形变信息图中提取了一条线的形变信息，得到每次测量的形变量累加值，并依次按时间显示。红色的曲线表示从2019年12月20日至2020年1月24日被提取的线的每个象元点形变量；黄色的曲线表示从2020年1月24日至2020年4月4日被提取的线的每个象元点形变量；蓝色的曲线表示从2020年4月4日至2020年5月9日被提取的线的每个象元点形变量；灰色的曲线表示从2020年5月9日至2020年7月18日被提取的线的每个象元点形变量；棕色的曲线表示从2020年7月18日至2020年8月22日被提取的线的每个象元点形变量，纵轴的单位为mm。剖面线下沉统计量如图12所示，剖面线的形变量如表2，图13所示。

Figure 12. Extracted section line

图12. 提取的剖面线

Table 2. Settlement statistics of section line

表2. 剖面线的下沉量统计

Figure 13. Schematic diagram of shape variables of section line

图13. 剖面线的形变量示意图

从上面所得的数据以及图表可知剖面线b从2019年12月20日至2020年1月24日的最大下沉量为30 mm，平均月下沉量为30 mm/月；从2020年1月24日至2020年4月4日最大下沉量为20 mm，平均月下沉量为10 mm/月；从2020年4月4日至2020年5月9日最大下沉量为20 mm，平均月下沉量为20 mm/月；从2020年5月9日至2020年7月18日最大下沉量为20 mm，平均月下沉量为10 mm/月；从2020年7月18日至2020年8月22日最大下沉量为20 mm，平均月下沉量为20 mm/月。

5. ALOS与ENVISAT数据提取地表沉降的对比分析

选取ALOS卫星影像上的从2020年1月13日到2020年2月28日拍摄00113-00228卫星影像和ENVISAT卫星影像上获得的从2020年1月24日到2020年4月4日所拍摄的00124-00404影像进行对比研究，因为二者在时间度上有所重合，所以可以用来进行对比研究。

5.1. ALOS影像与ENVISAT影像对采煤区的对比分析

对采煤区的对比分析。对比分析统计如表3所示。

Table 3. Comparative analysis of coal mining area

表3. 对采煤区的对比分析

可以看出ENVISAT卫星的灵敏度要比ALOS卫星高，因为ALOS卫星波长长，对小的地形地物监测不如ENVISAT卫星灵敏。ENVISAT卫星对小的地形地物的监测要精确。

5.2. 对采空区的地表下沉对比分析

在采空区中，选取两幅影像上采空区的一系列的坐标点进行对比分析，统计变化量如表4所示，变化量对比图如图14所示。

Table 4. Comparative analysis of goaf

表4. 对采空区的对比分析

Figure 14. Schematic diagram of shape variables of goaf

图14. 采空区的形变量示意图

5.3. 结合数据的分析以及误差来源分析

从示意图可以看出ALOS卫星影像数据和ENVISAT卫星影像数据得出的结果基本相同，符合实际情况。但是ENVISAT数据基本是差了35天，ALOS相差46天；ALOS波长长，小变形不灵敏，ENVISAT的C波段，波长短，小变形灵敏；ALOS分辨率高，为10米，后者25米，是这个分辨单元的强散射体的回波，两者的散射单元会存在区别。ENVISAT ASAR数据和ALOS PALSAR数据提取的视线向最大形变值分别180 mm左右和250 mm左右。PALSAR的Ｌ波段干涉测量在高植被覆盖和地形起伏较大地区具有一定应用潜力，与ENVISAT ASAR的Ｃ波段干涉对相比，Ｌ波段数据良好的相干性更有利于相位的获取和识别。

误差分析以及来源：因为ALOS卫星采用L波段，波长长，对小变形不灵敏，ENVISAT采用C波段，波长短，小变形灵敏。再者由于我研究的时间较短，没能进行进一步的研究，两幅图像未进行配准，选取的坐标点不是完全的重合，导致选点有误差，两幅对比影像上面的选点不是能完全相同，会有一定的误差。这都造成了ALOS卫星与ENVISAT卫星数据上的差异。

去除平地效应和地形效应后的干涉图，还要对轨道残差的影响进行剔除，本人在实验中未进行此项操作。不同时间实施的SAR对地观测所对应的大气水份含量和电离层状态可能是不同的，该附加的相位给最终的量测结果增加了不确定性。解决方案包括选择波长较长的SAR系统获取的数据，或者利用多次重复观测的干涉影像对结果取平均。

不同时间实施的SAR对地观测所对应的大气水份含量和电离层状态可能是不同的，该附加的相位给最终的量测结果增加了不确定性。解决方案包括选择波长较长的SAR系统获取的数据，或者利用多次重复观测的干涉影像对结果取平均。INSAR干涉处理过分依赖原始数据质量，对相干性差的区域，不能获得令人满意的干涉图，应进一步探求提高干涉质量的方法。

6. 结论

通过对InSAR技术的原理分析知道，InSAR技术可以实现具有全天候监测的能力，可覆盖大范围地区，且数据获取的速度较快；不仅可以提供用大多数方法难以监测到的地区的数据，而且还可以进行可与传统测量技术相比拟的高质量的地表变形测量；可提供时间跨度较大的SAR影像数据；花费比其他方法少的费用即可获取观测数据，这使得长期对地面沉降进行监测成为可能。ENVISATASAR数据的多极化、多入射角、多成像模式和高的间分辨率特点将会提供更加丰富的雷达遥感信息，将会在海洋科学、海洋石油工业、国土资源勘察、船只与尾迹探测、水运与海运管理和土地覆盖分类等领域得到有效和广泛的应用。ALOS卫星在农业、土地利用遥感制图、环境监测、灾害监测、林业和海洋方面的应用，ALOS卫星可以提供丰富的信息，可广泛用于科研、生产和生活领域。ALOS卫星还可以实现快速、准确、动态的数据监测。及时获取信息为有关部门的决策提供技术支撑。ENVISAT数据基本是差了35天，ALOS相差46天；ALOS波长长，小变形不灵敏，ENVISAT的C波段，波长短，小变形灵敏；alos分辨率高，10米，后者25米，是这个分辨单元的强散射体的回波，两者的散射单元会存在区别。在利用双轨差分监测地震形变场时，ENVISAT ASAR数据干涉对的垂直基线最好小于250 m，否则会因为垂直基线过长而导致失相关，从而不能用于地震形变场的监测。