1 引言
GB-InSAR(Ground-based Interferometric Synthetic Aperture Radar,地基干涉合成孔径雷达)是近些年受到广泛关注的形变测量技术,具有全天时,高分辨,高精度等优势,其雷达成像范围可以达到数平方公里[1]。GB-InSAR的测量原理与星载SAR相同,均是通过雷达对同一片区域进行多次观测,基于相位差分干涉技术实现地表的形变测量。星载SAR受到监测角度与访问权限的影响,具有较长的重访周期,限制了监测的实时性;相比而言,GB-InSAR可以实现全天时、全天候的监测,图像获取仅需几分钟,因而能够实现对观测目标的近实时动态连续监测[2]。
在差分干涉测量中,受到时间、空间去相关的影响,雷达图像中大量像素点的相位稳定性很低,如植被、水体等区域的散射特性在短时间内会发生较大的改变,对地表形变的准确测量带来很大的干扰[3]。现阶段广泛采用PS技术来实现差分干涉测量,该技术由意大利的Ferretti等学者提出。通过计算多幅雷达图像中,每一个像素点幅度序列的标准差与均值的比值,定义为幅度离差值,然后设定合理的幅度离差门限,即可实现PS点的选择[4]。该方法实现简单,在利用GB-InSAR系统进行长时间形变监测时,由于GB-InSAR系统的图像获取速度很快,每天可获取到上百幅雷达图像。为了保证形变测量的实时性和准确性,现阶段常采用PS分组选择的方式,每获取到一幅新图像,即结合前面一定数量的历史图像,采用幅度离差法进行一次PS点的选择,然后基于差分干涉处理,实现对当前新图像的形变测量。
对时序GB-InSAR图像的分析结果表明,在长时间的监测中,PS点会出现消失和新生的现象,对于不同的SAR图像组,像素点质量存在着较大差异。在形变反演的处理过程中,是基于PS点的形变结果,利用空间维插值的方法,对整个观测场景进行形变估计,如果沿用统一的门限,会使得分组PS点选择的数量波动较大,导致插值过程中出现误差,影响形变反演的精度[5]。
针对上述的问题,本文基于K-Means(K均值)聚类算法,先后以幅度序列和相关系数序列作为输入,进行两次分类,进而挑选出符合要求的高质量像素点。相比于门限方法,聚类方法可以适应幅度分布在不同场景和时间上的变化,因此在长时间的连续监测中,利用这种方法可以保持PS点数量和质量的稳定。本文分为两部分,首先介绍了实验信息与PS选择算法流程,然后分别用常规方法和改进方法对一处露天矿坑的1650幅SAR图像进行了分析,通过对比验证了本文方法的有效性和优势。
2 实验信息
2.1 实验场景介绍
本实验地点位于河北省迁安市马兰庄镇,是一处露天的开采矿坑(E118°36′23″,N40°06′44″),图1(a)所示为场景卫星图,其中弧线标识出了雷达监测区域,矩形框标注了雷达的布设位置,图1(b)为场景照片。雷达的测量距离约为400～900 m,测量视角范围为90°,矿坑的深度约为400 m,其边坡主要由硬质岩石构成,基本无植被覆盖[6]。
图1 实验场景信息
Fig.1 Experimental scene information
实验选用的雷达为MIMO干涉雷达,图2所示为MIMO雷达照片,发射信号为Ku波段的线性调频连续波,其波长λ约为18.5 mm,雷达采用16个发射天线和32个接收天线,合成具有512个相位中心的等效线性阵列[7]。相邻发射天线的间隔为λ/2=9.3 mm,相邻接收天线的间隔为8λ=14.8 cm,表1给出了MIMO雷达测量系统参数[8]。
表1 MIMO雷达系统参数
Tab.1 MIMO radar system parameters
技术指标参数值技术指标参数值载波频率16.2 GHz合成孔径长度4.752 m波长18.5 mm距离向分辨率0.375 m带宽400 MHz方位向分辨率1.948 mrad采样率25 MHz探测范围30～3000 m
图2 MIMO干涉雷达照片
Fig.2 MIMO interferometric radar photo
2.2 PS实时处理
2.2.1 PS门限选择
在利用GB-InSAR进行形变测量时,由于时间和空间去相关的影响,像素点的相位信息存在误差,会影响形变反演的精度。Ferretti等人提出了PS技术,以多幅SAR图像为输入,选取出幅度特征高度稳定的像素点作为PS点。该算法是以幅度离差作为选择指标[9]:
(1)
式(1)中,σA为像素点幅度序列的标准差,mA为像素点幅度序列的均值,通过设立合适的门限DT,满足DA<DT条件的像素点即为PS点[10]。但只有在高信噪比(SNR)的条件下,幅度离差才可以准确的表征相位稳定性,因此仅仅通过幅度离差选择是不准确的,实际处理中需要进一步结合幅度或相干系数门限进行PS选择[11]。
为了更好地说明PS选择对形变处理结果的影响,对马兰庄露天矿坑的30幅图像数据进行处理,结果如图3所示。其中图3(a)和图3(b)所示分别为雷达成像结果和幅度离差信息,可以看出,在中间的坡体区域,像素点的幅度离差主要分布在0～0.2之间,幅度稳定性较高;图3(c)是时间基线为3 min的一幅干涉相位图,由于时间基线较短,可以认为场景并未发生形变,大气相位的影响也很小,因此干涉相位应该接近于0,但岩坡周围出现了较为明显的误差相位,说明这些像素点的相位信息并不稳定。设置幅度离差门限为0.1,筛选出的PS点如图3(d)所示,坡体大部分区域较为稳定,仅中间的道路由于过往车辆和雷达视角等原因未被选出。
图3 场景实测数据
Fig.3 Scene measured data
2.2.2 SAR图像分组
本实验利用MIMO雷达对马兰庄露天矿坑进行了监测,连续获取了1650幅SAR图像,时间从2019年6月30日16时53分到2019年7月5日20时30分。为了减弱长时间基线对反演精度的影响,对SAR图像进行分组处理,由于组内的时间基线较短,因此每幅图像的PS点集合可以认为是稳定的。图4显示了对SAR图像进行PS分组选择的流程,其中N为每组图像的数目。当第N幅图像生成后,基于第1～N幅进行第1次PS选择;从第N幅图像后,每获取一幅新的图像,都结合前N-1幅进行PS选择。通过分组处理,PS选择算法的实时性有了较大的提升,但在形变反演的处理中,是通过PS点的形变量估计观测区域的形变信息,如果各组图像中PS集合的差异较大,会影响最终结果的精度,因此如何保证每组PS点数量和质量的稳定性是一个需要解决的问题。
图4 PS分组处理示意图
Fig.4 PS grouping process diagram
2.2.3 PS门限选择结果
GB-InSAR图像的获取时间较短且图像数量较多,在实际处理中为了减少运算量,可以每获取N幅图像进行一次PS选择。对河北省迁安市马兰庄矿坑的1650幅SAR图像进行分组PS实时处理,一共分为55组,每组30幅,设立幅度离差门限0.1。PS分组选择结果如图5所示,横坐标为日期轴,虚线表示各组PS点数量的均值。由图可知,PS分组数量随时间有着较强的周期性变化,其中第16组的PS点最多,达到了123215个;第23组的PS点最少,仅为2630个,分别称这两组图像为A组和B组。可以看出,各组PS点的数量波动很大,这对最终的形变反演结果会造成较大影响。
图5 PS点数量变化
Fig.5 Changes in the number of PSs
对A组和B组图像进行分析, 图6(c)和图6(d)是两组图像的PS选择结果,图6(e)和图6(f)则显示了这两组幅度离差的概率密度曲线和概率分布曲线。由图可知,A组的幅度离差指数明显比B组要小,因此如果采用固定的幅度离差门限,会造成两组PS选择结果的较大差异。
图6 两组图像幅度离差对比
Fig.6 Comparison of amplitude dispersion
3 改进方法与结果分析
3.1 改进方法
3.1.1 K-Means聚类原理
聚类就是按照某个特定标准(如距离准则)把一个数据集分割成不同的类或簇,与分类算法不同,聚类技术常被称为无监督学习[12]。K-Means聚类算法运用了迭代的思想,事先确定类别数k和聚类中心,通过计算每个样本和质心之间的欧氏距离,将样本归纳到最相似的类中,并重复这样的过程,直到质心不再改变,最终就确定了每个样本所属的类别[13]。
假定有n个数据样本x={xi},i=1,2,3,…,n,其中每个样本xi有m个维度的特征值,然后指定k个类簇中心{C1,C2,…,Ck},根据式(2)将每个样本标记为距离类别中心最近的类别:
(2)
其中||·||表示2阶范数,即向量的模长,根据分类结果更新聚类中心,可以表示为式(3):
(3)
不断地重复这两个步骤,直到聚类中心不再发生变化,最终生成k个簇,即k个分类结果,K-Means聚类算法的示意图如图7所示,设置k=4,其中图7(a)显示了数据样本点的分布,其中黑色点代表初始类簇中心;图7(b)给出了聚类后的结果,四种形状的点各为一簇,黑色点为更新之后的类簇中心。
图7 K-Means算法示意图
Fig.7 K-means clustering algorithm
3.1.2 PS两级选择
根据2.2.3节的结果可知,SAR图像的幅值具有时变性,在分组处理时,不同组图像的时间基线较长,场景像素点的幅值分布会有差异,因此会造成PS数量的不稳定。相比于传统的门限选择,K-Means聚类算法是根据样本特征和已知类别进行优化迭代,因此对于不同图像组中像素点的输入序列,该方法有着较强的自适应性。
本文提出了一种基于K-Means聚类的PS两级分类方法,其流程如图8所示,对于每一组图像的M个像素点集合,以他们的幅度时间序列Sk={sk1,sk2,…,skN}(其中k=1,2,…,M,N为每组图像数量)作为输入,划分为两类像素点,将高质量的像素点作为PS候选点。如果仅根据幅度信息进行分类,可能会导致最后的结果存在相位误差,因此为了提高形变反演的质量,需要对PS候选点进行二次分类。
图8 基于K-Means聚类的PS两级选择方法流程图
Fig.8 PS two-level selection method based on K-means clustering
相关系数是衡量相位误差的重要手段,为了解决上述问题,可以利用像素点的相关系数序列γk={γk1,γk2,…γkN-1}对PS候选点进行二次分类,在计算每个像素点的相关系数时,先以这个点为中心,在主辅图像中分别构造d1×d2的矩形窗,然后依据式(4)进行计算:
(4)
式(4)中,
代表在第i幅SAR图像中,以第k个像素点为中心的矩形窗所包含的像素点信息,计算相关系数时,选用的主辅图像是相邻的,因此每个相关系数序列的维度是N-1,利用相关系数信息进行二次分类后,将高质量的像素点选为最终的PS集合。
3.2 实验分析
现根据3.1节所提出的PS分类方法,对A组SAR图像进行处理,首先以每个像素点的幅度序列Sk={sk1,sk2,…,sk30}(其中k=1,2,…,M,M为像素点个数)作为输入进行第一次分类。根据K-Means聚类原理,在二分类问题中,每轮迭代运算会生成两个幅度中心{Ci1,Ci2,…,Ci30}(其中i=1,2),并不断更新直到中心不再变化,实验中对输入数据进行了25次迭代,将靠近更大幅度中心的像素点选为PS候选点,根据这种方法一共选出了78267个候选点,占总像素点比例约为6.32%。图9(a)表示了PS候选点和非PS点的平均幅度序列,由图可知,PS候选点的整体幅度明显更大。图9(b)显示了PS候选点的归一化幅值,大部分的矿坑边坡已被选为PS候选点。
在第二次分类前,以每两幅相邻的SAR图像为主辅图像,根据式(4)生成29幅相关系数图,然后根据第一次K-Means聚类的结果,将每个PS候选点的相关系数序列γk={γk1,γk2,…,γk29}(其中k=1,2,…,Mc,Mc为PS候选点个数)作为输入进行二次聚类,一共进行了13次迭代,将相关系数较大的一类像素点选为PS点,通过这种方法,一共选出了75480个PS点,占PS候选点数量的96.44%,PS点的相关系数如图10(b)所示。图10(a)表示了PS点与非PS点的平均相关系数序列,由图可以看出,PS点具有更好的相关性与稳定性。
图9 第一次K-Means聚类结果
Fig.9 The first K-means clustering result
图10 第二次K-Means聚类结果
Fig.10 The second K-means clustering result
分别用传统门限方法与改进方法,对A组和B组图像进行PS选择,表2显示了两种方法选出的PS点数量,可以看出当使用改进的PS选择方法处理时,两组图像的PS点数目更加接近,并且对于B组图像,改进方法选出的PS点数量有了显著提高。图11(a)～(d)分别显示出了A组和B组图像使用两种方法进行PS选择的干涉相位图。
表2 PS选择方法结果对比
Tab.2 PS selection method results comparison
PS选择方法A组PS点数量占像素点总数的比例B组PS点数量占像素点总数的比例传统门限方法1232159.95%26300.21%改进的两级选择方法754806.10%683215.52%
分析这两组图像PS点的质量。根据差分干涉的基本原理,PS点干涉相位Δφ的主值范围是[-π, π]。为了实现高精度的形变测量,需要对干涉相位进行解缠处理,并准确地估计和补偿大气相位分量,进而才能得到准确的形变相位。首先对于每组的30幅SAR图像,利用门限方法与改进方法分别选出对应的PS点集合,然后以第1幅为主图像,其他图像为辅图像,生成29幅干涉相位图,并基于PS点集合进行相位解缠和大气相位补偿,最终得到对应的形变相位序列。由于30幅SAR图像的时间跨度仅为80分钟左右,观测区域为稳定的岩质边坡,可以认为在观测时间内几乎没有发生形变,每个PS点经过大气补偿后的干涉相位主要由噪声分量构成,因此可以用大气补偿后的相位序列标准差来表示PS点的稳定性。
图12显示了A组和B组图像中PS点相位标准差的分布情况,其中图12(a)和(c)分别表示这两组PS点相位标准差的概率密度曲线;图12(b)和(d)所示为对应的概率分布曲线。由图可知,对于B组图像,使用门限方法和改进方法选出的PS点集合,其相位标准差存在着差异,统计出标准差在不同区间上的分布如表3所示。当雷达工作在Ku波段(工作频率为16.2 GHz),0.8 rad的相位标准差对应着1.179 mm的形变测量精度,由表可知,对于两个PS点集合,相位标准差小于0.8 rad的占比分别为98.48%和96.40%,因此如果以大气相位补偿后的标准差作为相位质量的评判依据,改进方法可以在保证相位质量的前提下,选出更多的PS点。
图11 实测数据结果
Fig.11 Measured data results
图12 PS点相位标准差
Fig.12 Phase standard deviation of PSs
表3 PS点相位标准差分布(B组)
Tab.3 Phase standard deviation distribution of PSs (Group.B)
标准差/rad门限方法数量占比/%改进方法数量占比/%<0.3149957.002180031.91<0.4214181.413979258.24<0.5242992.365278377.26<0.6254596.776018788.09<0.7258198.146405893.76<0.8259098.486586096.40
图13 PS分组选择结果对比
Fig.13 Comparison of PS group selection results
利用改进方法对全部的55组SAR图像进行PS选择,数量结果如图13所示。由图可知,每组PS点的数量稳定性有很大的提升,随时间的波动也更加轻微,通过实测数据处理结果的对比,进一步说明了改进的PS选择算法有着较强的自适应性,在长时间连续监测中,无需人为设定和调整门限。
4 结论
本文基于GB-InSAR技术,对形变监测中永久散射体(PS)的选取展开了研究,针对传统门限方法所存在的不足,提出了一种基于K-Means聚类原理的PS两级分类方法。在实验中用MIMO干涉雷达对一处露天矿坑进行观测,连续获取了多幅SAR图像数据,分别利用传统的幅度离差阈值选择法和改进方法,对SAR图像序列进行分组PS选择。实测数据结果表明,改进方法能够有效提高PS点数量和分组选择的稳定性,具有较强的自适应性。本文方法也存在着一些不足,如对于输入序列只是做了简单的二分类,当更换场景时,类别个数K可能需要重新设定;对植被较多的场景,PS点的质量会有一定下降;当像素点很多时,运算量会比较大等,这些问题也需要在后续进一步的研究。
[1] 龙四春, 陈鹏琦, 袁英, 等. GB_InSAR技术及其形变监测[J]. 测绘通报, 2015, 9: 1-5.
LONG Sichun, CHEN Pengqi, YUAN Ying, et al. Ground-based Synthetic Aperture Radar Interferometry and Its Deformation Monitoring[J]. Bulletin of Surveying and Mapping, 2015, 9: 1-5.(in Chinese)
[2] 邹进贵, 张士勇, 李琴, 等. 基于GBSAR的变形监测方法综述[J]. 测绘地理信息, 2016, 41(6): 5- 8.
ZOU Jingui, ZHANG Shiyong, LI Qin, et al. Review of Ground-Based SAR Interferometry for Deformation Measurement[J]. Journal of Geomatics, 2016, 41(6): 5- 8.(in Chinese)
[3] 曾涛, 邓云开, 胡程, 等. 地基差分干涉雷达发展现状及应用实例[J]. 雷达学报, 2019, 8(1): 154-170.
ZENG 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. (in Chinese)
[4] ALESSANDRO Ferretti, CLAUDIO Prati, FABIO Rocca. Permanent scatterers in SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8-20.
[5] 朱茂. 基于动态PS的地基合成孔径雷达高精度形变测量技术研究[D]. 北京: 北京理工大学, 2016.
ZHU Mao. High Precision Deformation Measurement Using Ground Based Synthetic Aperture Radar (GBSAR) Based on Dynamic Persistent Scatterer (PS) Technique[D]. Beijing: Beijing Institute of Technology, 2016.(in Chinese)
[6] 胡程, 邓云开, 田卫明, 等. 地基干涉合成孔径雷达图像非线性大气相位补偿方法[J]. 雷达学报, 2019, 8(6): 831- 840.
HU Cheng, DENG Yunkai, TIAN Weiming, et al. A Compensation Method of Nonlinear Atmospheric Phase Applied for GB-InSAR Images[J]. Journal of Radars, 2019, 8(6): 831- 840.(in Chinese)
[7] 邓云开. 地基差分干涉雷达形变测量理论与方法研究[D]. 北京: 北京理工大学, 2019.
DENG Yunkai. Researches on Theory and Method of Ground-Based Differential Interferometric Radar Measurement[D]. Beijing: Beijing Institute of Technology, 2019.(in Chinese)
[8] HU C, DENG Y, TIAN W, et al. Novel MIMO-SAR system applied for high-speed and high-accuracy deformation measurement[J]. The Journal of Engineering, 2019, 2019(20): 6598- 6602.
[9] SABINE Rödelsperger. Real-time Processing of Ground Based Synthetic Aperture Radar(GB-SAR) Measurements[D]. Darmstadt: Technische Universität Darmstadt, 2011.
[10] MONSERRAT M, CROSETTO, GUIDO Luzi. A review of ground-based SAR interferometry for deformation measurement[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 93: 40- 48.
[11] 陈强. 基于永久散射体雷达差分干涉探测区域地表形变的研究[D]. 成都: 西南交通大学, 2006.
CHEN Qiang. Detecting Regional Ground Deformation by Differential SAR Interferometry Based on Permanent Scatterers[D]. Chengdu: Southwest Jiaotong University, 2006.(in Chinese)
[12] 段明秀. 层次聚类算法的研究及应用[D]. 长沙: 中南大学, 2009.
DUAN Mingxiu. Research and Application of Hierarchical Clustering Algorithm[D]. Changsha: Central South University, 2009. (in Chinese)
[13] 刘涛, 尹红健. 基于半监督学习的K均值聚类算法研究[J]. 计算机应用研究, 2010, 27(3): 913-916.
LIU Tao, YIN Hongjian. Semi-Supervised Learning Based on K-Means clustering algorithm[J]. Application Research of Computers, 2010, 27(3): 913-916.(in Chinese)
