2020, Vol. 40
2. 首都师范大学资源环境与旅游学院，北京市西三环北路105号，100048
近30 a来，InSAR技术成功应用于地形测绘、形变监测、地质灾害隐患早期识别等领域，但大气效应会严重制约InSAR技术的高精度应用[1]，尤其是在太阳活动高峰期或高温潮湿区域[2]。大气延迟相位源于2次SAR成像时大气状况差异引起的电磁波传播延迟误差，往往会覆盖地表真实形变信号。为削弱InSAR中的大气效应，国内外学者已提出多种方法，主要包括：1)逐对分析法。Massonnet等[3]研究Landers地震震后形变时首次提出逐对分析法。通过对比不同时间间隔的干涉结果，大气信号可被定性识别分离，但逐对分析法无法定量计算大气延迟。2)干涉图叠加法。干涉图叠加法[4]通过对同区域不同SAR影像生成的N幅独立干涉图进行叠加平均，以达到降低大气效应的目的，然而该方法需用较多SAR数据才能达到较好效果。3)永久散射体法。永久散射体法[5]利用离散且在时间上稳定的自然反射器的相位值来估计形变信号。根据大气相位在时间域上高频、空间域上低频的特性，通过时间域低通滤波、空间域高通滤波将大气噪声和形变相位分离，但该方法需要大量SAR影像观测数据(一般超过15景)。4)外部辅助数据大气校正法。常用的外部数据包含地面气象站观测资料、地基GNSS台网观测数据[6]、高分辨率数值气象模型(如ECMWF[7])、星载空间辐射计提供的近红外水汽产品[8](如MODIS、MERIS、FY-1C)。
随着水汽产品的日益丰富以及空间分辨率和精度的提高，利用全球尺度的水汽产品去除InSAR大气延迟相位的方法备受关注[9]，其优势在于仅用2景SAR影像和2景水汽产品即可削弱甚至去除干涉图中的大气误差。但传统空间辐射计校正法在原理上忽略了沿水平方向的大气扰动和水汽异质性，会降低校正精度。为提高空间辐射计校正法精度，本文将对流层中水汽垂直分层效应集成到大气延迟相位计算模型中，提出一种改进的重轨InSAR大气校正方法，并用实例验证改进方法的可行性。
1 大气对InSAR测量的影响雷达信号穿过大气层时，传播路径会发生延迟和路径弯曲，主要表现为对流层延迟和电离层延迟[10]。雷达卫星沿太阳同步轨道运行，几乎在同一当地时成像，可尽可能保证2次SAR成像时电离层延迟之差接近于0。研究表明，对流层中大气水汽是造成对流层延迟的最主要原因，尤其是对X、C等短波波段，对流层大气延迟效应更为显著。
对于重轨InSAR，2次回波信号的干涉相位可表示为：
| $ \begin{array}{l} \;\;\;\;\;\;\;{\varphi _{{\rm{int}}}} = {\varphi _1} - {\varphi _2} = \\ \frac{{4{\rm{ \mathsf{ π} }}}}{\lambda }\left( {{R_1} - {R_2}} \right) + \frac{{4{\rm{ \mathsf{ π} }}}}{\lambda }\left( {\Delta {R_1} - \Delta {R_2}} \right) \end{array} $ | (1) |
式中，λ为雷达信号波长，R1、R2为SAR卫星到地面目标的斜距，ΔR1、ΔR2为主辅影像成像时刻的大气传播延迟，
重轨InSAR大气延迟相位引入的地形和地表形变相位误差分别为：
| $ {{\sigma _h} = \frac{\lambda }{{4{\rm{ \mathsf{ π} }}{B_ \bot }}}R{\rm{sin}}\theta {\sigma _\varphi }} $ | (2) |
| $ {{\sigma _{\Delta r, 2}} = \frac{\lambda }{{4{\rm{ \mathsf{ π} }}}}{\sigma _\varphi }} $ | (3) |
式中，σφ为干涉图中大气延迟相位，σh为大气延迟引起的高程误差，B⊥为垂直基线长度，σΔr, 2为二轨法InSAR地表形变相位误差。
2 空间辐射计监测大气校正方法空间辐射计提供的近红外水汽产品目前主要有MODIS、MERIS、FY-1C可降水汽(PWV)产品。天顶湿延迟(ZWD)表示水汽引起的天顶方向路径延迟，与可降水汽的关系为：
| $ {{\rm{ZWD}} = \mathit{\Pi } \cdot {\rm{PWV}} = \frac{{{R_0}{\rho _w}\left( {\frac{{{k_3}}}{{{T_m}}} + k{' _2}} \right)}}{{{{10}^6}{M_v}}}{\rm{PWV}}} $ | (4) |
式中，Π为无量纲系数(取值为6.0~6.5)，ρw为液态水密度，R0和Mv分别为普适水汽常数和水汽摩尔质量，k′2和k3为大气反射率常数，Tm为对流层加权平均大气温度。
对于侧视成像雷达，对流层水汽引起的斜距相位延迟为：
| $ {\Delta {\varphi _{{\rm{atm}}}} = \frac{{4\pi }}{{\lambda {\rm{cos}}\theta }}\left( {{\rm{ZW}}{{\rm{D}}_{{t_1}}} - {\rm{ZW}}{{\rm{D}}_{{t_2}}}} \right) = \frac{{4{\rm{ \mathsf{ π} }}}}{\lambda }\frac{{{\rm{ZPDDM}}}}{{{\rm{cos}}\theta }}} $ | (5) |
式中，λ和θ分别为雷达信号波长和入射角，t1、t2为2次水汽场的观测时间，ZPDDM为天顶路径延迟差分。
由式(5)可知，传统空间辐射计校正法在将天顶路径延迟转换为斜距向延迟时，直接将天顶路径延迟差分投影到斜距方向，显然未考虑对流层水汽沿水平方向上异质性对路径延迟的影响。本文尝试将对流层中水汽垂直分层集成到空间辐射计校正模型中，以优化重轨InSAR监测中的大气延迟相位。
3 改进的空间辐射计校正方法图 1为InSAR干涉测量中天顶湿延迟与斜距湿延迟的几何关系。InSAR为一种相对测量技术，水汽的高时空变化会严重影响雷达信号，由于水平方向上大气扰动和水汽异质性，雷达视线方向上的斜距湿延迟(SWD)不等于天顶湿延迟除以入射角余弦值，而为地面像元A到像元B之间所有像元上斜距湿延迟之和。
|
图 1 天顶湿延迟与斜距湿延迟的几何关系 Fig. 1 Geometric relationship between zenith wet delays and slant wet delays |
对流层中水汽含量占大气总体积1%~4%，近99%的水汽聚集于底层大气层中，一般在近地面含量最高，随着海拔的增高水汽含量逐渐减小。根据大气折射率沿高程方向的空间分布特征[11]，大气可看成若干个垂直层，每一层的大气折射率相同。因此，对流层中水汽同样呈垂直分层状态：海拔高度在2 km以下的水汽约占大气层中总水汽含量的1/2，约3/4的水汽含量集中在海拔4 km以下，12 km(对流层平均厚度)以下水汽含量占全部水汽总量的99%[12-13]。换言之，对流层中水汽可分为3个垂直层：层1、层2、层3，海拔高度分别为0~2 km、2~4 km、4~12 km，水汽含量分别占大气水汽总含量的50%、25%、25%(图 1)，每层的水汽密度在空间上近似一致，在每层所有像元上则表现为具有相同的大气延迟贡献率。因此，本文采用线性加权方程来描述3个垂直水汽层上的综合大气延迟：
| $ {{n_1} = 2{\rm{tan}}\theta /R, {n_2} = 2{\rm{tan}}\theta /R, {n_3} = 8{\rm{tan}}\theta /R} $ | (6) |
式中，n1、n2、n3分别为3个垂直水汽层上的像元个数，R为可降水汽产品的空间分辨率。
现解释水汽层1上所有像元在天顶湿延迟差分图上的权重确定过程。如上所述，层1中水汽含量占大气水汽总含量50%，则水汽层1中大气延迟量占整个对流层水汽延迟总量50%。水汽层1共分布有n1个像元，每个像元的大气延迟贡献率相等，则层1上每个像元的大气延迟贡献率为水汽含量除以像元个数，即1/2n1。将水汽层1上每个像元的大气延迟贡献率作为天顶湿延迟差分图对应像元上的权重值，则层1上每个像元的权重值为1/2n1。同理，水汽层2和水汽层3上每个像元在天顶湿延迟差分图上的权重值分别为1/4n2和1/4n3。联合式(6)，求得3个水汽层上每个像元的权重值(wi)为：
| $ {w_i} = \left\{ \begin{array}{l} R/\left( {4{\rm{tan}}\theta } \right), \left( {i = 1, \cdots , {n_1}} \right)\\ R/\left( {8{\rm{tan}}\theta } \right), \left( {i = {n_1} + 1, \cdots , {n_1} + {n_2}} \right)\\ R/\left( {32{\rm{tan}}\theta } \right), \left( {i = {n_1} + {n_2} + 1, \cdots , n} \right) \end{array} \right. $ | (7) |
此时斜距向大气延迟相位计算模型可改进为：
| $ \begin{array}{l} \Delta \varphi {' _{{\rm{atm}}}} = \frac{{4\pi }}{{\lambda {\rm{cos}}\theta }}\left( {\sum\limits_{i = 1}^n {{\rm{ZPDD}}{{\rm{M}}_i} \cdot {w_i}} } \right), \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {i = 1, 2, \ldots , n} \right) \end{array} $ | (8) |
式中，ZPDDMi和wi分别为像元i天顶路径延迟差和权重值。
4 改进的大气校正方法验证案例随着城市化迅猛发展，城市用水量急剧增加，导致北京地区开始出现地面沉降现象并呈加速趋势[14]。本文以欧空局提供的Envisat ASAR数据和与ASAR同时相的MERIS近红外水汽产品[15]作为数据源，研究北京地区在2008-12~2009-03期间地面沉降的空间分布状况，并验证改进的大气延迟相位校正模型在InSAR差分干涉处理中去除大气延迟误差的可行性。图 2为ASAR影像覆盖范围。
|
图 2 ASAR影像覆盖范围 Fig. 2 ASAR image coverage |
如表 1所示，本文使用条带式降轨VV极化ASAR数据，时间基线为105 d，垂直基线为150 m。基于GAMMA软件进行差分干涉处理，去除精轨矢量文件中轨道误差相位和30 m空间分辨率的SRTM DEM数据中地形相位。为使SAR影像保持较高的空间分辨率和达到去除斑点噪声的目的，采用1:5的距离向、方位向视数比，同时采用最小费用流法进行相位解缠。
|
表 1 Envisat ASAR数据 Tab. 1 Envisat ASAR data |
图 3为从差分干涉相位中去除地形相位和轨道误差后的解缠相位图，由图可见，地面沉降漏斗清晰连续，但靠近山区、西南部分区域(红圈标识)存在明显的大气相位，需进一步去除。
|
图 3 解缠后的干涉相位图 Fig. 3 The unwrapped interferogram |
本文选用MERIS水汽产品去除干涉相位图中的大气相位误差。MERIS水汽数据与ASAR数据同时相，为300 m全分辨率的MER_FRS_2P可降水汽产品(表 2)。
|
|
表 2 MERIS数据集 Tab. 2 MERIS data |
MERIS水汽反演算法对云层极其敏感，当存在云层时，MERIS反演的水汽值远小于真实值，并接近于0，如图 4(a)中黑色区域(云污染像素)所示。MERIS可提供用于识别水汽产品中云污染像素的云掩膜数据，本文借助云掩膜数据剔除MERIS水汽图中的云污染像素，并对剔除后的空像元值进行空间内插。为避免MERIS数据二次插值的误差累积，直接用含有云污染像素的MERIS水汽产品计算天顶湿延迟差，再对天顶湿延迟差进行单次内插，获得云污染像素填补后的天顶湿延迟差(图 5)。
|
图 4 MERIS可降水汽产品分布 Fig. 4 The distribution of MERIS PWV |
|
图 5 天顶湿延迟差分图分布 Fig. 5 The distribution of ZPDDM |
式(7)和式(8)中，MERIS水汽产品的空间分辨率R为0.3 km，ASAR雷达信号入射角θ为22.854 5°，波长λ为5.63 cm，可计算出3个垂直水汽层上每个像元在天顶路径延迟差分图上的权重值wi。再联合云污染像素内插后的天顶路径延迟差分图(图 5(b))，即可精确计算出改进算法的大气延迟相位Δφ′atm。
4.3 地面沉降监测结果对比与精度验证如图 6所示，未进行大气校正、传统空间辐射计监测法和改进的大气校正法提取的3个地面沉降场的最大沉降量均为-4 cm左右，沉降场的空间分布也较相似。对比3个地面沉降图中大气延迟分布可发现，在未进行大气校正的形变图中，北京昌平、怀柔、密云和河北廊坊地区存在明显的区域性大气延迟相位误差。经改进的大气校正法校正后，廊坊地区的大气误差已基本去除，昌平地区的大气误差改善明显，黄色形变区域大幅度减少为未发生形变的蓝色区域，怀柔和密云地区的大气误差也在一定程度上得到改善。因此从沉降场的空间分布上看，改进的大气校正法能有效减小大气延迟相位误差，可提取更精确的地面沉降结果。
|
图 6 3种方法提取的地面沉降结果对比 Fig. 6 Comparison of ground subsidence field extracted by three methods |
为定量对比3个形变结果，选取陆态网络中6个GNSS台站在SAR时间基线内沿垂直方向的形变量作为基准，分别对图 6中3个地面沉降结果进行验证。InSAR监测结果可反映雷达视线方向形变，本案例中将LOS向形变投影至地表垂直方向进行验证(表 3)。
|
|
表 3 InSAR监测结果与GNSS监测结果对比 Tab. 3 Comparison of InSAR monitoring results with GNSS monitoring results |
如图 7所示，北京地区地面沉降空间分布差异性较大，呈现明显的不均匀特性。东城区、西城区、海淀区、石景山区、丰台区、大兴区未发生地面沉降，而朝阳区、通州区、顺义区出现较大面积的地面沉降，形成2个明显的沉降漏斗：朝阳-通州沉降漏斗和朝阳来广营沉降漏斗，最大沉降位于朝阳来广营，累计沉降量达到4.5 cm。InSAR与GNSS监测结果对比发现，未进行大气校正的形变结果与GNSS监测结果差值最大达到1.127 cm，均方根误差为0.685 cm；传统的空间辐射计大气校正方法提取的形变结果与GNSS监测结果差值最大达到0.975 cm，均方根误差为0.603 cm；而改进的方法提取的形变结果与GNSS监测结果差值最大为0.554 cm，最小仅为0.071 cm，均方根误差为0.388 cm，说明2种方法的观测结果一致性较强，改进的InSAR大气延迟相位校正方法相比于未校正和传统的空间辐射计校正法具有更高的监测精度，能提取更准确、真实的地表形变信息。
|
图 7 地面沉降结果验证 Fig. 7 Verification of ground subsidence results |
本文分析了大气效应对InSAR干涉测量的影响，在传统空间辐射计大气校正方法的基础上，将3个垂直水汽层上每个像元的大气延迟贡献率作为天顶湿延迟差分图的权重，进而求出每个地面像元的斜距湿延迟，提出一种改进的重轨InSAR大气延迟相位校正方法。利用覆盖北京地区的ASAR数据和同时相的MERIS近红外水汽产品验证改进方法的可行性，并用陆态网络GNSS台站监测结果验证改进方法的精度，得出以下结论：
1) 利用覆盖北京地区的ASAR数据和MERIS水汽产品，可成功削弱ASAR干涉图中的大气延迟相位，证明本文提出的方法能有效校正重轨InSAR监测中的大气误差。
2) InSAR监测结果与GNSS台站监测结果对比表明，改进的大气校正方法的均方根误差最小，为0.388 cm，表明改进的大气校正方法具有更高的地表形变监测精度。
3) 2008-12~2009-03北京地区出现不均匀地面沉降现象，朝阳和通州部分地区存在明显的沉降漏斗，最大沉降量达到4.5 cm。
| [1] |
Goldstein R M. Atmospheric Limitations to Repeat-Track Radar Interferometry[J]. Geophysical Research Letters, 1995, 22(18): 2517-2520 DOI:10.1029/95GL02475
(  0) |
| [2] |
游新兆, 王琪, 乔学军, 等. 大气折射对InSAR影响的定量分析[J]. 大地测量与地球动力学, 2003, 23(2): 81-87 (You Xinzhao, Wang Qi, Qiao Xuejun, et al. Quantitative Estimation of Effect of Atmospheric Refraction on InSAR[J]. Journal of Geodesy and Geodynamics, 2003, 23(2): 81-87)
( 0) |
| [3] |
Massonnet D, Feigl K, Rossi M, et al. Radar Interferometric Mapping of Deformation in the Year after the Landers Earthquake[J]. Nature, 1994, 369(6477): 227-230 DOI:10.1038/369227a0
( 0) |
| [4] |
Sandwell D T, Price E J. Phase Gradient Approach to Stacking Interferograms[J]. Journal of Geophysical Research:Solid Earth, 1998, 103(B12): 30183-30204 DOI:10.1029/1998JB900008
( 0) |
| [5] |
Ferretti A, Prati C, Rocca F. Permanent Scatters in SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8-20 DOI:10.1109/36.898661
( 0) |
| [6] |
Yu C, Li Z H, Penna N T. Interferometric Synthetic Aperture Radar Atmospheric Correction Using a GPS-Based Iterative Tropospheric Decomposition Model[J]. Remote Sensing of Environment, 2018, 204: 109-121 DOI:10.1016/j.rse.2017.10.038
( 0) |
| [7] |
Tang W, Liao M S, Yuan P. Atmospheric Correction in Time-Series SAR Interferometry for Land Surface Deformation Mapping——A Case Study of Taiyuan, China[J]. Advances in Space Research, 2016, 58(3): 310-325 DOI:10.1016/j.asr.2016.05.003
( 0) |
| [8] |
Zeng Q M, Li Y, Li X F. Correction of Tropospheric Water Vapour Effect on ASAR Interferogram Using Synchronous MERIS Data[C]. IEEE International Geoscience and Remote Sensing Symposium, Barcelona, 2007
( 0) |
| [9] |
Li Z H, Muller J P, Cross P, et al. Assessment of the Potential of MERIS Near-Infrared Water Vapour Products to Correct ASAR Interferometric Measurements[J]. International Journal of Remote Sensing, 2006, 27(2): 349-365
( 0) |
| [10] |
Hanssen R H. Radar Interferometry: Data Interpretation and Error Analysis[M]. Dordrecht: Kluwer Academic Publisher, 2001
( 0) |
| [11] |
朱邦彦.InSAR对流层延迟校正及其在地表沉降监测中的应用研究[D].武汉: 武汉大学, 2017 (Zhu Bangyan. Research on Tropospheric Delay Correction of SAR Interferometry and Its Application in Land Subsidence Monitoring[D]. Wuhan: Wuhan University, 2017)
( 0) |
| [12] |
Högström U. Review of Some Basic Characteristics of the Atmospheric Surface Layer[M]//Garratt J R, Taylor P A. Boundary-Layer Meteorology 25th Anniversary Volume, 1970-1995. Dordrecht: Springer, 1996
( 0) |
| [13] |
刘晓阳, 毛节泰, 李成才, 等. 大气水汽总量的垂直分解[J]. 高原气象, 2007, 26(3): 453-459 (Liu Xiaoyang, Mao Jietai, Li Chengcai, et al. Estimating Moisture Profile Using Integrated Water Vapor[J]. Plateau Meteorology, 2007, 26(3): 453-459)
( 0) |
| [14] |
Chen M, Tomás R, Li Z H, et al. Imaging Land Subsidence Induced by Groundwater Extraction in Beijing(China) Using Satellite Radar Interferometry[J]. Remote Sensing, 2016, 8(6): 468 DOI:10.3390/rs8060468
( 0) |
| [15] |
宋小刚, 李德仁, 单新建, 等. 基于GPS和MODIS的ENVISAT ASAR数据干涉测量中大气改正方法研究[J]. 地球物理学报, 2009, 52(6): 1457-1464 (Song Xiaogang, Li Deren, Shan Xinjian, et al. Correction of Atmospheric Effect in ASAR Interferogram Using GPS and MODIS Data[J]. Chinese Journal of Geophysics, 2009, 52(6): 1457-1464)
( 0) |
2. College of Resource Environment and Tourism, Capital Normal University, 105 North-Xisanhuan Road, Beijing 100048, China