地基干涉合成孔径雷达图像非线性大气相位补偿方法

胡程 邓云开 田卫明 曾涛

引用本文:
Citation:

地基干涉合成孔径雷达图像非线性大气相位补偿方法

    作者简介: 胡 程(1981–) 男,湖南岳阳人,博士,北京理工大学,研究员。主要研究方向为新体制合成孔径雷达系统与信号处理、生物探测雷达系统与信息处理技术等;在IEEE Trans. on GRS, IEEE Trans. on AES, J-GR, Remote Sensing等国内外期刊发表SCI论文60余篇,国内外会议发表EI论文100余篇;授权发明专利30余项。E-mail: cchchb@163.com;邓云开(1992–) 男,河南周口人,博士,北京理工大学。主要研究为地基SAR高精度1维形变与3维形变实时测量算法,目前已发表SCI论文5篇,其他论文7篇。E-mail: yunkai_bit@foxmail.com.
    通讯作者: 邓云开, yunkai_bit@foxmail.com
  • 基金项目:

    国家自然科学基金(61427802, 61601031)

  • 中图分类号: TN95

A Compensation Method of Nonlinear Atmospheric Phase Applied for GB-InSAR Images

    Corresponding author: DENG Yunkai, yunkai_bit@foxmail.com ;
  • Fund Project: The National Natural Science Foundation of China (61427802, 61601031)

    CLC number: TN95

  • 摘要: 采用基于永久散射体(PS)技术的大气相位(AP)补偿方法对地基干涉合成孔径雷达(GB-InSAR) 干涉相位图进行分析时,部分图像的干涉相位随距离呈现出复杂的非线性,导致无法建立合理的多参数模型来模拟大气相位,常规的补偿方法不再适用。该文提出一种改进的GB-InSAR图像非线性大气相位补偿方法,首先采用常规方法对干涉相位图进行大气相位补偿,并根据PS点的补偿后相位序列的标准差,进行稳定PS点的选择,然后对稳定PS点进行子区域划分,通过反距离加权插值估计出所有PS点的大气相位,从而实现大气相位的有效补偿。采用该文所提方法,对460幅雷达图像进行了处理,相比于常规方法,可以有效地补偿干涉相位图中的非线性大气相位分量。基于若干高稳定参考点的对比结果表明,该方法减小了最大约1 rad的相位测量误差。
  • 图 1  场景信息

    Figure 1.  Scene information

    图 2  MIMO雷达系统照片

    Figure 2.  Photo of the MIMO radar system

    图 3  MIMO雷达图像与PS图

    Figure 3.  MIMO radar image and PS map

    图 4  长时间基线干涉相位图

    Figure 4.  Interferometric phase map with long temporal baseline

    图 5  短时间基线干涉相位图A和B

    Figure 5.  Interferometric phase maps of A and B with short temporal baselines

    图 6  形变PS点选择

    Figure 6.  Selection of deformation PS

    图 7  非线性大气相位补偿

    Figure 7.  Non-linear atmospheric phase compensation

    图 8  参考点相位曲线对比

    Figure 8.  Phase comparison curves of reference points

    图 9  累积相位图

    Figure 9.  Cumulative phase maps

  • [1] 刘斌, 葛大庆, 李曼, 等. 地基合成孔径雷达干涉测量技术及其应用[J]. 国土资源遥感, 2017, 29(1): 1–6. doi: 10.6046/gtzyyg.2017.01.01LIU Bin, GE Daqing, LI Man, et al. Ground-based interferometric synthetic aperture radar and its applications[J]. Remote Sensing for Land &Resources, 2017, 29(1): 1–6. doi: 10.6046/gtzyyg.2017.01.01
    [2] 曾涛, 邓云开, 胡程, 等. 地基差分干涉雷达发展现状及应用实例[J]. 雷达学报, 2019, 8(1): 154–170. doi: 10.12000/JR18115ZENG Tao, DENG Yunkai, HU Cheng, et al. Development state and application examples of ground-based differential interferometric radar[J]. Journal of Radars, 2019, 8(1): 154–170. doi: 10.12000/JR18115
    [3] 董杰, 董妍. 基于气象数据的地基雷达大气扰动校正方法研究[J]. 测绘工程, 2014, 23(10): 72–75. doi: 10.3969/j.issn.1006-7949.2014.10.017DONG Jie and DONG Yan. Atmospheric artifact compensation for deformation monitoring with ground-based radar[J]. Engineering of Surveying and Mapping, 2014, 23(10): 72–75. doi: 10.3969/j.issn.1006-7949.2014.10.017
    [4] IANNINI L and GUARNIERI A M. Atmospheric phase screen in ground-based radar: Statistics and compensation[J]. IEEE Geoscience and Remote Sensing Letters, 2011, 8(3): 537–541. doi: 10.1109/LGRS.2019.1590876
    [5] 黄长军, 夏红梅, 周吕. 基于GCP方法的地基InSAR大气扰动误差改正分析[J]. 测绘与空间地理信息, 2018, 41(10): 8–11. doi: 10.3969/j.issn.1672-5867.2018.10.003HUANG Changjun, XIA Hongmei, and ZHOU Lyu. Atmospheric disturbance error correction in GB-InSAR based on ground control point[J]. Geomatics &Spatial Information Technology, 2018, 41(10): 8–11. doi: 10.3969/j.issn.1672-5867.2018.10.003
    [6] 徐亚明, 周校, 王鹏, 等. GB-SAR构建永久散射体网改正气象扰动方法[J]. 武汉大学学报: 信息科学版, 2016, 41(8): 1007–1012, 1020. doi: 10.13203/j.whugis20140507XU Yaming, ZHOU Xiao, WANG Peng, et al. A method of constructing permanent scatterers network to correct the meteorological disturbance by GB-SAR[J]. Geomatics and Information Science of Wuhan University, 2016, 41(8): 1007–1012, 1020. doi: 10.13203/j.whugis20140507
    [7] NOFERINI L, PIERACCINI M, MECATTI D, et al. Permanent scatterers analysis for atmospheric correction in ground-based SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(7): 1459–1471. doi: 10.1109/tgrs.2005.848707
    [8] IGLESIAS R, FABREGAS X, AGUASCA A, et al. Atmospheric phase screen compensation in ground-based SAR with a multiple-regression model over mountainous regions[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(5): 2436–2449. doi: 10.1109/TGRS.2019.1590876
    [9] HU Cheng, DENG Yunkai, TIAN Weiming, et al. A PS processing framework for long-term and real-time GB-SAR monitoring[J]. International Journal of Remote Sensing, 2019, 40(16): 6298–6314. doi: 10.1080/01431161.2019.1590876
    [10] 张祥, 陆必应, 宋千. 地基SAR差分干涉测量大气扰动误差校正[J]. 雷达科学与技术, 2011, 9(6): 502–506, 512. doi: 10.3969/j.issn.1672-2337.2011.06.004ZHANG Xiang, LU Biying, and SONG Qian. Atmospheric disturbance correction in ground-based SAR differential interferometry[J]. Radar Science and Technology, 2011, 9(6): 502–506, 512. doi: 10.3969/j.issn.1672-2337.2011.06.004
    [11] HU Cheng, ZHU Mao, ZENG Tao, et al. High-precision deformation monitoring algorithm for GBSAR system: Rail determination phase error compensation[J]. Science China Information Sciences, 2015, 59(8): 082307. doi: 10.1007/s11432-015-5446-z
    [12] 张建萍, 刘希玉. 基于聚类分析的K-means算法研究及应用[J]. 计算机应用研究, 2007, 24(5): 166–168. doi: 10.3969/j.issn.1001-3695.2007.05.051ZHNAG Jianping and LIU Xiyu. Application in cluster′s analysis is analyzed in children development period[J]. Application Research of Computers, 2007, 24(5): 166–168. doi: 10.3969/j.issn.1001-3695.2007.05.051
    [13] 朱吉祥, 张礼中, 周小元, 等. 反距离加权法在区域滑坡危险性评价中的应用[J]. 水土保持通报, 2012, 32(3): 136–140. doi: 10.13961/j.cnki.stbctb.2012.03.009ZHU Jixiang, ZHANG Lizhong, ZHOU Xiaoyuan, et al. Application of inverse distance weighted method to regional landslide hazards assessment[J]. Bulletin of Soil and Water Conservation, 2012, 32(3): 136–140. doi: 10.13961/j.cnki.stbctb.2012.03.009
    [14] 刘作利, 刘景玉, 申修强, 等. 唐山马兰庄铁矿露天开采边坡变形监测的GB-InSAR技术[J]. 现代矿业, 2018(4): 165–170. doi: 10.3969/j.issn.1674-6082.2018.04.047LIU Zuoli, LIU Jingyu, SHEN Xiuqiang, et al. Deformation monitoring of the open-pit slope of Malanshan iron mine in Tangshan city based on GB-InSAR[J]. Modern Mining, 2018(4): 165–170. doi: 10.3969/j.issn.1674-6082.2018.04.047
    [15] TIAN Weiming, ZHAO Zheng, HU Cheng, et al. GB-InSAR-based DEM generation method and precision analysis[J]. Remote Sensing, 2019, 11(9): 997. doi: 10.3390/rs11090997
    [16] HU Cheng, WANG Jingyang, TIAN Weiming, et al. Design and imaging of ground-based multiple-input multiple-output synthetic aperture radar (MIMO SAR) with non-collinear arrays[J]. Sensors, 2017, 17(3): 598. doi: 10.3390/s17030598
  • [1] 王佳宁许小剑 . 二维线性与非线性海面的宽带散射特性仿真及分析. 雷达学报, doi: 10.12000/JR15053
    [2] 赵雨露张群英李超纪奕才方广有 . 视频合成孔径雷达振动误差分析及补偿方案研究. 雷达学报, doi: 10.12000/JR14153
    [3] 任笑真杨汝良 . 一种基于幅度和相位迭代重建的四维合成孔径雷达成像方法. 雷达学报, doi: 10.12000/JR15135
    [4] 陈曾平张炜承林钱强 . 宽带雷达ISAR 成像相位补偿新方法(英文). 雷达学报, doi: 10.3724/SP.J.1300.2013.13023
    [5] 钟雪莲向茂生郭华东陈仁元 . 机载重轨干涉合成孔径雷达的发展. 雷达学报, doi: 10.3724/SP.J.1300.2013.13005
    [6] 杜剑波李道京马萌 . 机载合成孔径激光雷达相位调制信号性能分析和成像处理. 雷达学报, doi: 10.3724/SP.J.1300.2014.13094
    [7] 李道京胡 烜 . 合成孔径激光雷达光学系统和作用距离分析. 雷达学报, doi: 10.12000/JR18017
    [8] 杜强宋耀良曹晓健 . 基于Hermite 插值滤波器的直接延时补偿超宽带波束形成技术研究. 雷达学报, doi: 10.3724/SP.J.1300.2013.13028
    [9] 吴一戎 . 多维度合成孔径雷达成像概念. 雷达学报, doi: 10.3724/SP.J.1300.2013.13047
    [10] 高敬坤邓彬秦玉亮王宏强黎湘 . THz全尺寸凸体粗糙目标雷达回波散射建模与成像仿真. 雷达学报, doi: 10.12000/JR17086
    [11] 于彬彬刘畅王岩飞 . 基于chirp-z 变换的斜视FMCW-SAR 非线性尺度变换算法. 雷达学报, doi: 10.3724/SP.J.1300.2012.20005
    [12] 詹学丽王岩飞王超李和平 . 一种用于合成孔径雷达的数字去斜方法. 雷达学报, doi: 10.12000/JR14117
    [13] 林世斌李悦丽严少石周智敏 . 平地假设对合成孔径雷达时域算法成像质量的影响研究. 雷达学报, doi: 10.3724/SP.J.1300.2012.20035
    [14] 王岩飞刘畅詹学丽韩松 . 无人机载合成孔径雷达系统技术与应用. 雷达学报, doi: 10.12000/JR16089
    [15] 李海英张珊珊李世强张华春 . 环境一号C 卫星合成孔径雷达相干性分析. 雷达学报, doi: 10.3724/SP.J.1300.2014.13060
    [16] 曾操梁思嘉王威徐青 . 基于频率步进信号的旋转式合成孔径雷达成像方法. 雷达学报, doi: 10.3724/SP.J.1300.2014.14043
    [17] 路满宋红军罗运华 . 基于调频连续波信号的圆弧式合成孔径雷达成像方法. 雷达学报, doi: 10.12000/JR16007
    [18] 曾涛 . 双基地合成孔径雷达发展现状与趋势分析. 雷达学报, doi: 10.3724/SP.J.1300.2012.20093
    [19] 杨汝良戴博伟李海英 . 极化合成孔径雷达极化层次和系统工作方式. 雷达学报, doi: 10.12000/JR16013
    [20] 金添 . 叶簇穿透合成孔径雷达增强成像方法. 雷达学报, doi: 10.12000/JR15114
  • 加载中
图(9)
计量
  • 文章访问数:  178
  • HTML浏览量:  94
  • PDF下载量:  20
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-08-12
  • 录用日期:  2019-11-08
  • 网络出版日期:  2019-11-25

地基干涉合成孔径雷达图像非线性大气相位补偿方法

    通讯作者: 邓云开, yunkai_bit@foxmail.com
    作者简介: 胡 程(1981–) 男,湖南岳阳人,博士,北京理工大学,研究员。主要研究方向为新体制合成孔径雷达系统与信号处理、生物探测雷达系统与信息处理技术等;在IEEE Trans. on GRS, IEEE Trans. on AES, J-GR, Remote Sensing等国内外期刊发表SCI论文60余篇,国内外会议发表EI论文100余篇;授权发明专利30余项。E-mail: cchchb@163.com;邓云开(1992–) 男,河南周口人,博士,北京理工大学。主要研究为地基SAR高精度1维形变与3维形变实时测量算法,目前已发表SCI论文5篇,其他论文7篇。E-mail: yunkai_bit@foxmail.com
  • ①. 北京理工大学信息与电子学院雷达技术研究所 北京 100081
  • ②. 北京理工大学卫星导航电子信息技术教育部重点实验室 北京 100081
基金项目:  国家自然科学基金(61427802, 61601031)

摘要: 采用基于永久散射体(PS)技术的大气相位(AP)补偿方法对地基干涉合成孔径雷达(GB-InSAR) 干涉相位图进行分析时,部分图像的干涉相位随距离呈现出复杂的非线性,导致无法建立合理的多参数模型来模拟大气相位,常规的补偿方法不再适用。该文提出一种改进的GB-InSAR图像非线性大气相位补偿方法,首先采用常规方法对干涉相位图进行大气相位补偿,并根据PS点的补偿后相位序列的标准差,进行稳定PS点的选择,然后对稳定PS点进行子区域划分,通过反距离加权插值估计出所有PS点的大气相位,从而实现大气相位的有效补偿。采用该文所提方法,对460幅雷达图像进行了处理,相比于常规方法,可以有效地补偿干涉相位图中的非线性大气相位分量。基于若干高稳定参考点的对比结果表明,该方法减小了最大约1 rad的相位测量误差。

English Abstract

    • 作为一种高精度的形变测量仪器,地基干涉合成孔径雷达(Ground-Based Interferometric Synthetic Aperture Radar, GB-InSAR)已经在形变监测领域得到了广泛的应用[1]。GB-InSAR通常是基于差分干涉测量技术,通过对同一位置、不同时刻获取的两幅SAR图像进行差分干涉处理,基于相位信息来实现形变测量,其一般工作在X或者Ku波段,形变测量精度可以达到毫米或者亚毫米量级。GB-InSAR测量误差的主要来源是大气相位,由于不同时刻气象条件的不同,电磁波在大气中传播的速度会发生改变,从而导致大气延时误差[2]

      为实现大气相位的补偿,国内外学者主要提出了3种解决方法。第1种的补偿方法是在观测场景内建立气象站,基于大气折射率模型,利用气象数据(温度、湿度、大气压)来对大气相位进行定量估计[3];第2种是在场景中人共布设或者选择出一些高度稳定的特征点,采用空间插值的方法实现对整幅图像的大气相位的补偿[4,5];第3种则是基于永久散射体 (Permanent Scatterer, PS)技术,根据大气相位的空间分布特征,建立描述大气相位的方程,估计大气相位参数,实现大气相位的补偿[6]

      在基于PS技术进行大气相位补偿时,首先需要选择出未处在形变区域的PS点,然后对大气相位进行合理的建模,并建立线性方程组,实现对大气参数的粗估计,此后剔除与模型偏差较大的PS点,逐步迭代实现大气参数的准确估计。该方法不需要气象参数及布设特征点,可以基于大量的PS点进行大气参数的估计,估计精度较高。一般情况下,大气在空间上均匀变化,可以将大气相位建模为随斜距线性变化的分量[7]。在地形陡峭的山区,大气在空间上非均匀变化,可以采用兼顾斜距和高程的多参数模型[8]

      基于PS技术的大气相位补偿方法,已经在GB-InSAR领域得到了广泛的应用,但还存在一些典型的问题。首先,该方法要求采用场景中的非形变PS点建立观测方程组;其次,气象条件一直在随时间改变,大气相位的时变性很强,在较差天气条件(降雨、降雪、强风等)下,大气在空间上非均匀变化,导致大气相位可能表现出复杂的空变性,无法建立合理的多参数模型来模拟大气相位。因此,常规的基于PS技术的补偿方法,在应用于时序GB-InSAR图像处理时,还存在较大的改进空间。

      针对上述问题,本文借鉴在场景中布设若干个稳定的地面控制点,通过空间维插值估计大气相位的思想。首先采用常规补偿方案对所有的干涉相位图进行大气相位补偿,并分析PS点的相位序列的标准差,设定门限选择出稳定PS点;然后基于K均值聚类算法 (K-means clustering algorithm, K-means)将稳定PS点划分出一定数量的子区域,将每一个子区域的中心点设定为控制点,采用反距离加权插值算法估计所有PS点的大气相位,从而实现非线性大气相位的补偿。

      本文详细介绍了补偿方法的实现流程,主要分为稳定PS点选择和空间维插值补偿两部分,然后分别采用改进方法和常规方法对460幅地基多输入多输出 (Multiple-Input Multiple-Output, MIMO)雷达图像进行了分析,对比验证了本文方法的有效性。

    • 在利用像素点的相位信息进行形变测量时,差分干涉相位的质量直接影响到形变测量的精度。因此,GB-InSAR差分干涉处理时,通常需要选择出一些高质量的像素点,即PS点,来进行形变分析。在GB-InSAR领域,广泛采用幅度离差法来进行PS点的选择,一个像素点的时序幅度序列的标准差与均值之比,即为幅度离差指数。通过设置合理的幅度离差门限,即可以实现PS点的选择[9]

      一个PS点的差分干涉相位$\Delta \varphi $可以建模为

      $ \Delta \varphi = {\varphi _{{\rm{defo}}}} + {\varphi _{{\rm{atm}}}} + {\varphi _{{\rm{noi}}}} + 2k\pi $

      其中,${\varphi _{{\rm{defo}}}}$为形变相位分量;${\varphi _{{\rm{atm}}}}$为两幅图像获取期间,由大气条件改变所导致的大气相位分量;${\varphi _{{\rm{noi}}}}$为噪声相位分量;$2k\pi $为相位模糊度,$k$为整数。由于相位周期性的影响,$\Delta \varphi $处在$ - \pi $$\pi $的范围内[10]。在进行大气相位补偿前,要对干涉相位图进行相位解缠,可以采用非均匀网格下的最小费用流算法,下文中采用$\Delta \varphi $指代解缠相位。

      大气相位${\varphi _{{\rm{atm}}}}$可以建模为

      $ {\varphi _{{\rm{atm}}}} = \frac{{4\pi }}{\lambda }\int_L {\Delta N\left( {{{r}},t} \right){\rm d}L} $

      其中,$\lambda $为信号的波长,$\Delta N$表示折射率的变化,其随着空间${{r}}$和时间$t$发生变化,$L$表示信号的传输路径。一般情况下,大气的空间同质性很好,可以假设$\Delta N$在空间${{r}}$上不发生变化,在时间$t$上随机变化,因此可以将${\varphi _{{\rm{atm}}}}$建模为随斜距线性变化的分量

      $ {\varphi _{{\rm{atm}}}} = \frac{{4\pi }}{\lambda }\Delta N \cdot R = \frac{{4\pi }}{\lambda }\left( {{\beta _0} + {\beta _1} \cdot R} \right) $

      其中,${\beta _{0}}$为常数分量,${\beta _1}$表示与斜距相关的线性系数,$R$表示雷达与目标点之间的距离。基于PS技术进行大气相位补偿时,首先建立线性方程组

      $\Delta {{\varPhi}} = {{X \beta}} + {{\varepsilon}}$

      其中,$\Delta {{\varPhi}} = \left[ {\begin{array}{*{20}{c}} {\Delta {\varphi _1}} \\ {\Delta {\varphi _2}} \\ \vdots \\ {\Delta {\varphi _n}} \end{array}} \right]$, $ {{X}} = \dfrac{{4\pi }}{\lambda }\left[ {\begin{array}{*{20}{c}} 1 & {{R_1}} \\ 1 & {{R_2}} \\ \vdots & \vdots \\ 1 & {{R_n}} \end{array}} \right]$, ${{\beta}} = \left[ {\begin{array}{*{20}{c}} {{\beta _0}} \\ {{\beta _1}} \end{array}} \right]$, ${{\varepsilon}} = \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}} \\ {{\varepsilon _2}} \\ \vdots \\ {{\varepsilon _n}} \end{array}} \right]$

      $\Delta {{\varPhi}} $$n$个PS点的解缠相位$\Delta {\varphi _1},\;\Delta {\varphi _2},\; ··· , $$\Delta {\varphi _n}$构成的$n \times 1$维向量,${{X}}$为常数1与$n$个PS点的斜距${R_1},\;{R_2}, ··· ,{R_n}$构成的$n \times 2$维矩阵,${{\beta}} $为待估计的$2 \times 1$维向量,${{\varepsilon}} $$n \times 1$维的随机误差向量,${\varepsilon _1},\;{\varepsilon _2},\; ··· ,{\varepsilon _n}$为各PS点的模型误差相位分量。采用最小二乘法对${{\beta}} $进行估计,可以得到

      $\hat {{\beta}} = {\left( {{{{X}}^{\rm{T}}}{{X}}} \right)^{ - 1}} \cdot {{{X}}^{\rm{T}}}\Delta {{\varPhi}} $

      其中,${\rm{T}}$表示矩阵的转置。大气相位的估计分量为

      $\Delta {{{\varPhi}} _{{\rm{APS}}}} = {{X}}\hat {{\beta}} $

      $\Delta {{\varPhi}} $$\Delta {{{\varPhi}} _{{\rm{APS}}}}$之间的差值即为补偿后相位。首先基于所有的PS点对${{\beta}} $进行估计,然后为提高估计的精度,剔除一些不可靠PS点。不满足式(7)准则的PS点即为不可靠PS点,$\Delta {T_{{\rm{atm}}}}$的取值在0.1~0.2 rad范围内。剔除不可靠PS点后,基于剩余的PS点对大气相位参数进行二次估计[11]

      $\left| {\Delta {{\varPhi}} - \Delta {{\hat {{\varPhi}} }_{{\rm{APS}}}}} \right| < \Delta {{{T}}_{{\rm{atm}}}}$

      除了最基本的线性斜距模型外,常用的参数模型还包括高阶斜距模型、斜距-方位角模型、斜距-高程模型等,如式(8)—式(10)。其中${\beta _0}$, ${\beta _1}$${\beta _2}$是各模型中待估计的未知参数,$R$, $\theta $$h$分别代表目标点的斜距、方位角和高程。对线性方程组式(4)修正后,即可以进行大气相位补偿

      $ \varphi _{{\rm{atm}}}^{(1)} = \frac{{4\pi }}{\lambda }\left( {{\beta _0} + {\beta _1} \cdot R + {\beta _2} \cdot {R^2}} \right)\ $

      $ \varphi _{{\rm{atm}}}^{(2)} = \frac{{4\pi }}{\lambda }\left( {{\beta _0} + {\beta _1} \cdot R + {\beta _2} \cdot \sin \theta } \right) $

      $ \varphi _{{\rm{atm}}}^{(3)} = \frac{{4\pi }}{\lambda }\left( {{\beta _0} + {\beta _1} \cdot R + {\beta _2} \cdot R \cdot h} \right) $

      基于PS技术的补偿方法,其有效补偿的条件,首先是要求所采用的参数模型可以准确模拟大气相位,其次是迭代估计参数时,可以有效剔除不可靠PS点。采用上述方法对一幅干涉相位图进行大气相位补偿后,一个PS点的补偿后相位$\Delta {\varphi _{{\rm{APC}}}}$可以表示为

      $\Delta {\varphi _{{\rm{APC}}}} = {\varphi _{{\rm{defo}}}} + \varphi _{{\rm{atm}}}^{\left( {\rm{e}} \right)} + {\varphi _{{\rm{noi}}}}$

      其中,$\varphi _{{\rm{atm}}}^{\left( {\rm{e}} \right)}$表示大气相位残余误差分量。考虑到式(11)中的3个相位分量彼此不相关,一个PS点的补偿后相位序列的标准差${\sigma _{{\rm{APC}}}}$可以表示为

      ${\sigma _{{\rm{APC}}}} = \sqrt {{{\left( {{\sigma _{{\rm{defo}}}}} \right)}^2}{\rm{ + }}{{\left( {\sigma _{{\rm{atm}}}^{\left( {\rm{e}} \right)}} \right)}^2}{\rm{ + }}{{\left( {{\sigma _{{\rm{noi}}}}} \right)}^2}} $

      其中,${\sigma _{{\rm{defo}}}}$, $\sigma _{{\rm{atm}}}^{\left( {\rm{e}} \right)}$${\sigma _{{\rm{noi}}}}$分别对应形变相位分量、大气相位残余分量和噪声相位分量的标准差。在图像的数目足够多时,通过合理的设定相位门限,将标准差大于该门限的PS点剔除,即相位序列中包含形变相位${\varphi _{{\rm{defo}}}}$或者噪声相位分量${\varphi _{{\rm{noi}}}}$较大的PS点,将剩余的PS点定义为稳定PS点。

    • 在选择出稳定PS点后,对于每一幅干涉相位图,分别进行大气相位补偿。由于稳定PS点的相位分量中包含噪声相位和大气相位,考虑到噪声相位在干涉相位图上随机变化,不具备空间相关性,而大气相位虽然会在整幅图像范围内发生变化,但在较小的距离范围内可以视为一个常数。因此如果对空间上较小距离范围内的所有PS点进行相位平均,则噪声相位可以得到很好的滤除。由于稳定PS点是非均匀的分布在整幅图像范围内,可以基于K-means算法对稳定PS点进行簇划分,将每一个簇定义为一个子区域。每一个簇的中心点作为控制点,并对该簇内的所有PS点进行相位平均,作为当前控制点的相位。

      K-means算法的实现原理是对于给定的样本集x,按照样本之间的距离大小,将样本集划分为K个簇,让各个簇内的点的距离尽可能的小,而让各个簇之间的点的距离尽可能的大[12]。假设将样本集x划分为K个簇$\left( {{C_1},{C_2}, ··· ,{C_K}} \right)$,各簇间的平方误差和E可以表示为

      $E = \sum\limits_{i = 1}^K {\sum\limits_{x \in {C_i}} {{{\left\| {x - {{{\mu}} _i}} \right\|}^2}} } $

      其中,$\left\| {\cdot} \right\|$表示2阶范数,即向量的模。${{{\mu}} _i}$是簇${C_i}$的均值向量,可以表示为

      ${{{\mu}} _i} = \frac{1}{{\left| {{C_i}} \right|}}\sum\limits_{x \in {C_i}} x $

      其中,$\left| {\cdot} \right|$表示1阶范数,即簇中点的数量。

      经过K-means划分后可以得到K个子区域,将每一个子区域中心点定义为控制点${\rm{C}}{{\rm{P}}_{{\rm{KM}}}}$,其相位${\bar \varphi _{{\rm{KM}}}}$为该子区域内所有PS点的相位均值。

    • 基于这K个控制点${\rm{C}}{{\rm{P}}_{{\rm{KM}}}}$的平均相位${\bar \varphi _{{\rm{KM}}}}$,进行空间维插值来估计所有PS点的大气相位,可以采用反距离加权插值算法。反距离加权插值算法是利用已知点与待插点之间的距离来定义加权因子,然后加权计算待插点的相位,距离越近加权比重越大。其计算公式为

      $ {Z_{\left( {x,y} \right)}} = {{\sum\limits_{i = 1}^n {\frac{{{Z_i}}}{{{{\left| {{d_i}} \right|}^\mu }}}} } / {\sum\limits_{i = 1}^n {\frac{1}{{{{\left| {{d_i}} \right|}^\mu }}}} }} $

      其中,${Z_{\left( {x,y} \right)}}$为插值结果,$\left( {x,y} \right)$为待插点的空间坐标,${Z_i}$为第$i\left( {i = 1,2, ··· ,n} \right)$个参考点的数值,$\left| {{d_i}} \right|$表示待插点与第$i$个参考点之间的空间距离,$\mu $表示加权因子的幂指数,一般取为2[13]

      将这些控制点视为一组离散点,可以基于Delaunay三角剖分法则来构建一个三角形网络。对于该网络内的任意一个三角形,其由3个控制点作为顶点,且外接圆中不包含其他控制点。在基于反距离加权插值算法来估计所有PS点的大气相位时,如果一个PS点处在某一个三角形内部,则以该三角形的3个控制点作为参考点,基于式(15)估计当前PS点的大气相位。如果一个PS点处在所有三角形的外部,则选择最近的3个控制点,同样基于式(15)来估计大气相位。由于每一个控制点的相位,均是通过对一个子区域内的所有PS点进行相位平均获得,则上述方法可以视为采用了大量的PS点作为参考点,空间维插值精度可以得到保证。

    • 本实验选定的区域为一露天开采矿坑(E$118^\circ 36'23''$, N$40^\circ 06'44''$),其位于河北省迁安市马兰庄镇。图1(a)所示为场景照片,其中红色椭圆所示为场景的形变区域。图1(b)所示为场景卫星图,矿坑正上方为椭圆形,其长轴约1050 m,短轴约680 m,其中黄色矩形代表雷达的布放位置,雷达的观测视角范围为$60^\circ $。矿坑边坡为典型岩质边坡,基本无植被覆盖,最大开采深度约400 m,边坡倾角为$38^\circ \sim 47^\circ $[14,15]

      图  1  场景信息

      Figure 1.  Scene information

      实验中采用一部MIMO雷达对该矿坑进行了监测,连续获取了460幅雷达图像,时间从2018年10月26日17:30~2018年10月27日11:30。该MIMO雷达采用16个发射天线构成两个密集子阵列,16个接收天线构成一个稀疏子阵列,如图2所示。其工作在Ku波段,波长为$\lambda = 1.86\;{\rm{cm}}$,等效合成孔径长度为1.138 m,角分辨率为$0.466^\circ $[16]

      图  2  MIMO雷达系统照片

      Figure 2.  Photo of the MIMO radar system

      图3(a)所示为该露天矿坑在极坐标系下的成像结果,边坡区域内像素点的幅值主要分布在–30~0 dB的范围内。基于幅度离差法进行PS点的选择,设置幅度离差门限0.15,幅度门限–25 dB,可筛选出61764个PS点,如图3(b)所示。可以看出,坡体上大部分像素点的幅度稳定性很高,仅中间的道路上由于雷达观测视角的原因,未能有效选择出PS点。

      图  3  MIMO雷达图像与PS图

      Figure 3.  MIMO radar image and PS map

      以第1幅图像为主图像,最后一幅图像为辅图像,获取差分干涉相位图,如图4(a)所示,将该相位图反投到3维地形图上,如图4(b)所示。借助于3维地形,可以更加准确地确定形变区域的发生位置。图4(c)所示为PS点的干涉相位图,显然大部分PS点的相位在0 rad左右,而红色椭圆标识出的一块区域,其干涉相位与其他区域的明显不同,该区域为发生了较大形变的形变区域。图4(d)所示为PS点的干涉相位随其斜距变化的分布图,如果将大气相位建模为随斜距线性变化的相位分量,其中红色实线所示为线性大气相位估计结果,显然存在很大的误差。因此,对较长时间基线的干涉相位图进行大气相位补偿时,如果不能合理的剔除形变PS点,将会严重影响到对大气参数的估计。在处理时序GB-InSAR图像时,为提高大气相位补偿的准确度,进行差分干涉的两幅雷达图像之间的时间基线不宜过大。

      图  4  长时间基线干涉相位图

      Figure 4.  Interferometric phase map with long temporal baseline

    • 分析这460幅时序MIMO雷达图像时,采用相邻的两幅图像构成一个干涉图像对,则可以获取459幅干涉相位图。每一幅干涉相位图的时间基线仅为2~3 min,通常在这么短的时间内,大气条件的变化很小,相应的大气相位误差分量也很小。图5(a)图5(b)所示为干涉相位图A和其相位散点图,所有PS点的干涉相位均在0 rad左右。在采用线性斜距模型补偿大气相位时,最大的补偿分量不超过0.05 rad。但在部分时间段,大气在空间上不再是均匀变化,导致大气相位呈现出非线性变化。图5(c)图5(d)所示为干涉相位图B和其相位散点图,显然PS点的干涉相位变化比较复杂,采用线性斜距模型补偿大气相位时,最大会带来约0.5 rad的补偿误差,相应的形变测量误差约为0.74 mm,会严重影响到MIMO雷达的形变测量精度。即使是采用其他模型对相位图B进行分析,依然无法有效的补偿大气相位。

      图  5  短时间基线干涉相位图A和B

      Figure 5.  Interferometric phase maps of A and B with short temporal baselines

      对459幅干涉相位图进行相位解缠,并采用线性斜距模型进行大气相位补偿,然后计算每一个PS点的补偿后相位序列的标准差。图6(a)所示为所有PS点的相位标准差图,显然形变PS点的标准差很大,设置0.3 rad的标准差门限,筛选出4341个PS点。这些PS点的分布情况图如图6(b)所示,除红色椭圆形变区域内的形变PS点外,还有少量的形变区域外的噪声PS点也被选择,如紫色椭圆所标识出的部分点。

      图  6  形变PS点选择

      Figure 6.  Selection of deformation PS

      经过相位标准差筛选后,得到57423个稳定PS点,设定每个子区域中的PS点平均数量为200个,则可以划分出287个子区域。以图5(c)所示的干涉相位图B为例来说明本文所提方法,图7(a)中红色方形点标识出了每一个子区域中心,即287个控制点,以干涉相位图B作为背景,图7(b)为局部放大图,其中每一个黑色多边形,代表各个子区域中PS点的最小外接多边形。之后对于各个子区域中的PS点进行相位平均,获取控制点的平均相位,采用反距离加权法估计所有PS点的大气相位,如图7(c)所示,可以很直观地看出,估计图中PS点的相位随其斜距呈现出非线性变化,和图5(c)高度相似。图7(d)所示则为补偿结果,所有PS点的补偿后相位均在0 rad左右,空变性的大气相位得到了有效的补偿。

      图  7  非线性大气相位补偿

      Figure 7.  Non-linear atmospheric phase compensation

      采用本文所提非线性补偿方法和常规的线性斜距模型补偿方法,对这459幅相位图分别分析,然后从不同斜距处选出5个幅度离差指数最小的PS点作为参考点,来对比说明本文方法的有效性。图8(a)所示为采用本文方法获取的累积相位图,由459幅补偿后的干涉相位图累加获取,该幅图像可以用来反映整个监测周期内,场景中PS点的相位变化情况。形变区域PS点的相位变化达到了约–3 rad,非形变区域PS点的相位变化在±1 rad范围内,直观上很难看出非形变区域PS点的相位是否存在空变性。图8(a)中红色方形点标识出的点1~点5为非形变区域的参考点。图8(b)所示为这5个参考点的时序相位变化曲线,点1~点5的相位变化趋势高度相似,这是由于大气相位具有较高的空间相关性。图8(c)图8(b)所示分别为采用本文方法和常规方法对点1~点5的补偿结果。常规方法的补偿结果中,点1的相位随着时间逐渐增大,这说明了大气相位残余误差随着时间逐渐积累了起来。本文方法的补偿结果中,各点的相位序列变化更为随机。

      图  8  参考点相位曲线对比

      Figure 8.  Phase comparison curves of reference points

      图9(a)图9(b)所示为采用常规方法和本文方法获取的累积相位图,为了更直观地对比场景中PS点的相位变化情况,两幅图像中均忽略了形变区域。可以很直观地看出,常规方法的补偿相位图中,PS点的相位依然呈现出明显的空变性,尤其是图像中500~600 m范围内的区域,该区域PS点的相位达到了1 rad,相应的形变量为1.48 mm。由于大气相位补偿误差随着时间积累起来的原因,这部分区域很容易被误认为在雷达监测周期内,也出现了形变。

      图  9  累积相位图

      Figure 9.  Cumulative phase maps

    • 本文提出了一种适用于时序GB-InSAR图像非线性大气相位补偿的改进方法,解决了常规方法中多参数模型无法准确模拟大气相位,从而无法对部分干涉相位图进行有效补偿的问题。采用本文方法与常规方法,对一露天矿坑的460幅时序地基MIMO雷达图像分别进行了处理,对比了参考点的时序相位曲线和累积相位图,验证了本文方法的有效性。

      本文方法还存在着一些不足。首先本文方法要求先采用常规方法对时序干涉相位图进行补偿,基于PS点的补偿后相位序列的标准差大小来选择稳定PS点,但相位标准差门限是人为设定的;其次如果部分PS点仅在少量图像中发生了形变,这些PS点的相位标准差可能较小,无法有效地筛选出;最后本文方法目前是应用于对时序GB-InSAR图像的事后大气相位补偿,实际形变监测中更关注的是实时处理,需要进一步研究实时形变处理中的非线性大气相位补偿方法。

参考文献 (16)

目录

    /

    返回文章
    返回