翰墨点睛：国投坐标
栏头合书：惠民、乔杨
国土名片网钤印：王子廬
插播导读：国土名片网小编芬儿、新宇、张扬、李浩
GB-SAR影像坐标到三维地形坐标转换方法
《长江科学院院报》2018年6期
， (1.苏州科技大学 环境科学与工程学院，苏州 215009；2.武汉大学 a.精密工程与工业测量国家测绘地理信息局重点实验室; b.测绘学院，武汉 430079)
1 研究背景
地基合成孔径雷达(Ground-based Synthetic Aperture Radar，GB-SAR)干涉测量技术是近十多年间发展起来的地面主动微波遥感探测技术。它主要应用3项核心技术：①脉冲压缩技术，如步进频率连续波技术(SFCW)、线调频连续波技术(LFMCW)等，在获取距离域高分辨率的同时又确保了较高的接收机功率和回波信号强度；②合成孔径雷达技术，相比真实孔径的地基雷达，合成孔径能在较短时间内采集雷达辐射场区的二维影像图；③干涉测量技术，通过比较变形前后的观测相位提取变形信息[1-3]。星载SAR受限于卫星重访周期和传感器侧视视角，在地表局部小区域的监测应用中不够灵活。与星载SAR相比，GB-SAR结构简单、安装便捷，它可以为监测目标建立特定的几何场景，并根据目标变形特征灵活确定时间基线。在GB-SAR连续影像序列的数据处理中，由于没有空间基线，无需考虑轨道基线误差、地平相位以及地形相位误差影响。经过十多年的研究发展和应用实践，GB-SAR业已在滑坡、支护边坡、危岩体、水工结构、露天煤矿和冰川移动监测等方面展现了其特点和优势[4-9]，在星载InSAR和传统大地测量监测技术难以满足时空分辨率和精度需求时能够得到较好的应用。
GB-SAR系统依据目标到雷达中心的斜距和偏离雷达天线波束中心线的角度来分辨不同的监测目标。其成像投影方式不同于传统地形图的正射投影以及摄影测量的中心投影，在雷达影像的二维平面坐标系中对变形目标及其位置的识别一般采取人工判读的方式进行。这种方式要求分析人员具备一定的经验，否则容易造成变形目标和区域的错误识别。因此，将雷达二维平面坐标系统和常用的三维地形坐标系联系起来便非常有意义,不仅有利于提高变形体识别的正确性，也有助于长期趋势性变形监测、多种监测手段的融合以及监测成果直接或间接地转化应用。本文深入研究了GB-SAR影像坐标到地方三维坐标变换方法，成果有利于实际GB-SAR监测应用中对变形的三维解译和分析，以更好地理解变形发展的时空规律。
2 相似变换参数估计坐标变换及其误差分析
星载SAR传感器的目标区域满足天线远场近似条件，SAR影像坐标系可视作规则格网。而GB-SAR作用距离相对较近且辐射视场距离跨度大，形成了特殊的扇形格网坐标系。GB-SAR所探测的变形是监测目标实际变形在雷达视线向(Light of Sight，LOS)的投影，各像元以雷达中心点为中心，以到雷达中心的斜距为半径沿圆弧线投影到成像投影面上。GB-SAR的投影方式不仅决定了距离向和方位向分辨单元的划分形式，同时也决定了GB-SAR影像坐标系统变换的方式方法。
为简化分析，进行如下假设：①小区域内高斯-克吕格投影不存在投影变形；②雷达二维平面坐标系和高斯平面坐标系在同一个平面内且轴系转角为0，仅存在平移参数X0和Y0；③高程与高斯平面坐标系构成地面三维坐标系，与XY轴构成右手坐标系。对于GB-SAR内部的几何关系，选定与高斯平面坐标系平行的平面作为雷达成像投影面，则目标点到雷达中心的相对高度就是目标点到该平面的高度，三维地形坐标系与雷达平面坐标系的几何关系如图1所示。
图1 三维地形坐标系与GB-SAR影像平面坐标系的几何关系Fig.1 Geometric relationship between 3D terrain coordinate system and coordinate system of GB-SAR imaging plane
根据上述假设，可以计算得到雷达二维平面坐标系中第i行、第j列的像元Tij在高斯平面坐标系下的坐标为
(1)
式中：θij为相对成像投影面的高度角；(Sij,βj)为目标点在雷达二维平面坐标(xij,yij)的极坐标形式，即
利用相似变换参数估计的模型，由雷达二维平面坐标计算高斯平面坐标，即
(3)
经过一定变形，计算缩放因子m。
(4)
代入式(1)，得
(5)
假设雷达中心线平行于高斯平面，雷达对同样平行于高斯平面的一个目标平面进行观测，雷达中心高出该目标平面距离Hij，给定不同的相对高度Hij，计算缩放因子m。
从图2中可以看出，随着距离的增加，缩放因子m逐渐减小，减小的速率由快到慢，直至趋近于0。由上述分析可知，由相似变换参数估计方法算得的m值大小固定，并不适用于整个雷达辐射区域的坐标变换。该方法仅适用于距离雷达较远、范围较小的局部区域，并且相对于雷达中心的高度变化不能过于剧烈。
图2 缩放因子m的变化规律Fig.2 Variation rule of scaling factor m
图3 坐标直接变换的误差分析Fig.3 Error analysis for direct coordinate transformation
如果不考虑雷达影像投影方式，只进行平面转角旋转和原点平移，将雷达二维平面坐标系转换至高斯平面坐标系中，坐标值产生的点位误差实际上就是雷达沿圆弧线的投影位置与向高斯平面的正射投影之间的距离差值ΔSij，ΔSij与斜距和高度角之间的关系如图3所示。易知，该方法在目标区域分布距离雷达中心线较近时产生的点位误差较小，例如高度角在5°以内，直接坐标变换方法产生的点位误差在4 m以内。但使目标区域均在与雷达中心线高度角[-5°, 5°]之内的要求是十分苛刻的，在实际中往往难以满足。
3 外部高程信息辅助GB-SAR影像坐标到地方三维坐标转换
由上述分析可知，GB-SAR成像特殊的投影方式导致影像坐标相对实际平面坐标发生偏移，且不同高度和位置的地形点成像后偏移的大小是不一样的。实现GB-SAR影像坐标到地方三维地形坐标高精度的转换计算，需要利用外部高程数据(如测区数字地面模型(DTM)/数字高程模型(DEM)、三维地形点云数据等)进行。在三维地形坐标系下的高程数据可依据GB-SAR投影原理计算为雷达投影面的平面坐标，进而利用其与GB-SAR影像像元坐标的匹配关系实现逆向求解像元三维坐标，主要包括以下几个步骤。
3.1 轴系水平旋转
在确定雷达中心的三维地形坐标后，坐标轴系的朝向不一致，因此首先需要进行GB-SAR二维影像坐标系到三维地形坐标系的水平旋转角度的计算。图4显示了GB-SAR影像平面坐标系和地方平面坐标系之间的相互关系。
图4 坐标系轴系关系与水平角度旋转Fig.4 Horizontal angle rotation of coordinate axes
图4中x-y平面是GB-SAR的雷达二维平面坐标系；E-N平面是地方高斯平面坐标系；θ是雷达二维平面坐标系下目标点到雷达中心线方向的偏角；αx-y是x-y坐标系下目标点的极角，αE-N是E-N坐标系下目标点的极角，而βE-N即为我们需要计算的轴系水平转角。通常情况下，轴系水平转角的计算利用公共点进行[10]。但实际监测工作中难有足够小的自然地表可以作为控制点。一般需要布设角反射器并测量出该角反射器的三维空间坐标，同时在雷达影像中识别出该角反射器相应像元的雷达平面坐标。
3.2 三维地形坐标的极坐标化
假设雷达影像中关注区域内一像元点P，其雷达平面坐标为(xP,yP)，相应地表三维坐标系坐标为(EP,NP,ZP)。P点的雷达平面坐标(xP,yP)也可以表示为极坐标形式，如式(6)。
(6)
像元三维空间关系的恢复实际上就是求取雷达影像中的P点在三维地形点中的对应点。为了达到目的，通常先将三维地形数据按照式(7)转化为二维平面上的极坐标形式，从而可以将三维地形数据的极坐标数据与雷达像元的极坐标直接进行比较，这一过程称为三维坐标的极坐标化。
(7)
式中：rTerrain，θTerrain分别为地形点极坐标的极径和极角；(E,N,Z)为地形点三维坐标；(Eradar,Nradar,Zradar)为雷达中心在三维地形坐标系下的坐标。
3.3 GB-SAR成像投影面上像元与地形点的匹配
对三维地形数据极坐标转化后，在其极坐标数据中按照最小距离准则寻找P(rP,θP)的对应点，并将该点的三维坐标作为像元点P(xP,yP)对应的三维空间坐标系下的坐标，如图5。为保证结果的可靠性，在提取像元的三维信息时，亦可在目标点云中选择距离最小的前几个点，计算其三维坐标的平均值作为像元的三维坐标。
图5 极坐标化点云与雷达像元的匹配Fig.5 Coordinate polarization of point cloud and matching with GB-SAR pixel
4 试验案例分析
隔河岩水利枢纽位于湖北省长阳县境内长江支流清江干流上，上距恩施市207 km，下距高坝洲水电站50 km。工程开发任务以发电为主，兼有防洪、航运等综合效益。枢纽工程由混凝土重力拱坝、泄水建筑物、右岸引水式水电站和左岸垂直升船机组成。水电站位于右岸，引水式地面厂房，4条直径为9.5 m的隧洞接直径8 m的压力钢管，分别接至4台300 MW的混流式水轮机组。隧洞采用预应力混凝土衬砌。厂房和压力钢管开挖形成170 m的高边坡，采用混凝土局部置换、设置预应力锚束及加强山体排水等措施处理。
本试验主要对坝体与坝区右岸发电厂房南侧高边坡进行监测。雷达设备安装于大坝下游左岸、距坝体约1 km处，导轨安装基本处于水平状态，雷达视线中心朝向坝体右侧，坝体与右岸高边坡均处于雷达最优视场范围内，如图6。在进行连续变形监测试验之前，首先利用高精度全站仪获取了导轨空间坐标。在导轨两端固定刻画位置处安置固定高度的棱镜基座，分别测量出两端棱镜中心的三维坐标。经计算雷达中心三维坐标为(2 200.608 1，1 109.037 4，106.044 6)(单位：m)。
图6 坝区平面坐标系与雷达平面坐标系Fig.6 Plane coordinate system of dam area and GB-SAR plane coordinate system
表1 测站坐标系与地方坐标系控制点对应关系Table 1 Coordinates of control points in scanner coordinate system and dam area coordinate system
4.1 三维激光扫描仪边坡高程提取
试验中利用地面激光三维扫描仪RIEGL VZ400采集了高边坡及其邻近地物的点云数据，如图7。
图7 三维激光扫描仪地形点云数据采集Fig.7 Acquisition of terrain point cloud data using 3D laser scanner
直接获取的三维点云数据仍然处于扫描仪三维坐标系下，为将地面三维激光扫描仪的点云坐标转换到地方三维坐标系，在边坡上的几个控制点上布设反射片作为控制点，如图7(b)，本次试验中共获取了4个较为理想的控制点，其扫描仪三维坐标和坝区三维坐标列于表1。
点云数据预处理时，将边坡表面的植被剔除，仅保留边坡坡面三维坐标数据。在得到控制点2套坐标之后，通过应用大转角的七参数坐标变换方法，得到边坡点云在地方坝区三维坐标系下的坐标，如图8所示。
图8 坝区坐标系下边坡点云的三维坐标Fig.8 Three-dimensional coordinates of point cloud of a slope in dam area
4.2 轴系水平转角的计算与三维地形坐标的极坐标化
在完成空间坐标系的变换之后，三维地形坐标归算至以雷达中心为原点的水平面上，如图9(a)。由于坝区坐标系的横轴和纵轴定义刚好与雷达平面坐标系的横轴和纵轴反向，所以需要对平移计算后的x,y坐标值做反号处理，便得到局部坐标系下的地形点云，如图9(b)。
图9 坝区坐标系、局部坐标系和雷达坐标系下的地形点云坐标Fig.9 Coordinates of terrain point cloud in dam area coordinate system, local coordinate system,and GB-SAR coordinate system
坝区坐标系和局部坐标系中，电厂厂房结构清晰可辨。此时局部坐标系与雷达平面坐标系的纵轴之间只相差一个数值为锐角的水平转角。三维坐标的极坐标化后的位置只与目标点到原点的斜距和相对于过原点的水平面的高度有关，而与y轴的朝向没有关系。点云数据中特征线由大量三维坐标形成，因此可以较为精确地确定直线极坐标化后的方位角，而建筑结构在雷达影像中所成的像同样占据多个分辨单元，也能够较为精确地确定其方位角。因此本文利用场区特征直线计算轴系水平转角，为0.177 142 rad。利用坐标旋转公式计算极坐标化后的坐标顺时针旋转0.177 142 rad后的坐标，如图9(c)。
4.3 GB-SAR影像像元高程提取
在将三维坐标极坐标化并完成轴系水平角旋转之后，便可以从得到的雷达平面坐标系中寻找GB-SAR影像各像元的对应点。如前文所述，试验中采用到待求像元距离最小作为对应三维地形点目标选择的准则，直接在极坐标化后的坐标系中寻找与雷达影像像元坐标距离最近的目标点，将该点相应的极坐标化前的三维坐标作为像元的新坐标。至此便完成了GB-SAR影像像元的高程提取工作，只要导轨未发生变动，监测获取的所有影像(例如信号强度图、相干图和变形分布图等)都可以按照上述工作求得的关系映射到局部三维表面模型上。
选取高边坡区域的雷达影像，热信噪比(Thermal Signal-to-Noise Ratio，TSNR)阈值设为10 dB，剔除影像中部分信号较弱的区域。对经阈值处理后的TSNR影像坐标值计算相应的坝区三维坐标XYH并绘制到高边坡的三维表面模型上。为去除不在模型区域的像元，在计算前设定淘汰距离为3 m，超过该距离的像元不会出现在三维模型上，如图10。由于试验中高边坡整体偏离雷达中心线的角度有些大，只有高边坡右侧部分区域的雷达反射信号较为理想。第三级马道边缘较为明显，并直接连接4个引水隧道外壁，具有相对较高的TSNR值。该区域GB-SAR影像像元与实际地表结构符合较好。
图10 GB-SAR TSNR图到三维表面模型的映射Fig.10 Mapping of GB-SAR TSNR graph to 3D surface model
5 结 论
(1)由于GB-SAR成像过程投影方式的特殊性，各像元在影像坐标系中的位置相比于实际地形坐标发生径向偏移，使得仅靠人工识图方式对变形区域进行识别与划分容易产生偏差或错误。为更加直接和有效地解译变形监测成果，有必要将雷达影像坐标系统与常用坐标系统统一起来。
(2)利用相似变换参数估计方法进行GB-SAR影像坐标变换模型不严密，在雷达监测视场有一定跨度时易产生较大误差。精确的GB-SAR影像坐标到地方三维坐标变换必须顾及各像元相对于雷达中心的高度信息。本文研究的顾及外部高程信息的GB-SAR影像坐标三维变换方法在理论上能够降低坐标转换误差，但目前缺少对转换精度的直接定量分析。后期笔者考虑引入角反射器阵列以定量分析和评价GB-SAR影像坐标系统的转换精度。
参考文献：
[1] PIERACCINI M, FRATINI M, PARRINI F,etal. High-speed CW Step-frequency Coherent Radar for Dynamic Monitoring of Civil Engineering Structures[J]. Electronics Letters，2004,40(14):907-908.
[2] 王 鹏,周 校.地基SAR干涉测量原理及其形变监测应用研究[J].测绘信息与工程,2012,37(4):22-25.
[3] RÖDELSPERGER S. Real-time Processing of Ground Based Synthetic Aperture Radar (GB-SAR) Measurements [D]. Darmstadt, Germany: Technische Universität Darmstadt, 2011.
[4] MARGOTTINI C, ANTIDZE N, COROMINAS J,etal. Landslide Hazard, Monitoring and Conservation Strategy for the Safeguard of Vardzia Byzantine Monastery Complex, Georgia[J]. Landslides, 2015, 12(1): 193-204.
[5] ATZENI C, BARLA M, PIERACCINI M,etal. Early Warning Monitoring of Natural and Engineered Slopes with Ground-based Synthetic-aperture Radar[J]. Rock Mechanics and Rock Engineering, 2015, 48(1): 235-246.
[6] 邱志伟,岳建平,汪学琴.地基雷达系统IBIS-L在大坝变形监测中的应用[J].长江科学院院报,2014,31(10):104-107.
[7] TARCHI D, ANTONELLO G, CASAGLI N,etal. On the Use of Ground-based SAR Interferometry for Slope Failure Early Warning: the Cortenova Rock Slide (Italy) [M]∥Landslides. Berlin：Springer Berlin Heidelberg，2005：337-342.
[8] TRAGLIA F D, NOLESINI T, INTRIERI E,etal. Review of Ten Years of Volcano Deformations Recorded by the Ground-based InSAR Monitoring System at Stromboli Volcano: A Tool to Mitigate Volcano Flank Dynamics and Intense Volcanic Activity[J]. Earth-Science Reviews, 2014,139:317-335.
[9] SEVERIN J, EBERHARDT E, LEONI L,etal. Development and Application of a Pseudo-3D Pit Slope Displacement Map Derived from Ground-based Radar[J]. Engineering Geology, 2014, 181: 202-211.