时序合成孔径雷达干涉测量的多级化解算

胡凤鸣, 伍吉仓

(同济大学 测绘与地理信息学院,上海 200092)

上海磁悬浮全长30 km,是世界第一条商用的高架磁悬浮运营线.磁悬浮工程在建设和运营期间对沿线的地面沉降控制要求严格,每一跨径的沉降差异须小于3 mm[1]. 上海浦东地区第四纪松散土层厚300余米,极易受地下水开采以及工程活动的影响诱发地面沉降[2].为保证磁悬浮的安全运营,需要进行高频高密度高精度的沉降监测.传统的精密水准测量由于空间密度较低,测量时间间隔大,难以满足磁悬浮沿线沉降监测需求.

合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)是一种新兴的空间对地观测技术,能够全天候、大范围、高精度获取地球表面地形和形变信息.针对传统InSAR技术时空间失相干以及大气延迟影响[3]提出的多时序InSAR(MT-InSAR)可以有效提高测量精度.研究表明,通过分析处理高质量的相干点可以获得亚毫米级的形变监测精度[4].该技术已经被成功应用于地壳形变[5]、城市地表[6-7]以及大型结构体的监测中[8-9].X波段的TerraSAR-X/TanDEM-X以及COSMO-Skyed系统获取的高分辨率SAR数据使得该技术可以获取更加精细的地表信息[10-11],这对于城市区域内大型结构体健康监测具有重要意义.

MT-InSAR技术仅使用在时间序列上始终保持高相干性的散射体,依据干涉相位的时空特性分离不同的分量,从而获取高精度的形变速率.因此选取点的质量直接影响最终的解算结果,选取高相干点的方法可以分为基于点和弧段两类,前者包括归一化振幅离差[12-13],记为NAD,以DA表示,信噪比[14],相干性[15],相位分析[16],极大似然法[17],相干矩阵特征值分析[18]等,这些方法通过依据不同原则分析散射体在时间序列上相位的稳定性来选取高相干点.后者首先进行点分析得到较多的备选点集合[19-20],之后分析弧段质量以及点之间的连通性得到高相干点.选取高质量的相干点并保证空间点密度之间保持平衡是MT-InSAR技术的一大难点.宽松的阈值可以得到高密度的观测结果,但是测量精度随之降低,反之可以保证精度,但是观测结果却过于稀疏.因此需要研究针对不同质量的相干点需采用不同的解算方式.目前的解决方案分为两步,利用较小阈值选取较高质量相干点,解算结果作为一级参考网,保证解算精度;之后扩大阈值,构成二级网求解一般质量相干点,以保证空间点密度[21-22].但是现有的空间点分级方法没有考虑空间点密度的不均匀性,不同级相干点连接方式以及精度评价方法也有待完善,并且计算效率不高,使用原始分辨率数据的大范围处理通常难以实施.

本文首先阐述了基于弧段选取的多级化MT-InSAR技术处理方法,包括基本观测量,参考网和附加网的求解和精度评价方法,并对每一级网络中相干点的选取以及求解策略作了详细讨论.利用C波段的ENVISAT、Sentinel-1A以及X波段的TerraSAR-X数据对上海磁悬浮及其沿线区域开展沉降监测,分析不同数据的适用性,并与精密水准数据验证,为城市大型结构体健康监测提供参考.

1 多级化MT-InSAR技术

1.1 基本观测量

在MT-InSAR技术中,基本观测量为双差相位观测值.假设在时间序列上有m幅已配准的合成孔径雷达(synthetic aperture radar,SAR)数据,选取其中一幅作为主影像后可以得到m-1幅干涉相位图.干涉图中任意两个像元i与j构成弧段,去除平地相位后,第k幅干涉图中弧段的相位差可以表示为[23]

(1)

式中:Δhi,j和Δvi,j分别为点i与j之间的高程改正之差和形变速率差;ni,j,k是相位模糊度;Bcik和tk为垂直基线和时间基线;Ri为像元点i与卫星间的斜距;λ为雷达波长;θi为入射角;Ci,j为相位常数,与主影像的大气延迟有关;ei,j,k为相位噪声.对式(1)相位解缠后使用最小二乘估计即可得到参数解和对应的参数精度[23].

在基于弧段选取的MT-InSAR技术中,初始相干点的选取包含大量的噪声点,最终可靠相干点的筛选依赖于选取高质量弧段解.弧段质量一般采用时序相干性来衡量[9].由于模型中采用了线性形变速率的假设,因此加入线性速率的假设检验进一步筛选弧段.此二者构成文献中的弧段可接受条件[20].

基于本节得到的弧段解进行多级化相干点构网求解.由于城市区域的相干点之间相关性较低,因此采用归一化振幅离差方法进行相干点选取,该方法在选取相干点时仅考虑相干点的时间序列噪声,并通过选取不同的阈值来控制选取相干点的质量,选点的阈值越低,获取点的精度越高,但是点的数目越少.数据处理的主要流程见图1.

图1 多级化相干点处理流程

Fig.1 Flowchart of processing coherent points with multi-level strategy

1.2 参考网

参考网的求解精度直接决定了最终结果的精度,构成该网络的相干点应当保证有足够高的精度.因此采用归一化振幅离差法选取相干点时使用尽可能低的阈值.考虑到大气效应的影响,求解弧段不宜过长[13],参考网中相干点在空间内应尽可能均匀分布.依据这两个准则,参考网相干点选取的阈值使用0.25,是经典永久散射体干涉测量方法中高相干点的选取阈值[12].这一步求解的目标是将尽可能多的点连接起来,采用二次构网策略.首先使用Delaunay三角网连接所有点,筛选弧段后,对孤立的点集采用渐进式星形网构成新的弧段.详细处理步骤可参见文献[20].

(2)

式中:B1为构成参考网弧段相关的设计矩阵;Δv1为求解的所有弧段形变速率之差;(·)+是求解广义逆运算符;P1为弧段解的权阵,可采用如下定义:

(3)

参考网的形变速率估计精度计算方式为

(4)

式中:σ1为参考网的单位权中误差.

当测区范围较大或是处理高分辨率数据时,法方程系数矩阵阶数较大,广义逆求解缓慢或是由于矩阵过大而难以求解,因此需要减少构成参考网的相干点数目.若使用降低阈值减少相干点的方式会导致相干点在空间分布不均匀,造成最终结果的精度在空间不均匀分布.此外,附加网在缺少参考网的相干点的区域无法生成弧段,从而造成相干点的遗漏.因此,本文采用基于点密度的分层抽样方法来对空间点进行稀疏,在保持相干点均匀分布的同时,减少空间点密度较大区域的相干点数目.抽样后得到的相干点组成参考网,余下的点使用附加网进行求解.

1.3 附加网

采用宽松阈值下的归一化振幅离差方法选取相干点,使用星形网将这些点与参考网中的点连接构成新的弧段.构网时仍然采用渐进式策略,每次将待求点与若干个参考网中的点连接,循环处理直到该点与参考网之间有弧段连接,否则剔除该点.联合参考网中相干点的参数对新求解的弧段进行积分可以得到附加网中相干点的形变速率:

(5)

式中:B2为构成附加网弧段相关的设计矩阵;Δv2为求解的所有弧段相对形变速率之差;P2为权阵,与P1为相同定义方式;k为拉格朗日函数的系数矩阵;v2是附加网点待求解的形变速率;G为限制矩阵,定义如下:

Gv2=v1

(6)

附加网的形变速率的精度计算方式为

(7)

式中:σ2为附加网的单位权中误差,Ql定义如下:

(8)

其中,B11和B22是B2的子矩阵,定义如下:

B2=[B11|B22]

(9)

求解参考网后,附加网可以使用任意多次以获取更多相干点的信息.求解附加网时,不同相干点的处理是相互独立的.当区域较大时,可以将整个区域分为独立的若干块,对每块区域内的相干点依次处理,从而减少运算的内存开销,适用于大范围的数据处理.

2 数据处理结果

本次实验采用2007~2018年间的31景C波段ENVISAT、53景X波段TerraSAR-X和51景C波段Sentinel-1A影像,详细数据参数见表1.不同卫星平台数据均完全覆盖磁悬浮沿线区域.由于上海地势较为平坦,且楼房众多,处理中不再使用地面高程数据进行初始高程相位改正,仅改正平地相位.地面点高程引起的干涉相位在求解中通过时间序列方法进行分离.

表1 SAR数据主要参数

以高分辨率TerraSAR-X数据为例,验证多级化数据方法的有效性.构成参考网的相干点采用经典的永久散射体[12]判断方法进行选取,DA阈值取为0.25,该阈值可以保证选出的相干点在整个空间均有分布.选取初始相干点数目为472 929个,点分布参见图2a,从图中可以发现相干点在磁悬浮沿线呈现明显的线性分布规律,磁悬浮周边房屋区域也存在着高密度的相干点.初始相干点数目多,难以直接求解,需要降低参考网相干点数目.

根据1.2节中的处理策略,首先使用降低阈值法,将初始选点的阈值更改为0.1,可以得到17 680个相干点,点分布可参见图2b.从该图可以发现,相干点在空间的分布很不均匀,磁悬浮沿线几乎不存在相干点.在多级化数据处理中,估算精度随着分级的进行而逐步递减,若磁悬浮沿线的相干点只存在于附加网中,那么这些点的精度相对于整个测区是偏低的.因此降低阈值选取参考网的相干点是不合理的.

使用基于空间点密度的分层抽样方法来合理选择参考网的规模.基于图2a的初始相干点集,先计算相干点的空间密度.取100 m为搜索半径,对每个点统计搜索窗口内的相干点数目,从而得到各点的空间相干点密度.为简化运算,将空间点密度量化为等间距的5级密度,点密度分布见图3.由图3可以发现,磁悬浮沿线相干点密度远低于周边房屋区域.获取空间点密度后进行分层抽样,依据空间点密度将空间点分为5层,每层抽样时采用等间距随机选点,空间密度大的区域采用较大的采样间隔,反之则采用较小的采样间隔.在本次试验中,5层抽样的采样间隔比为1∶2∶3∶4∶5.设置目标抽样点数目为15 000,可以得到5层的采样间隔依次是14,28,42,56,70.进而从472 929个初始相干点中抽样得到17 187个相干点,这与降低阈值法得到的相干点数目相当,分层抽样得到的相干点分布见图2c.对比图2b可以发现,相干点在空间的分布更加均匀,并且较好地保留了磁悬浮沿线的相干点.

a 初始选取相干点

b 降低选点阈值法

c 基于点密度的分层抽样法

图2 参考网相干点分布

Fig.2 Distribution of coherent points in the reference network

图3 空间点五级密度分布

求得参考网中相干点的参数后,扩大DA阈值至0.4,选取中等质量相干点,并利用附加网求解分层抽样余下的相干点以及中等质量的相干点,为保证高质量相干点的精度,将附加网分为两次进行,分层抽样余下的相干点构成一级附加网,中等质量相干点构成二级附加网,各级网获取的相干点数目见表2.

表2 各级网络获取的相干点数目

Tab.2 The number of accepted coherent points in each network

参考网一级附加网二级附加网初选点数17 155470 1392 210 867接受点数15 813435 2341 711 388通过率/%92.292.677.4

基于弧段解的时间序列分析方法与选点策略无关,在数据处理的流程中通过网络连通性分析来去除噪声点,以保证最终结果的可靠性.初始选点质量越高,经网络连通性分析后可接受的相干点越多.由表2可以发现,参考网与一级附加网中的相干点是通过同一DA阈值选取的,因此有着相似的通过率,二级附加网的相干点由于选点阈值较为宽松,噪声点多,因此可接受相干点的通过率低于前者.

Sentinel-1A与ENVISAT数据的处理使用类似的多级化处理方式.不同数据集得到的磁悬浮沿线及其周边区域的视线向形变速率如图4所示.为了尽可能突出沉降结果的细节信息,所有结果均保持SAR图像的原始分辨率,未做多视处理.

不同卫星数据在磁悬浮轨道及其周边区域得到的相干点数目依次为:ENVISAT 129 182个,Sentinel-1A 356 663个,TerraSAR-X 2 165 435个.由图4b和4c可知,Sentinel-1A与TerraSAR-X得到的结果中,磁悬浮沿线均存在大量呈线性排列的相干点,而在ENVISAT的结果中(图4a),磁悬浮沿线相干点极少.Sentinel-1A与ENVISAT均为中等分辨率C波段数据,且Sentinel-1A数据的方位向分辨率低于ENVISAT数据.由此可以得知磁悬浮沿线上相干点的数目与距离向分辨率有关,而与方位向分辨率关系较弱.ENVISAT数据由于距离向分辨率较低,在磁悬浮沿线难以形成有效的点目标,不能获取磁悬浮轨道上足够的形变信息.

对比表1中不同传感器的测绘幅宽,尽管TerraSAR-X数据分辨率较高,但是单景数据覆盖范围仅为Sentinel-1A数据的1/30.ENVISAT数据与Sentinel-1A数据有着相似的空间分辨率,但是Sentinel-1A数据的覆盖范围是ENVISAT数据的5倍,且相干点的密度更高.因此Sentinel-1A数据的可用性更好.

此外,选取一块形变显著的区域(图4中白色方框区域),详细研究磁悬浮周边区域沉降对轨道的影响.不同卫星数据获取的该区域沉降速率见图5.图5中表示沉降速率的色条与图4一致,白线表示磁悬浮轨道位置.从图5中可以发现,在2007~2011年间,磁悬浮周边区域沉降较大,最大年沉降速率大于20 mm·年-1(图5a);2011~2015年间,该区域趋于稳定,但是仍然存在沉降速率大于15 mm·年-1的点(图5b);2016~2018年间,该区域大多数点保持稳定,沉降速率小于10 mm·年-1(图5c).不同数据得到的结果中沿轨道的线相干点密度有很大差异,

a ENVISAT

b TerraSAR-X

c Sentinel-1A

Fig.4 Subsidence velocity along the line of sight within the maglev zone

但是这些点在不同时段均保持稳定,年平均沉降速率小于3 mm·年-1.因此周边区域的沉降并未影响磁悬浮轨道的安全.

a ENVISAT

b TerraSAR-X

c Sentinel-1A

图5 沉降较大区域放大图

Fig.5 Magnified view of the selected area with large subsidence

3 结果验证

利用实测水准数据验证本文多级MT-InSAR结果的可靠性.精密水准包含2009~2014年间5期数据,与TerraSAR-X数据获取的时间段一致.该区域共26个水准点(BM1~BM26),位于磁悬浮轨道以及周边区域,具体位置见图6.

图6 磁悬浮区域水准点分布

Fig.6 Distribution of levelling points within the maglev zone

水准测量基于地面单个点,而SAR图像中的单个像元反映了地面某一块区域的综合信息,二者不能直接对比.为了获取MT-InSAR结果中与水准点对应的沉降信息,以每个水准点为中心,选取30 m×30 m的搜索窗口,以搜索窗口内所有相干点的加权平均值作为与水准点对应的结果.MT-InSAR结果表明磁悬浮轨道与周边区域呈现不同的沉降信息,在使用搜索窗口时必须将轨道与周边区域相干点分开处理.因此水准数据的验证分为两组,磁悬浮轨道上点以及轨道周边区域的点.

使用缓冲区法分离磁悬浮轨道上点,选取10 m缓冲区半径,在磁悬浮沿线建立缓冲区,位于缓冲区内部的点视为磁悬浮轨道相干点.图7表示磁悬浮轨道上水准点与MT-InSAR相干点的形变速率对比结果,二者结果差异的平均值是-0.020 4 mm·年-1,标准差是0.628 3 mm·年-1.图8表示磁悬浮周边区域内水准点与MT-InSAR相干点的形变速率对比结果,二者差异的均值为-0.414 4 mm·年-1,标准差是1.170 8 mm·年-1.MT-InSAR的结果中,磁悬浮轨道上点与水准测量结果符合度较高.而磁悬浮周边区域中相干点难以与水准点直接对应,因此符合度较低.

图7 磁悬浮轨道上水准点与MT-InSAR形变速率对比

Fig.7 Comparison between deformation velocities derived from levelling and from MT-InSAR along the maglev track

图8 磁悬浮周边区域水准点与MT-InSAR形变速率对比

Fig.8 Comparison of deformation velocity between levelling and MT-InSAR nearby the maglev track

选取4个水准点BM6、7、11、12,对比水准与MT-InSAR获取的时间序列累积形变量,如图9所示.由图可以发现,MT-InSAR得到的形变时间序列与水准结果有着良好的一致性.并且水准结果表明BM11与BM6在2013年10月至2014年期间出现了明显的抬升现象,这与MT-InSAR获取的结果有着相似的形变趋势.2013年开始出现的抬升现象与上海市实行的限制地下水开采以及地下水回灌有关[24].

a 水准点11

b 水准点6

c 水准点12

d 水准点7

图9 水准点与MT-InSAR时间序列形变量对比

Fig.9 Comparison between time series displacement derived from levelling and from MT-InSAR

4 结语

本文介绍了一种多级化的MT-InSAR数据处理方法,提出了一种基于相干点空间密度的抽样方法,综合考虑选取相干点的密度以及空间分布,得到优化的分级MT-InSAR处理方法.并利用该方法获取了2007~2018年间的上海磁悬浮及其周边区域沉降信息,不同数据结果均显示磁悬浮轨道保持稳定,而周边区域存在不同程度的沉降,随着时间的累积,沉降速率是逐渐变小的.采用精密水准数据对沉降监测结果进行了精度验证,两者在时间以及空间上吻合度较高.MT-InSAR技术得到的形变速率可以满足高精度的地面形变监测要求.

此外,本文还分析了不同SAR影像数据分辨率对于磁悬浮沿线相干点获取的影响,发现相干点的密度与SAR数据的距离向分辨率有着密切关系,而方位向分辨率影响较小.对比上一代C波段中等分辨率ENVISAT数据,开源的Sentinel-1A数据有着更大的覆盖范围,并且凭借较高的距离向分辨率,可以获取更多的相干点,在城市地表以及大型结构体形变检测中有着突出的优势.