基于高分三号SAR数据的城市建筑高分辨率高维成像

毕 辉 金 双 王 潇 李 勇 韩 冰 洪 文

①(南京航空航天大学电子信息工程学院 南京 211106)

②(南京工业大学计算机科学与技术学院 南京 211816)

③(中国科学院空天信息创新研究院 北京 100094)

1 引言

与传统光学观测手段不同，合成孔径雷达(Synthetic Aperture Radar,SAR)具有全天时、全天候的工作能力，因而在国土资源勘测、自然灾害监测等领域得到广泛应用[1]。然而传统SAR成像只能获取目标的方位-距离二维图像，无法准确反映目标的三维散射特征，一定程度上影响了SAR图像的进一步应用。层析合成孔径雷达(Synthetic Aperture Radar Tomography,TomoSAR)是SAR成像技术进一步扩展。它将合成孔径原理延伸至高程向，可基于多景二维复图像获得目标的方位-距离-高程信息，进而实现三维成像[2]。1999年，Reigber等人[3,4]首次展示了机载TomoSAR成像结果，实现了基于谱估计方法的森林区域三维成像。2006年，Fornaro和Serafino[5]基于长时间基线星载SAR数据，开展了TomoSAR成像星载实验，证实了谱估计技术应用于高程向重构的可行性与有效性。2011年，Reale等人[6]证实了高分辨率数据和先进干涉处理技术相结合可以更好地重建建筑物三维结构。2013年，Shahzad和Zhu[7]提出了一种全新的建筑物立面重建方法，获取了拉斯维加斯百乐宫酒店的三维点云，证实了TomoSAR点云在构建动态城市模型方面的优越性。2017年，Wang等人[8]提出了迭代重加权的交替方向乘子算法，用于实现快速TomoSAR成像。2018年，Wang和Zhu[9]提出了一种基于核主成分分析的TomoSAR成像方法，使用极小代价即可分离同一方位-距离分辨单元中沿高程向分布的多个散射体。2019年，秦斐等人[10]针对TomoSAR成像中高程向分辨率较低、建筑物叠掩、提取建筑物目标特征效率较低等问题，提出了一种基于机器学习的建筑物目标识别和提取算法，提高了观测目标的特征提取效率，并通过机载阵列三维SAR实验数据验证了方法的有效性。差分层析合成孔径雷达(Differential Synthetic Aperture Radar Tomography,D-TomoSAR)是TomoSAR的进一步扩展。它基于多基线观测数据可实现对目标的方位-距离-高程-时间四维成像，不仅解决了TomoSAR成像中的高度错位与模糊问题，还可高精度获取目标的形变信息。D-TomoSAR概念由Lombardini于2003年首次提出[11]。2007年，Fornaro等人[12]证明了D-TomoSAR可作为传统永久散射体监测技术的有效替代方案，实现对大场景形变的有效监测。2008年，Fornaro等人[13]提出了一种可分离干扰散射体相关时间序列的技术，并应用于罗马地区的D-TomoSAR成像中，证明了该技术可用于城市复杂场景的形变监测。2009年，Zhu等人[14]基于多景TerraSAR-X数据，实现了对拉斯维加斯的D-TomoSAR成像，获取了城市四维雷达地图。2010年，Fornaro等人[15]使用ERS数据对罗马市Grotta Perfetta地区进行了D-TomoSAR成像，结果显示差分层析成像方法可以有效区分同一分辨单元中的不同散射体，克服了差分干涉的局限性，进一步提高了对城市基础设施形变的监测能力。2011年，Zhu和Bamler[16]提出了一种“时间扭曲”方法，并基于TerraSAR-X数据实现了城市区域的D-TomoSAR成像，有效获取了城市建筑的线性运动和季节性运动速率。2015年，Siddique等人[17]将D-TomoSAR方法和永久散射体干涉技术相结合，基于50幅TerraSAR-X图像获得了巴塞罗那中高层建筑物立面散射体的时空反演结果。2020年，Wang和Liu[18]提出了一种广义D-TomoSAR成像系统模型和一种基于拟极大似然的成像算法，同时反演出了建筑物线性运动、季节性运动等多个形变运动速率。目前，TomoSAR和D-TomoSAR成像技术已在城市建筑物和基础设施的三维重建和长期形变监测方面展现了极大的应用潜力。

压缩感知(Compressive Sensing,CS)是一种重要的稀疏信号处理技术，它可以使用较少样本实现对稀疏信号的高质量恢复[19–21]。CS-TomoSAR成像的前提是观测场景的高程向分布稀疏，而城市区域主要为人造建筑，其高程向分布都满足稀疏性条件。因此，CS在城市区域三维、四维成像中具有广阔的应用前景。起初，由于现代米级分辨率星载SAR系统轨道限制，TomoSAR成像高程向分辨率远低于方位向和距离向，因而迫切需要超分辨算法来解决这一问题，2010年，Zhu和Bamler[22]介绍了一种基于CS理论的TomoSAR成像方法，相比于传统谱估计方法，该算法实现了对目标高程向分布的超分辨重构。2012年，Zhu和Bamler[23]将所提出的SL1MMER算法应用到TerraSAR-X星载数据处理中，获得了拉斯维加斯百乐宫酒店的高分辨率TomoSAR成像结果，并证明了该算法具有超分辨能力。2015年，Weiss等人[24]提出了一种适用于TomoSAR的自适应CS算法，准确识别了同一分辨单元中的两个散射体位置。2017年，Li等人[25]研究了基于SPICE的TomoSAR成像方法，并利用8幅TerraSAR-X条带影像，实现了对内蒙古根河市某建筑的高精度三维重建。2010年，Zhu和Bamler[26]将CS技术应用在D-TomoSAR成像中，证明了在高程向多散射体分离上CS技术相比于传统谱估计方法的优越性，同时指出CS可自动识别散射体数量，非常适用于星载SAR系统的三维、四维成像。2010年，Zhu和Bamler[27]基于SL1MMER算法重构获得了拉斯维加斯会议中心建筑群的D-TomoSAR成像结果，展示了CS技术在四维成像方面的优势和能力。2014年，Leng等人[28]将最小绝对收缩和选择算子CS算法应用于建筑区域，并展示了巴塞罗那的高程重建和形变监测结果。

高分三号卫星是我国首颗分辨率达到1 m的C频段多极化SAR卫星，于2016年8月10日在中国太原卫星发射中心由长征四号丙运载火箭发射升空[29]。高分三号卫星是“国家高分辨率对地观测系统重大专项”中唯一的民用微波遥感成像卫星，具有高分辨率、大成像幅宽、多成像模式、长寿命运行等特点，可实现全天时、全天候的全球海洋与陆地监测[30]。目前高分三号已成功应用于高精度测绘、自然灾害监测等多个领域[31,32]。然而，由于设计之初未考虑后续高维成像应用，现有高分三号获取的SAR图像存在有一定的空间、时间去相干问题，对其进一步应用于干涉SAR、差分干涉SAR、TomoSAR、D-TomoSAR等存在一定挑战。2019年，余博等人[33]基于高分三号数据，对河南省登封市周围地区进行了干涉测量实验，通过与哨兵一号卫星数据结果进行对比分析，验证了高分三号的干涉能力以及可提取地表形变信息的能力。2021年7月，黄震等人[34]基于高分三号SAR数据进行了干涉测量实验，成功提取了观测区域的数字高程模型。

本文基于7景高分三号SAR复图像数据，利用CS技术，开展了TomoSAR和D-TomoSAR成像实验研究，获取了北京市雁栖湖周围建筑的高分辨率三维、四维SAR图像，实现了建筑物的高质量三维重建以及高精度形变监测，为后续基于高分三号SAR数据的干涉系列应用及多维高分辨率成像提供了技术支撑。

本文后续结构如下：第2节主要介绍TomoSAR,D-TomoSAR成像模型，并给出了上述两个模型的CS求解方案；第3节介绍了本文所使用的高分三号SAR复图像数据集；第4节基于仿真数据，开展了TomoSAR,D-TomoSAR成像实验，证明了CS技术在高分辨三维成像及高精度形变监测方面的有效性；第5节基于7景高分三号SAR复图像数据，对北京雁栖湖周围建筑进行了TomoSAR,D-TomoSAR成像研究，获取了代表性建筑和大观测区域的三维、四维雷达图像；第6节对文章进行了总结和展望。

2 成像模型

2.1 TomoSAR成像模型

TomoSAR利用对同一场景观测获取的多幅配准的二维SAR复图像(多基线观测数据)在高程向上进行孔径合成，以获得高程向上的分辨能力，进而重构目标的三维散射信息[3,4,35,36]。TomoSAR成像几何如图1所示。设共有N条基线用于数据获取，令bn(n=1,2,...,N)表示高程向孔径分布，对于一个选定的方位-距离分辨单元，第n幅SAR图像对应的聚焦测量值可以表示为

图1 TomoSAR成像几何Fig.1 TomoSAR imaging geometry

式中，ξn=-2bn/(λr)表 示高程向频率，其中λ为波长，r为斜距；γ(s)表 示沿高程向s的复反射函数；Δs为 高程向跨度。沿高程向s对高程向复反射函数γ(s)进行离散化，则式(1)中的成像模型可近似表示为

式中，L为高程向离散化点数；g=[g1,g2,...,gN]T表示测量值向量；R=exp(-j2πξnsl)为根据Tomo-SAR成像几何所构建的观测矩阵；γ=[γ(s1),γ(s2),...,γ(sl)]T表示高程向离散复反射函数，其中sl(l=1,2,...,L)为离散高程向分布。从式(2)可以看出，TomoSAR成像模型可视为对γ(s)不规则采样的离散傅里叶变换。因此，一个SAR测量值可以看作目标复反射函数沿高程向的一个谱参数。对于非参数化谱分析问题，高程向理论分辨率ρs依赖于高程孔径大小 Δb，在高程向采样密集的情况下，ρs可由式(3)进行计算。

2.2 D-TomoSAR成像模型

相比于TomoSAR,D-TomoSAR在三维的基础上多了一个时间方向的维度。它沿高程向和形变速度向合成两个孔径，进而获取被观测目标的高程和形变速率的联合分辨率，实现四维成像[11]。对于N个复图像而言，当高程孔径位置为bn、时间基线为tn时，第n幅图像的聚焦测量值可以表示为

近年来，随着我国社会经济的迅速发展，我国建筑工程项目数量日益增加，而工程测量是保证建筑工程项目正常进行的前提。传统的工程测量必须要花费很多时间、人力和物力，并且不能确保测量精准度。地面三维激光扫描技术具有高效率、高精准度等优点，将其在工程测量中应用，能够明显提升工程测量精准度。

式中，ηn=-2πtn/λ表 示形变速度频率；V(s)表示形变速率。式(4)中的模型也可写为

式中，Δv表 示被观测目标的形变速率跨度；δ(·)是与形变项相关的谱分布。令aγ(s,v)=γ(s)δ(v-V(s))，则式(5)可以写成

式(6)中的模型可以视为a(s,v)在高程-形变平面的二维傅里叶变换。因此，其在立面轴上的投影为反射率剖面γ(s)[37]。将式(6)中的s和v离散化后，D-TomoSAR成像模型可以表示为

其中，g=[g1,g2,...,gN]T表 示测量向量；R=exp(-j2π·(ξnsl+ηnvq))为D-TomoSAR成像观测矩阵；sl(n=1,2,...,L)为 离散高程向分布；vq(q=1,2,...,Q)为 离散形变向分布；γ由离散化的a(s,v)组成。若时间孔径大小 Δt，则形变分辨率ρv可由式(8)进行计算

2.3 CS重建

城市中被观测目标主要为人造建筑，其高程向分布通常都是稀疏的，即每个方位-距离分辨单元中的散射体个数有限。因此，当测量矩阵R满足有限等距性质条件时，面向式(2)和式(7)中的模型，本文通过解决如下的最优化问题分别实现基于CS的TomoSAR和D-TomoSAR成像。

式中，β为正则化参数，与噪声水平和样本数目有关。CS算法可以在短时间内从获取的样本数据中实现高质量信号恢复[38,39]。基于该算法的优势，本文采用CS算法进行TomoSAR和D-TomoSAR成像。

3 高分三号数据集简介

本文所使用的高分三号数据集总共包含7景复图像，具体参数如表1所示。该数据集的7景图像的空间基线孔径大小约为1417 m，时间基线跨度是从2018年6月到2019年9月，共464 d，表2给出其时空基线具体参数，其时空基线分布情况如图2所示。本文以2019年3月1日获取的SAR图像为主影像，其余6景为辅图像。辅图像时空基线位置是相对于主影像计算得到的。本文第4节将基于表1的参数进行点目标仿真实验，第5节将基于该高分三号数据集进行实验。

图2 高分三号数据集时空基线分布图Fig.2 Spatial-temporal baseline distribution of GF-3 dataset

表1 高分三号数据集参数Tab.1 Parameters of GF-3 dataset

表2 高分三号数据集时空基线参数Tab.2 Spatial-temporal baseline parameters of GF-3 dataset

4 仿真实验

基于表1中的实验参数，本文设定高程向分布有两个散射体，模拟产生了7景仿真数据，添加信噪比为20 dB的噪声，并分别使用经典谱估计方法和CS算法对高程向进行TomoSAR和D-TomoSAR成像[40]。本节将展示3种经典谱估计算法即Beamforming (BF)[41],Adaptive beamforming (Capon)[42]和Multiple signal classification (MUSIC)[43,44]的仿真结果，用于与CS算法重构结果进行比较，以说明CS算法在TomoSAR,D-TomoSAR成像中的优势。图3为高程向两个散射点的TomoSAR成像结果，横坐标为高程向位置分布，纵坐标为散射点的幅度值。本文设置两个散射体的距离分别为11 m和50 m。由图3可以看出，当距离为11 m时，由于小于高程向理论分辨率，3种谱估计算法的成像结果均存在分辨率较低和模糊严重等问题，无法准确分离两个散射点，造成重建失败；而CS算法仍可有效识别两个散射体，实现高程向的超分辨成像。当距离为50 m时，可以看出，BF算法重构结果具有严重的、不规则的旁瓣；Capon和MUSIC算法相较于BF算法，重构结果具有更低的旁瓣，提升了TomoSAR高程向重构质量；而相较于3种谱估计算法，CS则可以更加有效地抑制旁瓣和噪声，进一步提升了高程向散射体的可分辨能力。图4为高程向3个散射点的TomoSAR成像结果，3个散射点的间隔是20 m，与两个散射体的成像结果相似，可以看出CS算法对于多散射体分离也具有很好效果。

图3 高程向两个散射点的TomoSAR成像结果(左图：两个散射点之间的距离为11 m；右图：两个散射点之间的距离为50 m)Fig.3 TomoSAR reconstructed reflectivity profiles of two scattering points along the elevation direction (left image:the distance between two scattering points is 11 m;right image:the distance between two scattering points is 50 m)

图4 高程向3个散射点的TomoSAR成像结果(3个散射点之间的间隔为20 m)Fig.4 TomoSAR reconstructed reflectivity profiles of three scattering points along the elevation direction (the distance between three scattering points is 20 m)

图5 D-TomoSAR仿真结果(两个散射体高程位置为–10 m,10 m；散射体形变速率分别为4毫米/年、–7毫米/年)Fig.5 D-TomoSAR simulation results (elevation position of two scatters are–10 m and 10 m;deformation velocity of two scatters are 4 mm/year and–7 mm/year,respectively)

5 高分三号实验结果

5.1 代表性建筑物反演

本节将展示基于此高分三号SAR数据集的Tomo-SAR和D-TomoSAR重建结果。首先选取观测场景内两处代表性建筑进行了TomoSAR和D-TomoSAR成像，分别为某生态农业公司和雁栖湖会展中心。图6(a)所示为北京某生态农业公司的光学图像，红色虚线框出的区域是所关注的目标。图6(b)为该生态农业公司的二维SAR图像。图7(a)为TomoSAR重建结果，获取了区域内5栋建筑的高度信息。由图7(a)可以看出，5个独立建筑的高度均为35 m，这与实际该5栋完全相同的建筑高度相符，验证了CS-TomoSAR成像技术的有效性，说明其可以用于高分三号SAR数据的三维高精度成像中。图7(b)为基于CS算法的D-TomoSAR重建结果，可以看出，5栋建筑均有不同程度的形变，这可能与建筑材料随着季节变化的热胀冷缩紧密相关，在日常监测中应该有所关注，防止出现相关危险。为更加直观地展示所关注场景，图8给出了生态农业公司的三维点云，更加准确地反映出了建筑物的三维散射结构。

图6 生态农业公司Fig.6 Ecological agricultural company

图7 生态农业公司高程图及形变速率图Fig.7 Elevation and deformation velocity maps of ecological agricultural company

图8 生态农业公司三维点云图Fig.8 3-D point cloud of ecological agricultural company

图9(a)为北京雁栖湖国际会展中心的光学影像，该建筑真实高度约30 m。图9(b)为雁栖湖国际会展中心的SAR图像。图10(a)给出了CS-TomoSAR成像结果，并在图11给出了其三维点云图。由图10(a)和图11可以看出，重建结果基本反映了建筑物的三维结构，重构的建筑物高度与实际相吻合，这也说明了本文基于高分三号SAR数据，已可实现对复杂建筑物较高精度的三维重建。图10(b)展示了基于CS技术的雁栖湖国际会展中心D-TomoSAR重建结果。由图10(b)可以看出，该建筑的形变大约在–10毫米/年～10毫米/年之间，且建筑的左半部分和右半部的线性形变速率正好相反，这反映了该建筑下方地面可能处于一边抬高一边塌陷的变化之中，应当着重关注。

图9 北京雁栖湖国际会展中心Fig.9 Beijing Yanqi lake international convention and exhibition center

图10 北京雁栖湖国际会展中心高程图及形变速率图Fig.10 Elevation and deformation velocity maps of Beijing Yanqi lake international convention and exhibition center

图11 北京雁栖湖国际会展中心三维点云图Fig.11 3-D point cloud of Beijing Yanqi lake international convention and exhibition center

5.2 大场景三维、四维成像

本节展示了基于高分三号数据集的大场景TomoSAR和D-TomoSAR成像结果。图12分别给出了所选取的顶秀美泉小镇区域的光学影像和SAR图像。图13(a)给出了基于CS算法的TomoSAR成像结果，重构图像显示该区域内建筑高度均在15 m至20 m之间，准确反映了小镇内均为4～6层居民楼建筑的实际情况。图13(b)给出了基于CS算法的D-TomoSAR成像结果。其显示小镇左下角区域和右上角区域的建筑形变约为–10毫米/年，可能由于这两块区域周围的某些施工导致地面沉降，小镇中心的形变速率约为5毫米/年，总体来说比较稳定。该结果表明，基于高分三号SAR数据，已经可以实现大场景的高质量三维重建以及高精度形变监测，进一步验证了高分三号SAR卫星在城市感知与监测中的应用潜力。

图13 顶秀美泉小镇区域的高程图及形变速率图Fig.13 Elevation and deformation velocity maps of Dingxiumeiquan town

6 结论

本文开展了基于高分三号SAR复图像数据的TomoSAR和D-TomoSAR成像研究，获取了北京雁栖湖地区两处代表性建筑的三维、四维成像结果，并给出了大面积观测场景的三维、四维雷达图像。研究结果显示了我国高分三号SAR卫星应用于TomoSAR和D-TomoSAR成像方面的潜力，为后续拓展高分三号的干涉系列应用领域提供了技术支撑。

后续本团队将继续开展相关研究工作，收集高分三号多个区域、多景数据信息，同时探索全新的高分辨TomoSAR,D-TomoSAR成像算法，实现面向复杂城市场景的大范围三维、四维成像，进一步挖掘高分三号在干涉系列应用中的巨大潜力。