一种基于最大后验框架的聚类分析多基线干涉SAR高度重建算法

斯奇 王宇 邓云凯 李宁 张衡

引用本文:
Citation:

一种基于最大后验框架的聚类分析多基线干涉SAR高度重建算法

    作者简介: 王   宇(1980–),男,河南人,现为中国科学院电子学研究所研究员,博士生导师,研究方向为SAR系统设计与信号处理技术。E-mail: yuwang@mail.ie.ac.cn;.
    通讯作者: 斯奇, si_qi0616@163.com
  • 基金项目:

    国家自然科学基金优秀青年基金(61422113),国家万人计划青年拔尖人才,中科院百人计划

A Novel Cluster-Analysis Algorithm Based on MAP Framework for Multi-baseline InSAR Height Reconstruction

    Corresponding author: Si Qi, si_qi0616@163.com ;
  • Fund Project: The National Natural Science Foundation of China (61422113), The National Ten Thousand Talent ProgramYoung Top Notch Talent Program, The Hundred Talents Program of the Chinese Academy of Sciences

  • 摘要: 多基线干涉SAR能有效减小由目标高度急剧变化和较大噪声干扰带来的不利影响,可以获取比单基线干涉SAR更精确的地表数字高程模型(DEM)。传统的基于最大似然估计(ML)的多基线高度重建算法在通道数目较少情况下重建结果不佳,基于最大后验估计(MAP)的多基线高度重建算法存在运行时间较长的问题,针对以上问题,该文提出了一种基于最大后验框架的聚类分析高度重建算法(CABMAP)。该算法首先利用了ML估计法得到粗略的DEM,以此为基础在每次迭代过程中利用聚类分析(CA)判断出邻域内的噪声像素,并通过计算后验概率完成重建,此外采用了一种改进措施提高精度。这样,既保留了ML估计法运行速度快的特征,又具有MAP估计法精度高的优点。经实验验证,该算法精度较好且运行效率较高。
  • 图 1  多基线干涉SAR系统几何关系示意图

    Figure 1.  Multi-baseline InSAR system geometric relationship

    图 2  改进的CABMAP算法流程框图

    Figure 2.  Optimized CABMAP algorithm flow chart

    图 3  仿真数据集

    Figure 3.  Simulation dataset

    图 4  算法结果

    Figure 4.  Result of several algorithms

    图 5  不同聚类次数下CABMAP算法结果

    Figure 5.  Result of CABMAP algorithm with different clustering iterations

    图 6  巴塞罗那实测数据集

    Figure 6.  Real dataset of Barcelona

    表 1  多基线干涉SAR系统仿真参数

    Table 1.  Multi-baseline InSAR system simulation parameters

    参数 数值
    景中心斜距(km) 500
    平台高度(km) 433
    视角(°) 30
    基线角(°) 5
    信号波长(cm) 3.1
    信号带宽(MHz) 100
    基线1长度(m) 199.794
    基线2长度(m) 133.196
    基线3长度(m) 79.918
    信噪比(dB) 30
    DEM网格间距(m) 1.5×1.5
    下载: 导出CSV

    表 2  算法性能对比(Isolation Peak仿真数据)

    Table 2.  Algorithm performance comparison (Simulation dataset: Isolation Peak)

    算法模型 时间(s) 精度(归一化均方误差)
    鲁棒性CRT 2.012318 1.2295
    CA 0.868110 0.9139
    CANOPUS 3.393945 0.8986
    ML估计 4.198716 3.8801
    ML估计后均值滤波 4.682451 2.0460
    MAP估计 229.851504 0.0117
    CABMAP 22.921332 0.0041
    改进的CABMAP 24.109097 0.0028
    下载: 导出CSV

    表 3  信噪比和精度关系

    Table 3.  Relationship between SNR and accuracy

    算法 信噪比(dB)
    5 10 15 20 25 30 35
    鲁棒性CRT 3.5594 3.1339 2.5860 2.0128 1.5663 1.2295 1.0121
    CA 1.0152 0.9763 0.9495 0.9290 0.9133 0.9139 0.9025
    CANOPUS 0.9311 0.9243 0.9177 0.9138 0.9019 0.8986 0.8919
    ML估计 4.8221 4.6486 4.4204 4.1489 4.0443 3.8801 3.7875
    ML估计后均值滤波 3.0314 2.7886 2.5026 2.2465 2.1602 2.0460 1.9831
    MAP估计 1.2184 0.8829 0.4422 0.1285 0.0229 0.0117 0.0111
    CABMAP 0.4145 0.2423 0.0768 0.0171 0.0064 0.0041 0.0030
    改进的CABMAP 0.4075 0.2254 0.0574 0.0110 0.0045 0.0028 0.0021
    下载: 导出CSV

    表 4  不同聚类迭代次数下CABMAP算法精度

    Table 4.  The accuracy of CABMAP algorithm for different clustering iterations

    聚类迭代次数N 精度(归一化均方误差)
    0 0.2378
    1 0.0081
    2 0.0041
    3 0.0168
    4 0.0240
    5 0.0299
    下载: 导出CSV

    表 5  不同阈值点数下CABMAP算法精度

    Table 5.  The accuracy of CABMAP algorithm for different threshold points

    阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 精度(归一化均方误差)
    2 0.0158
    3 0.0144
    4 0.0107
    5 0.0047
    6 0.0041
    7 0.0044
    下载: 导出CSV

    表 6  不同高程差下CABMAP算法精度

    Table 6.  The accuracy of CABMAP algorithm for different elevation difference

    高程差 $\Delta h$ (m) 精度(归一化均方误差)
    10 0.0042
    20 0.0041
    30 0.0043
    40 0.0043
    50 0.0043
    60 0.0045
    下载: 导出CSV

    表 7  实测数据集参数

    Table 7.  Real dataset parameters

    主图像获取时间
    (年-月-日)
    副图像获取时间
    (年-月-日)
    时间基线(d) 空间垂直基线(m)
    2009-02-20 2009-02-09 11 151.5
    2009-02-20 2009-01-29 11 86.6
    下载: 导出CSV

    表 8  算法运行时间

    Table 8.  Run time of the algorithms

    算法 时间(s)
    鲁棒性CRT 0.320831
    CA 0.518862
    CANOPUS 2.009810
    ML估计 0.161702
    ML估计后均值滤波 0.164832
    MAP估计 11.617357
    CABMAP 5.175075
    改进的CABMAP 5.253087
    下载: 导出CSV
  • [1] Bamler P and Hartl P. Synthetic aperture radar interferometry[J]. Inverse Problems, 1998, 14(4): R1–R54. DOI: 10.1088/0266-5611/14/4/001
    [2] Rosen P A, Hensley S, Joughin I R, et al. Synthetic aperture radar interferometry[J]. Proceedings of the IEEE, 2000, 88(3): 333–382. DOI: 10.1109/5.838084
    [3] Richards M A. A beginner’s guide to interferometric SAR concepts and signal processing [AESS Tutorial IV][J].IEEE Aerospace and Electronic Systems Magazine, 2007, 22(9): 5–29. DOI: 10.1109/MAES.2007.4350281
    [4] Karout S. Two-dimensional phase unwrapping[D]. [Ph. D. dissertation], Liverpool John Moores University, 2007.
    [5] 张妍, 冯大政, 曲小宁. 枝切法与曲面拟合结合的InSAR相位展开算法[J]. 西安电子科技大学学报(自然科学版), 2012, 39(5): 47–53. DOI: 10.3969/j.issn.1001-2400.2012.05.009Zhang Yan, Feng Dazheng, and Qu Xiaoning. Hybrid phase unwrapping algorithm combining branch-cut and surface-fitting for InSAR[J]. Journal of Xidian University(Natural Science), 2012, 39(5): 47–53. DOI: 10.3969/j.issn.1001-2400.2012.05.009
    [6] Krieger G, Hajnsek I, Papathanassiou K P, et al. Interferometric synthetic aperture radar (SAR) missions employing formation flying[J]. Proceedings of the IEEE, 2010, 98(5): 816–843. DOI: 10.1109/JPROC.2009.2038948
    [7] 庞蕾, 张继贤, 范洪冬. 多基线干涉SAR测量技术发展与趋势分析[J]. 电子学报, 2010, 38(9): 2152–2157Pang Lei, Zhang Jixian, and Fan Hongdong. Progress and tendency of multibaseline Synthetic Aperture Radar interferometry technique[J]. Acta Electronica Sinica, 2010, 38(9): 2152–2157
    [8] 侯丽英, 林赟, 洪文. 干涉圆迹SAR的目标三维重建方法研究[J]. 雷达学报, 2016, 5(6): 538–547. DOI: 10.12000/JR16009Hou Liying, Lin Yun, and Hong Wen. Three-dimensional reconstruction method study based on interferometric circular SAR[J]. Journal of Radars, 2016, 5(6): 538–547. DOI: 10.12000/JR16009
    [9] 师君, 张晓玲, 韦顺军, 等. 基于变分模型的阵列三维SAR最优DEM重建方法[J]. 雷达学报, 2015, 4(1): 20–28. DOI: 10.12000/JR14136Shi Jun, Zhang Xiao-ling, Wei Shun-jun, et al. An optimal DEM reconstruction method for linear array Synthetic Aperture Radar based on variational model[J]. Journal of Radars, 2015, 4(1): 20–28. DOI: 10.12000/JR14136
    [10] 靳国旺, 张红敏, 徐青, 等. 多波段InSAR的CRT相位解缠方法[J]. 西安电子科技大学学报(自然科学版), 2011, 38(6): 97–102. DOI: 10.3969/j.issn.1001-2400.2011.06.015Jin Guowang, Zhang Hongmin, Xu Qing, et al. Phase unwrapping algorithm with CRT for multi-band InSAR[J]. Journal of Xidian University(Natural Science), 2011, 38(6): 97–102. DOI: 10.3969/j.issn.1001-2400.2011.06.015
    [11] Yuan Zhihui, Deng Yunkai, Li Fei, et al. Multichannel InSAR DEM reconstruction through improved closed-form robust Chinese remainder theorem[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(6): 1314–1318. DOI: 10.1109/LGRS.2013.2238886
    [12] Yu Hanwen, Li Zhenfang, and Bao Zheng. A cluster-analysis-based efficient multibaseline phase-unwrapping algorithm[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(1): 478–487. DOI: 10.1109/TGRS.2010.2055569
    [13] Liu Huitao, Xing Mengdao, and Bao Zheng. A cluster-analysis-based noise-robust phase-unwrapping algorithm for multibaseline interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(1): 494–504. DOI: 10.1109/TGRS.2014.2324595
    [14] Pascazio V and Schirinzi G. Estimation of terrain elevation by multifrequency interferometric wide band SAR data[J]. IEEE Signal Processing Letters, 2001, 8(1): 7–9. DOI: 10.1109/97.889635
    [15] Pascazio V and Schirinzi G. Multifrequency InSAR height reconstruction through maximum likelihood estimation of local planes parameters[J]. IEEE Transactions on Image Processing, 2002, 11(12): 1478–1489. DOI: 10.1109/TIP.2002.804274
    [16] Ferraiuolo G, Pascazio V, and Schirinzi G. Maximum a posteriori estimation of height profiles in InSAR imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2004, 1(2): 66–70. DOI: 10.1109/LGRS.2003.822882
    [17] Li S Z. Modeling image analysis problems using markov random fields[J]. Handbook of Statistics, 2003, 21: 473–513. DOI: 10.1016/S0169-7161(03)21015-4
    [18] Ferraioli G, Shabou A, Tupin F, et al. Multichannel phase unwrapping with graph cuts[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(3): 562–566. DOI: 10.1109/LGRS.2009.2021165
    [19] Deledalle C A, Denis L, Ferraioli G, et al.. Combining patch-based estimation and total variation regularization for 3D InSAR reconstruction[C]. Proceedings of IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Milan, Italy, 2015: 2485–2488.
    [20] Ferraiuolo G, Meglio F, Pascazio V, et al. DEM reconstruction accuracy in multichannel SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 191–201. DOI: 10.1109/TGRS.2008.2002644
    [21] Saquib S S, Bouman C A, and Sauer K. ML parameter estimation for markov random fields with applications to bayesian tomography[J]. IEEE Transactions on Image Processing, 1998, 7(7): 1029–1044. DOI: 10.1109/83.701163
  • [1] 黄瑞杜小勇胡卫东 . OFDM雷达多目标运动参数的近似最大似然估计. 雷达学报, 2018, 7(4): 507-513. doi: 10.12000/JR17116
    [2] 李强范怀涛 . 基于辅助数字高程模型的方位多通道SAR相位失配校正方法. 雷达学报, 2019, (): 1-8. doi: 10.12000/JR19009
    [3] 刘华有郑明洁张衡王宇秦小芳 . 一种有效的机载双频干涉地形高程重建方法. 雷达学报, 2018, 7(4): 475-486. doi: 10.12000/JR18013
    [4] 秦小芳张华春张衡王宇刘华有 . TerraSAR-X/TanDEM-X升降轨双基干涉模式获取DEM方法研究. 雷达学报, 2018, 7(4): 487-497. doi: 10.12000/JR18022
    [5] 张斌韦立登胡庆荣李爽 . 基于四阶累积量的机载多基线SAR谱估计解叠掩方法. 雷达学报, 2018, 7(6): 740-749. doi: 10.12000/JR18087
    [6] 刘忠胜汪丙南向茂生陈龙永 . 柔性基线抖动对机载干涉SAR 性能影响分析. 雷达学报, 2014, 3(2): 183-191. doi: 10.3724/SP.J.1300.2014.14015
    [7] 匡辉杨威王鹏波陈杰 . 多方位角多基线星载SAR三维成像方法研究. 雷达学报, 2018, 7(6): 685-695. doi: 10.12000/JR18073
    [8] 赵军香梁兴东李焱磊 . 一种基于似然比统计量的SAR相干变化检测. 雷达学报, 2017, 6(2): 186-194. doi: 10.12000/JR16065
    [9] 朱莹张弓张劲东 . 基于DCS 的统计MIMO 雷达信号模型及参数估计. 雷达学报, 2012, 1(2): 143-148. doi: 10.3724/SP.J.1300.2012.20016
    [10] 乔明潘舟浩刘波李道京 . 毫米波三基线InSAR 通道泄漏误差分析和补偿方法研究. 雷达学报, 2013, 2(1): 68-76. doi: 10.3724/SP.J.1300.2013.13008
    [11] 李学仕孙光才邵鹏吴玉峰邢孟道 . 基于数字阵列雷达的同时多模式SAR 成像体制研究. 雷达学报, 2014, 3(4): 480-489. doi: 10.3724/SP.J.1300.2014.13113
    [12] 刘玉琪易建新万显荣程丰饶云华龚子平 . 数字电视外辐射源雷达多旋翼无人机微多普勒效应实验研究. 雷达学报, 2018, 7(5): 585-592. doi: 10.12000/JR18062
    [13] 杨海峰谢文冲王永良 . 不同发射波形下多机协同探测雷达系统信号建模与性能分析. 雷达学报, 2017, 6(3): 267-274. doi: 10.12000/JR16142
    [14] 吴琴霞梁兴东李焱磊周良将 . 残余运动误差对机载多波段SAR 图像配准的影响分析. 雷达学报, 2015, 4(2): 209-216. doi: 10.12000/JR14065
    [15] 谭云华王李波李廉林 . 一种抑制探地/穿墙成像多径虚假目标的新型概率模型:数值研究(英文). 雷达学报, 2015, 4(5): 509-517. doi: 10.12000/JR15066
    [16] 张强万显荣傅饶云华龚子平 . 基于CDR 数字音频广播的外辐射源雷达信号模糊函数分析与处理. 雷达学报, 2014, 3(6): 702-710. doi: 10.12000/JR14050
    [17] 朱洪伟唐小明何友 . 基于航迹无关模型的传感器系统误差可观测性分析. 雷达学报, 2013, 2(4): 454-460. doi: 10.3724/SP.J.1300.2013.13068
    [18] 浮瑶瑶柳彬张增辉郁文贤 . 基于词包模型的高分辨率SAR 图像变化检测与分析. 雷达学报, 2014, 3(1): 101-110. doi: 10.3724/SP.J.1300.2014.13134
    [19] 黄传禄晁坤毛云志 . 空间谱估计中误差自校正方法研究. 雷达学报, 2014, 3(5): 518-523. doi: 10.3724/SP.J.1300.2014.13147
    [20] 陈功伯李勇陶满意 . 基于信号的环视SAR 成像参数估计方法. 雷达学报, 2013, 2(2): 203-209. doi: 10.3724/SP.J.1300.2013.20073
  • 加载中
图(7)表(8)
计量
  • 文章访问数:  442
  • HTML浏览量:  44
  • PDF下载量:  406
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-04-01
  • 录用日期:  2017-04-25
  • 刊出日期:  2017-12-28

一种基于最大后验框架的聚类分析多基线干涉SAR高度重建算法

    通讯作者: 斯奇, si_qi0616@163.com
    作者简介: 王   宇(1980–),男,河南人,现为中国科学院电子学研究所研究员,博士生导师,研究方向为SAR系统设计与信号处理技术。E-mail: yuwang@mail.ie.ac.cn;邓云凯(1962–),男,湖北人,现为中国科学院电子学研究所研究员,博士生导师,研究方向为星载合成孔径雷达系统设计。E-mail: ykdeng@mail.ie.ac.cn
  • ①. 中国科学院电子学研究所   北京   100190
  • ②. 中国科学院大学   北京   100039
基金项目:  国家自然科学基金优秀青年基金(61422113),国家万人计划青年拔尖人才,中科院百人计划

摘要: 多基线干涉SAR能有效减小由目标高度急剧变化和较大噪声干扰带来的不利影响,可以获取比单基线干涉SAR更精确的地表数字高程模型(DEM)。传统的基于最大似然估计(ML)的多基线高度重建算法在通道数目较少情况下重建结果不佳,基于最大后验估计(MAP)的多基线高度重建算法存在运行时间较长的问题,针对以上问题,该文提出了一种基于最大后验框架的聚类分析高度重建算法(CABMAP)。该算法首先利用了ML估计法得到粗略的DEM,以此为基础在每次迭代过程中利用聚类分析(CA)判断出邻域内的噪声像素,并通过计算后验概率完成重建,此外采用了一种改进措施提高精度。这样,既保留了ML估计法运行速度快的特征,又具有MAP估计法精度高的优点。经实验验证,该算法精度较好且运行效率较高。

English Abstract

    • 干涉合成孔径雷达(Interferometric Synthetic Aperture Radar, InSAR)可以利用两幅SAR图像的复数据提取出相位信息以获取地表的高程信息,是对传统SAR成像技术的重要拓展[13]。传统单基线InSAR技术利用相位解缠的方法求得干涉相位后,根据轨道参数及成像几何关系计算出相应的地形高程,即数字高程模型(Digital Elevation Model, DEM)。但是该方法建立在伊藤假设的基础上[4],即假定任意相邻两点的干涉相位差的绝对值小于 ${π} $ ,其本质是要求被测绘地区具有空间连续性,但在复杂场景区域(高山、城市等),该假设通常不成立,此时相位解缠算法很难获得令人满意的结果[5]。为解决这一问题,学者们提出了多基线干涉SAR技术,该技术可以有效地解决复杂地形的相位缠绕问题,从而提升高度重建精度[69]

      目前,有很多利用多基线干涉SAR数据进行高度重建的技术得到了发展。理论上,当基线长度满足互质条件时,中国余数定理(Chinese Reminder Theorem, CRT)可被用于相位解缠[10],但是该算法噪声鲁棒性较差,实际应用中效果不佳。为解决这一问题,袁志辉等人提出了一种闭式解的鲁棒性CRT算法[11],得到了更高精度的DEM。通过研究多基线数据特有的聚类特性,于翰雯等人提出了一种基于聚类分析(Cluster-Analysis, CA)的多基线干涉SAR相位解缠算法[12],该算法运算速度较快但易受噪声干扰。为了提升CA算法的噪声鲁棒性,保铮等人提出了一种基于密度聚类的噪声鲁棒性相位解缠算法(Cluster-Analysis based Noise-Robust Phase-Unwrapping algorithm, CANOPUS)[13]。近年来,通过求取概率极值完成高度重建的算法被广泛地研究,这些算法主要分为两种:较早的一种是基于最大似然(Maximum Likelihood, ML)估计的算法[14],这类算法可以在较短时间内重建出高度,但是其性能较依赖于干涉图的数目、工作频率、基线值和相干系数[15];另一类是基于最大后验(Maximum A Posterior, MAP)估计的算法[16],这类算法利用马尔可夫随机场统计分布模型来描述目标高度的先验分布,通过估计每个像素对应的超参数实现目标高度的估计。相比ML估计法,MAP估计法能够较精确地估计出高度,但是会产生运算复杂度变高、内存增加的问题[17]。随后,又出现了一些MAP估计的改进算法[18,19]。虽然基于MAP估计的方法比基于ML估计的方法精度高,但由于该方法需要重复迭代估计每个像素对应的超参数和高度值,因而存在运行时间长的缺陷,基于ML估计的方法则不存在这一缺点。

      针对ML估计法精度不高、MAP估计法速度较慢的问题,本文提出了一种基于最大后验框架的聚类分析(Cluster-Analysis Based on Maximum A Posterior, CABMAP)高度重建算法。该算法首先利用ML估计法得到粗略的DEM后开始迭代,利用聚类分析判断噪声像素和非噪声像素后进行超参数估计,接着通过计算后验概率更新目标像素点的高度并以此为基础进行下一次迭代,重复上述步骤,直到满足阈值条件后终止算法得到最终重建结果。和ML估计法相比,CABMAP算法考虑了目标像素的先验分布,具有MAP估计法精度高的优势,并通过聚类分析区分了噪声像素和非噪声像素,比MAP估计法进一步地提升了重建精度;和MAP估计法相比,CABMAP算法以ML估计法得到的粗略DEM为基础,减少了MAP估计法重复迭代估计的次数,因此具有ML估计法运行速度快的优势。本文第2节简要介绍了基于ML估计的高度重建算法和基于MAP估计的高度重建算法,第3节详述CABMAP及其改进措施,第4节将分别给出仿真和实测实验的结果和分析,最后一节给出结论。

    • 本文提出的高度重建算法结合了ML估计法和MAP估计法的优点,因此在详细描述本文提出的算法之前,先简要介绍ML估计法和MAP估计法。

    • 多基线干涉SAR系统的几何示意图如图1所示,图中, ${A_0},{A_1}, ·\!·\!· ,{A_N}$ 为垂直于航迹方向的天线相位中心, $\theta $ 为主天线的侧视角, ${r_0},{r_1}, ·\!·\!· ,{r_N}$ 分别为地面某点P ${A_0},{A_1}, ·\!·\!· ,{A_N}$ 的斜距, ${y_p}$ 是地距,hP点的高度值, ${A_1}, ·\!·\!· ,{A_N}$ ${A_0}$ 分别形成稳定的干涉基线 ${B_1}, ·\!·\!· ,{B_N}$ ,对应的水平基线角分别为 ${\alpha _1},{\alpha _2}, ·\!·\!· ,{\alpha _N}$

      图  1  多基线干涉SAR系统几何关系示意图

      Figure 1.  Multi-baseline InSAR system geometric relationship

      在多基线干涉SAR系统中,令 ${\varphi _l}(p)$ 为干涉图中像素pl条基线对应的缠绕相位 $( - {π} \le {\phi _l}(p) \le {π} )$ , ${\phi _l}(p)$ 为像素pl条基线对应的真实相位,那么缠绕相位和真实相位满足:

      ${\phi _l}(p) = 2{π} {K_l} + {\varphi _l}(p),\ \ \{ {K_l}|{K_l} \in N\} $

      在基于ML估计的高度重建算法中,高程重建问题被转化成了一个寻找使似然函数最大化的高程值h(p)的问题,该似然函数可写成:

      $ \hat h(p) = \mathop {\arg \max }\limits_{{h}} F\left( {\varPhi \left( p \right)|h\left( p \right)} \right)$

      其中, $\hat h(p)$ 是像素p的估计高度, $F\left( {\varPhi \left( p \right)|h\left( p \right)} \right)$ 为观测相位 $\varPhi \left( p \right)$ 关于真实高度 $h\left( p \right)$ 的概率密度。对多基线干涉SAR系统, $F\left( {\varPhi \left( p \right)|h\left( p \right)} \right)$ 可以近似看成K个独立的缠绕干涉相位观测值的乘积[14]

      $F\left( {\varPhi \left( p \right)|h\left( p \right)} \right) = \prod\limits_{k = 1}^K {f({\phi _k}(p)|h(p))} $

      对于K=1的特殊情况,式(3)可以写成[14]

      $\begin{aligned} f(\phi (p)|h(p)) = & \frac{1}{{2{π} }}\frac{{1 - |\gamma {|^2}}}{{1 - |\gamma {|^2}{{\cos }^2}\left(\phi (p) - \displaystyle\frac{{4{π} {B_ \bot }}}{{\lambda {R_0}\sin \theta }}h(p)\right)}}\\ & \cdot \left\{ {1 + \frac{{|\gamma |\cos \left(\phi (p) - \displaystyle\frac{{4{π} {B_ \bot }}}{{\lambda {R_0}\sin \theta }}h(p)\right){{\cos }^{ - 1}}\left[ - |\gamma |\cos \left(\phi (p) - \displaystyle\frac{{4{π} {B_ \bot }}}{{\lambda {R_0}\sin \theta }}h(p)\right)\right]}}{{{{\left[1 - |\gamma {|^2}{{\cos }^2}\left(\phi (p) - \displaystyle\frac{{4{π} {B_ \bot }}}{{\lambda {R_0}\sin \theta }}h(p)\right)\right]}^{1/2}}}}} \right\} \end{aligned}$

      需注意,基于ML估计的重建算法速度较快,但是抗噪性能较差,且易受工作频率、基线数目和相干系数的影响[20]

    • 基于MAP估计的高度重建算法认为目标高度是一个随机变量,并可用一个先验分布模型进行建模。通过引入关于地形高度的先验分布模型,寻找使后验概率最大化的高度 $h(p)$ ,从而达到高度重建的目的。假设目标区域的高度为一幅图像,所有像素的高度组成随机矢量H,而该随机矢量的一个具体实现形式为 ${h} = {[h(1),h(2),·\!·\!·,h(N \times M)]^{\rm T}}$ ,即其中的一个样本。如果对随机矢量H采用某种先验统计模型进行建模,且使某种样本h出现的概率达到最大,就可以使用最大后验估计的方法进行高程重建。马尔可夫随机场模型可以描述该矢量集合的分布,这种模型考虑了每个像素高度关于它的一组邻近像素高度的条件分布,可以有效地描述图像的局部特征,任何类型的马尔可夫随机场均可以用吉布斯分布的形式来表示[17]

      $g({h};{σ} ) = \frac{1}{{Z({σ} )}}\exp \{ - \varepsilon ({h};{σ} )\} $

      其中, $Z({σ} ) = \displaystyle\int {\exp \{ - \varepsilon ({h};{σ} )\} {\rm d}{h}} $ 是一个归一化因子,称为划分函数, $\varepsilon ( \cdot )$ 描述了相邻像素间的关系,称为能量函数, ${σ} $ 是未知参数的集合,称为超参数矢量。选取高斯马尔可夫随机场模型描述高程图像的特征[16,21]

      $\varepsilon ({h};{σ} ) = \sum\limits_{p = 1}^{N \times M} {\sum\limits_{j \in {C_p}} {\frac{{{{(h(p) - h(j))}^2}}}{{2\sigma _{pj}^2}}} } $

      其中, ${C_p}$ 为像素p对应的邻域系统[17],一般采用其周围的8个像素即可,而 ${\sigma _{pj}}$ 是超参数矢量 ${σ} $ 中的元素,描述了高度图像的局部空间变化特征。根据贝叶斯定律,多基线干涉SAR的后验概率密度函数可以写成[21]

      $\begin{aligned} f({h}|{Φ} ) = & F({Φ} |{h}) \cdot g({h};{σ} )\\ = & \biggr\{ \prod\limits_{p = 1}^{N \times M} {\prod\limits_{k = 1}^K {f({\phi _k}(p)|h(p))\biggr\} \frac{1}{{Z({σ} )}}} } \\ & \cdot \exp \{ - \varepsilon ({h};{σ} )\} \end{aligned}$

      其中, $F({Φ} |{h})$ 是最大似然项。可将高程重建问题转换成后验概率估计问题:

      $ \begin{array}{l} {\hat {h}} \!\! & = \mathop {\arg \max }\limits_{{h} \in {H}} f({h}|{Φ} ) = \mathop {\arg \max }\limits_{{h} \in {H}} [\ln f({{h}}|{Φ} )]\\ & \!\!\!\! = \mathop {\arg \max }\limits_{{h} \in {H}} \left\{ {\sum\limits_{p = 1}^{N \times M} {\sum\limits_{k = 1}^K {\ln f({\phi _k}(p)|h(p))} } } \right.\\ & \!\!\!\!\!\!\!\!\!\!\!\!\! \left. { - \sum\limits_{p = 1}^{N \times M} {\sum\limits_{j \in {C_p}} {\frac{{{{(h(p) - h(j))}^2}}}{{2\sigma _{pj}^2}}} } } \right\} \end{array} $

      有多种方法可以求解式(8),如模拟退火算法[17],也可随机抽样后使用条件迭代模(Iterative Conditional Modes, ICM)算法进行求解[16]。相比基于ML估计的方法,基于MAP估计的方法重建精度更高[20],但该方法需要重复迭代估计每个像素对应的超参数和高程值,因此算法速度较慢。

    • 基于ML估计的方法运行速度较快但精度较低,而基于MAP估计的方法精度较高、运行速度较慢。结合两种方法的优点,本文提出了CABMAP算法,该算法将ML估计法得到的粗略的地形高度作为初始高度,这样可以减少MAP估计法中重复迭代估计每个像素对应的高程值和超参数的次数,节约了算法的运行时间。将ML估计法得到的结果作为初始高度进行迭代,每次迭代中,利用聚类分析方法判断出噪声像素后利用MAP框架计算后验概率,取最大后验概率对应的高程值为该次迭代的结果,并以该结果为基础进行下一次迭代,直到迭代次数满足设定的阈值后结束算法,得到最终的高度重建结果。

      基于MAP估计的方法利用了目标像素的邻域内所有的像素辅助重建[16],但是当其邻域中存在被噪声严重污染的像素时,目标像素的重建会受到严重干扰。为解决该问题,CABMAP算法采用了聚类的思想,对目标像素的邻域进行聚类处理,判断并剔除目标像素邻域内的噪声像素,有效地防止了噪声像素对结果的干扰,提高了重建精度。算法的具体步骤如下:

      步骤1  设置算法的结束条件,当迭代次数达到N时结束;设置邻域内能聚成一类的高程差阈值为 $\Delta h$ ,邻域中能聚为一类的像素点个数阈值为 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 。将第i次迭代中的重建高度记为 ${{h}^{(i)}}$ ,超参数记为 ${{σ} ^{(i)}}$ 。首先利用ML估计出粗略的地形高度 ${{h}^{(0)}}$

      步骤2  对目标像素p,判断是否为噪声像素,若与其邻域 ${C_p}$ 中某个像素q的高程差不超过 $\Delta h$ ,则p, q判为一类。考查邻域中同类像素点的个数,如果不小于 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ ,则像素p不是噪声像素,反之为噪声像素。下式中, ${\varOmega _{{\rm{noise}}}}$ 表示噪声像素的集合:

      $\left\{ \begin{array}{l} \sum\limits_{q \in {C_p}} {f({h_p} - {h_q}} ) < {\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}, \ {\rm{ }}p \in {\varOmega _{{\rm{noise}}}}\\ \sum\limits_{q \in {C_p}} {f({h_p} - {h_q}} ) \ge {\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}, \ {\rm{ }}p \notin {\varOmega _{{\rm{noise}}}} \end{array} \right.$

      $f( \cdot )$ 为:

      $f(x) = \left\{ \begin{array}{l}\!\! 1, \ {\rm{ }}x \le \Delta h\\ \!\! 0, \ {\rm{ }}x > \Delta h \end{array} \right.$

      其中, ${\varOmega _{{\rm{noise}}}}$ 表示噪声像素的集合, $f({h_p} - {h_q}) = 1$ 表示像素p和像素q对应的高度差不超过 $\Delta h$ ,此时像素p和像素q判为一类, $f({h_p} - {h_q}) = 0$ 表示像素p和像素q对应的高度差超过 $\Delta h$ ,此时像素p和像素q不属于同一类。 $\sum\limits_{q \in {C_p}} {f({h_p} - {h_q}} )$ 表示目标像素p邻域内和p属于同一类的像素点个数,比较该值和 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ ,若小于 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ ,目标像素p则被判为噪声像素,否则为非噪声像素。

      步骤3  对第i次迭代,如果像素p不是噪声像素,超参数 $\sigma {(p)^{(i)}}$ 可通过下式来估计:

      $\sigma {(p)^{(i)}} = \frac{{\displaystyle\sum\limits_{j \in {C_p}} {{{\left(h{{(p)}^{(i)}} - h{{(j)}^{(i)}}\right)}^2}} }}{{|{{\tilde C}_p}|}}$

      反之,如果像素p是噪声,超参数 $\hat \sigma {(p)^{(i)}}$ 可表示为:

      $\sigma {(p)^{(i)}} = \frac{{\displaystyle\sum\limits_{j \in {C_p}} {{{\left(h{{(p)}^{(i)}} - h{{(j)}^{(i)}}\right)}^2}} }}{{|{C_p}|}}$

      其中, $|{\tilde C_p}|$ 是非噪声目标像素邻域内同类像素的个数, $|{C_p}|$ 是噪声像素邻域内像素的个数。基于MAP的估计法中,估计超参数时没有区分噪声像素和非噪声像素,即在目标像素邻域内统一采用式(12)进行估计,这种做法会使被噪声严重污染的像素干扰目标像素的重建。CABMAP算法首先利用步骤2区分了噪声像素和非噪声像素,针对非噪声的目标像素,利用其邻域内的同类像素进行超参数估计,这样避免了噪声像素的干扰,而对噪声像素,则保留传统MAP估计法的处理方式,即采用式(12)进行超参数估计。

      步骤4  估计出第i次迭代的高程。如果像素p是噪声像素,使用下式进行估计:

      $ h{(p)^{(i + 1)}} = \mathop {\arg \max }\limits_{{h} \in {H}} \left\{ {\ln F \left({Φ} (p)|h{{(p)}^{(i + 1)}}\right) - {{\sum\limits_{j \in {C_p}} {\left[ {\frac{{\left(h{{(p)}^{(i + 1)}} - h{{(j)}^{(i + 1)}}\right)}}{{\sqrt 2 \left(\displaystyle\frac{{\sigma {{(p)}^{(i + 1)}} + \sigma {{(j)}^{(i + 1)}}}}{2}\right)}}} \right]} }^2}} \right\} $

      如果像素p不是噪声像素,则采用:

      $h{(p)^{(i + 1)}} = \mathop {\arg \max }\limits_{{h} \in {H}} \left\{ {\ln {F_{MCh}}\left({Φ} (p)|h{{(p)}^{(i + 1)}}\right) - {{\sum\limits_{j \in {{\tilde C}_p}} {\left[ {\frac{{(h{{(p)}^{(i + 1)}} - h{{(j)}^{(i + 1)}}}}{{\sqrt 2 \left(\displaystyle\frac{{\sigma {{(p)}^{(i + 1)}} + \sigma {{(j)}^{(i + 1)}}}}{2}\right)}}} \right]} }^2}} \right\}$

      对噪声像素,式(13)中邻域系统是 ${C_p}$ ,由目标像素的邻域像素构成;对非噪声像素,式(14)中邻域系统是 ${\tilde C_p}$ ,由 ${C_p}$ 中和目标像素同类的像素组成。在此步中,利用了噪声像素和非噪声像素的超参数,并根据噪声像素和非噪声像素不同的邻域特征进行高度重建,分别求式(13)、式(14)的极值从而得到此次迭代的高度重建结果。

      步骤5  迭代次数加1,判断是否达到迭代阈值,如果达到结束算法,否则返回步骤2继续执行。

      为了考查CABMAP算法的运行效率,需要量化分析该算法的时间复杂度。考虑在ML估计法中,令 ${{H}} = \{ {h_1},{h_2}, ·\!·\!· ,{h_L}\} $ 表示像素对应的候选高度值,干涉图大小为 ${\rm{row}} \times {\rm{col}}$ ,则ML估计法的时间复杂度为 $O(L \times {\rm{row}} \times {\rm{col}})$ , CABMAP算法ML估计法的结果上,进行聚类处理后计算后验概率,相比ML估计法,CABMAP算法增加了聚类步骤,令聚类迭代次数为N,则算法时间复杂度为 $O(L \times {\rm{row}} \times {\rm{col}} + N \times L \times {\rm{row}} \times {\rm{col}})$ ,当N取较小时,CABMAP算法运行时间和ML估计法近似。

    • 在CABMAP算法中,噪声像素的邻域中每个像素都对其重建有贡献,这种措施利用了噪声像素邻域中所有的像素进行均衡地重建,但却隐含着下列问题:当噪声像素的缠绕相位被严重污染时,后验概率中最大似然项的准确性会大大降低,此时即使利用了噪声像素的邻域特征作为先验项进行弥补,也无法准确地重建高度。为了解决该问题,本文在CABMAP算法的基础上加入了一种改进措施。图2是改进的CABMAP算法的流程图,其中虚线左边是CABMAP算法的流程,虚线右边是改进措施的流程。

      图  2  改进的CABMAP算法流程框图

      Figure 2.  Optimized CABMAP algorithm flow chart

      具体的改进措施如下:

      步骤1  设置算法的终止条件,当迭代次数达到M时结束。将CABMAP算法得到的高度作为初始的高度 ${{h}^{(0)}}$

      步骤2  对第i次迭代,判断像素p是否为噪声像素,具体判断方法仍如式(9)。若p是噪声像素,则采用均值滤波的方式对其进行高度重建,若p不是噪声像素,则保持高度不变。

      ${h_p}^{\!\!(i)} = \left\{ \begin{array}{l}\!\! \frac{1}{{|{{\tilde C}_p}|}}\sum\limits_{j \in {{\tilde C}_p}} {h_j^{(i - 1)}} , \ {\rm{ }}p \in {\varOmega _{{\rm{noise}}}}\\ h_p^{(i - 1)}, \quad\quad\quad \ \, {\rm{ }}p \notin {\varOmega _{{\rm{noise}}}} \end{array} \right.$

      其中, ${\tilde C_p}$ 表示邻域中能够聚成一类的像素, ${\varOmega _{{\rm{noise}}}}$ 为噪声像素的集合。 $h_p^{(i)}$ 为像素pi次迭代的高度。

      步骤3  迭代次数加1,判断是否达到阈值,如果达到则结束,否则转到步骤2继续执行。

    • 本部分采用了仿真数据和实测数据进行实验,从算法运行时间、精度两个方面将本文提出的算法与其它几种算法进行对比,验证了该算法的可行性。另外,本文所有算法的实现均采用Matlab编程,并在内存为4 GB, CPU主频为3.2 GHz的主机上运行。

    • 本文采用了一组接近真实地形的仿真数据进行实验,该仿真数据是根据美国Isolation Peak国家公园的真实数字高程得到的。仿真参数如表1所示,3条基线对应的高度模糊数分别为21.4 m, 32.1 m和53.5 m。图3(a)图3(b)图3(c)分别为仿真得到的3条不同基线对应的缠绕干涉相位图。图3(d)为Isolation Peak国家公园的真实DEM(458×157像素)。图3(d)右侧有一条垂直高度差达到137 m的悬崖,对应的干涉相位差分别为12.8 ${π} $ , 8.5 ${π} $ 和5.1 ${π} $ ,不满足相位连续性假设,不能使用传统的单基线相位解缠的方法对整个区域进行高程重建。

      参数 数值
      景中心斜距(km) 500
      平台高度(km) 433
      视角(°) 30
      基线角(°) 5
      信号波长(cm) 3.1
      信号带宽(MHz) 100
      基线1长度(m) 199.794
      基线2长度(m) 133.196
      基线3长度(m) 79.918
      信噪比(dB) 30
      DEM网格间距(m) 1.5×1.5

      表 1  多基线干涉SAR系统仿真参数

      Table 1.  Multi-baseline InSAR system simulation parameters

      图  3  仿真数据集

      Figure 3.  Simulation dataset

      图4(a)图4(h)分别是鲁棒性CRT, CA, CANOPUS, ML, ML估计后均值滤波、MAP估计、CABMAP及改进措施的高度重建结果。下面分别从精度和运行时间两个方面对上述算法的性能进行分析。采用归一化的均方误差来衡量精度,定义为: $\bar \varepsilon = \parallel {\hat {h}} - {h} {\parallel^2}/\parallel{h}{\parallel^2}$ ,其中h为各个像素的真实高度所组成的矩阵, ${\hat {h}}$ 为对应的估计值, ${\left\| \cdot \right\|^{\rm{2}}}$ 表示取范数的平方。在30 dB的信噪比水平下,这8种方法的归一化均方误差分别为1.2295, 0.9139, 0.8986, 3.8801, 2.0460, 0.0117, 0.0041, 0.0028,可以发现,采用本文方法得到的精度效果最好。其次,从运行时间来看,CA算法的速度最快,只要0.868110 s, MAP估计法则最慢,需要229.851504 s, CABMAP算法运行时间是22.921332 s,由于减少了重复估计每个像素的超参数次数,CABMAP算法运行速度快于MAP估计法,各算法的运行时间和精度具体见表2。从以上两个方面的对比分析,本文的算法在较短的时间内大大提高了重建精度,适用于多基线数据的处理。

      算法模型 时间(s) 精度(归一化均方误差)
      鲁棒性CRT 2.012318 1.2295
      CA 0.868110 0.9139
      CANOPUS 3.393945 0.8986
      ML估计 4.198716 3.8801
      ML估计后均值滤波 4.682451 2.0460
      MAP估计 229.851504 0.0117
      CABMAP 22.921332 0.0041
      改进的CABMAP 24.109097 0.0028

      表 2  算法性能对比(Isolation Peak仿真数据)

      Table 2.  Algorithm performance comparison (Simulation dataset: Isolation Peak)

      图  4  算法结果

      Figure 4.  Result of several algorithms

      为了探究噪声对重建精度的影响,分别对仿真数据加入不同程度的噪声后进行实验。仍然采用归一化均方误差作为精度指标。从表3中,在信噪比相同时,改进的CABMAP算法精度最高,其次是CABMAP算法。随着信噪比的增加,各算法的精度均会增加,对比各类算法的增速,CABMAP算法及改进措施的精度增速较快。综合以上两个方面,本文的算法可以在不同噪声水平的数据上均表现出良好的性能。

      算法 信噪比(dB)
      5 10 15 20 25 30 35
      鲁棒性CRT 3.5594 3.1339 2.5860 2.0128 1.5663 1.2295 1.0121
      CA 1.0152 0.9763 0.9495 0.9290 0.9133 0.9139 0.9025
      CANOPUS 0.9311 0.9243 0.9177 0.9138 0.9019 0.8986 0.8919
      ML估计 4.8221 4.6486 4.4204 4.1489 4.0443 3.8801 3.7875
      ML估计后均值滤波 3.0314 2.7886 2.5026 2.2465 2.1602 2.0460 1.9831
      MAP估计 1.2184 0.8829 0.4422 0.1285 0.0229 0.0117 0.0111
      CABMAP 0.4145 0.2423 0.0768 0.0171 0.0064 0.0041 0.0030
      改进的CABMAP 0.4075 0.2254 0.0574 0.0110 0.0045 0.0028 0.0021

      表 3  信噪比和精度关系

      Table 3.  Relationship between SNR and accuracy

      本文提出的CABMAP算法需要设置迭代次数N、高程差 $\Delta h$ 和阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 3个参数,这些参数将会影响算法的高度重建结果,为了探究各参数对重建精度的影响,下面给出不同参数组合下的重建结果。

      首先研究聚类迭代次数N对精度的影响,采用八邻域结构的邻域系统,分别固定阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}} =$ 6和高程差 $\Delta h = 20\;{\rm{m}}$ ,变化迭代次数N图5(a)图5(f)分别是不经过聚类处理及聚类迭代次数分别是1~5时的重建结果。在不同的聚类迭代次数下,CABMAP的精度如表4所示。由表2, ML估计法的精度为3.8801, MAP估计法的精度为0.0117,由表4,如果不进行聚类处理,直接以ML估计法得到的结果为基础,采用MAP框架重建得到的精度为0.2378,即其精度虽然比ML估计法高,但低于MAP估计法。而通过聚类步骤判断出邻域内的噪声像素后再采用MAP框架进行重建,实验结果表明精度得到了显著提升。同时,随着迭代次数的增加,精度呈现先增后降的趋势,当N=2时可以取得最好的重建结果。

      图  5  不同聚类次数下CABMAP算法结果

      Figure 5.  Result of CABMAP algorithm with different clustering iterations

      聚类迭代次数N 精度(归一化均方误差)
      0 0.2378
      1 0.0081
      2 0.0041
      3 0.0168
      4 0.0240
      5 0.0299

      表 4  不同聚类迭代次数下CABMAP算法精度

      Table 4.  The accuracy of CABMAP algorithm for different clustering iterations

      为了考查阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 对重建精度的影响,固定迭代次数N=2和高程差 $\Delta h = 20\;{\rm{m}}$ 。从表5中,归一化均方误差会随着阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 的增加呈现先降后增的趋势,当 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}} = 6$ 时可以取得最好的重建结果。

      阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 精度(归一化均方误差)
      2 0.0158
      3 0.0144
      4 0.0107
      5 0.0047
      6 0.0041
      7 0.0044

      表 5  不同阈值点数下CABMAP算法精度

      Table 5.  The accuracy of CABMAP algorithm for different threshold points

      固定迭代次数N=2和阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}} = 6$ ,考查高程差 $\Delta h$ 对重建精度的影响。从表6中,随着高程差 $\Delta h$ 的变化,归一化均方误差呈现先减小后增大的趋势,当高程差 $\Delta h = 20\;{\rm{m}}$ 时,可以取得最好的重建结果。

      高程差 $\Delta h$ (m) 精度(归一化均方误差)
      10 0.0042
      20 0.0041
      30 0.0043
      40 0.0043
      50 0.0043
      60 0.0045

      表 6  不同高程差下CABMAP算法精度

      Table 6.  The accuracy of CABMAP algorithm for different elevation difference

      通过不同参数情况下精度的量化分析,聚类迭代次数N对重建结果的影响较明显,当N偏大时,时间消耗变大且精度下降,当N偏小时会导致重建结果未收敛,从而精度下降,从实验结果来看,选择较小的迭代次数可以取得良好的重建结果,一般将迭代次数设为2~4。阈值点数 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 和高程差 $\Delta h$ 同样也会影响重建精度,当 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 设置偏大时,目标像素被判为噪声像素的概率变大,存在将非噪声像素错判成噪声像素的问题,当 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ 设置偏小时,存在将噪声像素错判成非噪声像素的问题;当 $\Delta h$ 设置偏大,则目标像素邻域内所有的像素被判为聚成一类,无法剔除邻域内的噪声像素,如果 $\Delta h$ 设置偏小,则邻域内的像素都被判断为噪声像素,也会使重建精度下降;实际中不同的噪声水平和数据特征均会影响 ${\rm{Hpt}}{{\rm{s}}_{{\rm{th}}}}$ $\Delta h$ 的设置,为了得到较好的重建结果,需要给不同的参数设定步长,考查不同的参数组合从而得到较好的重建结果。

    • 实测实验采用的数据是西班牙巴塞罗那地区某建筑物场景的真实干涉SAR图像。该数据集有两幅干涉相位图,由TerraSAR-X卫星以重复轨道的方式获得。两条基线对应的高度模糊数分别为36.64 m, 64.12 m。在形成干涉条纹前,对SAR图像做了多视处理,距离向5视,方位向4视。参数如表7

      主图像获取时间
      (年-月-日)
      副图像获取时间
      (年-月-日)
      时间基线(d) 空间垂直基线(m)
      2009-02-20 2009-02-09 11 151.5
      2009-02-20 2009-01-29 11 86.6

      表 7  实测数据集参数

      Table 7.  Real dataset parameters

      图6(a)图6(b)分别是基线1和基线2的干涉相位图,由于已经去除平地效应,所以看不见明显的干涉条纹,图6(c)图6(d)分别是基线1和基线2的相干系数图,可见大部分区域的相干性比较良好,因此适合做干涉处理。

      图  6  巴塞罗那实测数据集

      Figure 6.  Real dataset of Barcelona

      图7(a)图7(h)分别是鲁棒的CRT, CA, CANOPUS, ML, ML估计后均值滤波、MAP, CABMAP及改进措施的高度重建结果,图7(i)是改进的CABMAP重建结果的3维图,图7(j)是场景的光学图。首先从运行时间方面来分析,CABMAP和MAP估计法的运行时间分别是5.175075 s, 11.617357 s,相比MAP估计法,CABMAP具有较快的处理效率,具体见表8。从精度方面来分析,对比场景的光学图,鲁棒性CRT的结果中场景下部区域没被正确地重建出;CA在较短的运行时间内正确地估计出部分场景的高程信息,但仍有一些区域估计地不准确,如图7(b)中的黑框区域。CANOPUS在CA的基础上,考虑了像素的位置信息,图7(c)中除了较小的黑框区域没被正确地建,整体上比CA的精度更高。ML估计法能在较短的运行时间内估计出精度不错的结果,但在图7(d)中,黑红两块区域产生了较明显的高度突变,对比光学图,这两块区域的估计失真。在ML估计法的结果上加入均值滤波措施后,会将部分孤立的噪声像素过滤掉,但当噪声像素连接成片后,也无法过滤掉这些噪声像素,如图7(e)中的黑框区域。在MAP估计法的结果中,没有出现较明显的高程突变,但建筑物的边缘被模糊化,损失了细节信息。CABMAP的精度较高,对比图7(d)图7(g)中的红框区域,CABMAP有效地减少了ML估计法中的高程突变点,但仍有部分无法消除,如图7(g)中的黑框区域。改进的CABMAP进一步消除了高程突变点,图7(i)是改进的CABMAP的高程3维图,可以看出该算法在实测数据上也表现出较好的精度。综合来看,本文提出的算法在实测数据处理时表现出较高的处理效率和精度,因此具有良好的实用价值。由于数据有限,该实测数据的不模糊高度大于场景高度差,但能部分验证本文算法,后期条件允许,将进一步验证分析本文算法。

      算法 时间(s)
      鲁棒性CRT 0.320831
      CA 0.518862
      CANOPUS 2.009810
      ML估计 0.161702
      ML估计后均值滤波 0.164832
      MAP估计 11.617357
      CABMAP 5.175075
      改进的CABMAP 5.253087

      表 8  算法运行时间

      Table 8.  Run time of the algorithms

    • 本文提出了CABMAP算法,并利用了一种改进措施进行优化。采用了美国Isolation Peak国家公园的仿真数据和TerraSAR-X的巴塞罗那实测数据进行了实验验证。验证结果表明CABMAP算法重建精度高、运行时间短,并且能有效重建被噪声污染区域的高程信息。但同时,该算法只利用了SAR图像的干涉相位和相干系数等信息进行场景的高度重建。未来的工作可以考虑在改进的CABMAP基础上,利用SAR图像的幅度信息提升高度重建精度。

参考文献 (21)

目录

    /

    返回文章
    返回