一种基于形态学梯度的反演城区建筑物沉降的方法

摘要

本发明涉及一种基于形态学梯度的反演城区建筑物沉降的方法，该方法考虑到建筑物在高分辨率的COSMO SkyMED雷达影像表现出强反射特性以及强相干性，首先基于波谱统计特性的相干系数阈值方法和基于后向散射强度的幅度阈值方法从高分辨率的SAR影像中提取出高质量的候选点目标，然后采用小基线集（SBAS）方法解算得出点目标的沉降信息；对InSAR解算得到的结果设置高程阈值，粗选建筑物上的点目标；利用形态学梯度方法提取的SAR影像中建筑物边界，并且利用提取的建筑物边界进行粗选点目标中非建筑点去除，最终实现建筑物点目标的自动提取。该方法提高了建筑物沉降信息识别的准确性和可靠性。

说明书

技术领域

本发明涉及一种基于形态学梯度的反演城区建筑物沉降的方法，属于遥感技术应用领域。

背景技术

随着我国经济的高速发展，城市化进程的加快，建筑事业也快速发展。这些高大建筑在给人们生活带来便利的同时也存在许多安全隐患。为了保证建筑物的安全，需要对建筑物进行沉降观测。

以InSAR时序分析技术为代表的方法大大增强了InSAR技术在城市地面沉降监测的能力。随着卫星遥感技术的发展，高分辨率的卫星影像如TanDEM-X，TerraSAR-X和Cosomo-SkyMed雷达影像相比较中低分辨率的雷达影像能够更好地捕获建筑物时空沉降特性，使得InSAR时序分析技术进行城市建筑物的精细形变监测成为可能。但是在实际应用中还存在一些问题，如不能够较好地将建筑物与背景地物沉降信息分离出来。针对上述问题，本发明采用了形态学梯度方法以SAR强度图中的建筑物为研究对象，简化图像数据，保持图像基本的形状特征，能够较好的提取建筑物的边缘。然后再结合InSAR时序分析技术监测的沉降结果将建筑物沉降信息从背景地物中分离出来，对于城市建筑物沉降监测具有重要的应用价值。

发明内容

本发明要解决的技术问题是：克服现有技术手段的缺点，提出一种基于形态学梯度的高分辨InSAR时序分析技术反演城区建筑物沉降的方法，有效的从高分辨率的InSAR城区形变监测结果中提取建筑物沉降值，大大提高了城区建筑物InSAR精细化监测水平。

为了解决上述技术问题，本发明提出的技术方案是：一种基于形态学梯度的反演城区建筑物沉降的方法，包括以下步骤:

第一步、选取高分辨率SAR数据；

第二步、获取SAR影像平均强度图，并且对高分辨率的SAR数据进行InSAR时序分析解算，获得地理坐标系下的地物目标点的高程改正值和地物目标点的地表形变速率信息；

第三步、提取高程改正值大于预设高程阈值的地物目标点，作为建筑物粗选点；

第四步、使用形态学梯度算法提取SAR影像平均强度图上的建筑物；

第五步、采用边缘检测算子对第四步提取的建筑物进行边缘检测，提取出建筑物边缘；

第六步、将第五步提取的建筑物边缘与第三步获得的建筑物粗选点相结合，识别出建筑物点目标，建筑物点目标的地表形变速率信息为最终的建筑物沉降信息。

本发明“一种基于形态学梯度的反演城区建筑物沉降的方法”，还具有如下改进：

1、所述第二步中，InSAR时序分析方法采用的小基线集SBAS方法(小基线集SBAS方法可参见文献"A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms,"IEEE T Geosci Remote,vol.42,pp.1377-1386,2004.和"An Overview of the Small BAseline Subset Algorithm:a DInSAR Technique for Surface Deformation Analysis,"Pure Appl Geophys,vol.164,no.4,pp.637-661,2007.)，该方法的具体实施步骤如下：

步骤S1、利用SBAS方法将所获得的SAR数据组成若干个集合；

步骤S2、对步骤S1集合内干涉对进行配准，对配准后的主辅影像复共轭相乘得到差分干涉图，采用自适应滤波去除差分干涉图的噪声相位，使用SRTM DEM去除差分干涉图的地形相位，采用最小费用流法进行干涉相位解缠；

步骤S3、采用基于波谱统计特性的相干系数阈值方法和基于后向散射强度的幅度阈值方法相结合提取永久散射体目标点；

步骤S4、对点目标的干涉相位经过回归分析，迭代求解地物高程改正值和地表形变速率信息；

步骤S5、将雷达坐标转化到地理坐标系下。

2、所述步骤S3中，基于波谱统计特性的相干系数阈值方法计算公式如下：

$\frac{1}{M}\underset{i=1}{\overset{M}{\Sigma }}{\gamma }_{i,j}\ge {\gamma }_T$

式中M为影像的个数，γ_i表示第j个像元在第i幅干涉图中的相干系数，γ_T为设定的相干系数阈值；

基于后向散射强度的幅度阈值方法的计算公式为：

D_i≤D_threshold

式中，D_i为点目标i的幅度离散指数，D_threshold为幅度离散指数阈值，D_threshold取的范围为0.25—2.00，即可选择符合要求的高相干的永久散射体候选点。

3、所述的第四步中,形态学梯度算法提取SAR影像强度图上建筑物，该方法的具体步骤如下：

步骤T1、SAR影像中像元梯度值计算，依据建筑物和周围地物过渡处的灰度值发生突变，采用梯度来表示灰度值突变的程度，将其作为建筑物提取的依据，采用形态学梯度为依据的高分辨SAR影像上建筑物边缘信息提取的算法，

根据如下公式计算像元的形态学梯度，

$d_f=(f\oplus b)-f$

式中d_f代表形态学梯度值，f代表待处理的SAR影像平均强度图中的像元强度值，b代表结构元素，代表膨胀运算，通过对像元点膨胀运算后的灰度值减去原来的灰度值得到该像元的梯度值；

步骤T2、对SAR影像进行条件腐蚀，如果SAR影像中像元k的梯度值d_f大于阈值t_threshold，则以像元k为中心，结构元素c的大小为尺寸对邻域内的像元做腐蚀，若所述邻域内的某像元灰度值大于像元k的灰度值，则该像元经过腐蚀后的灰度值不变，否则该像元经过腐蚀后的灰度值等于像元k的灰度值；

步骤T3、对上述处理得到的结果进行条件膨胀，用结构元素c做膨胀，如果像元q经腐蚀后的灰度值大于原始灰度值，则以像元q为中心结构元素c的大小为尺寸对邻域内的像元做膨胀，若所述邻域内的某像元灰度值小于像元q的灰度值，则经过膨胀后的灰度值不变，否则该像元经过膨胀后的灰度值等于像元q的灰度值；

步骤T4、计算每一个像元的原始灰度值与膨胀后的灰度值的差，如果大于阈值t_threshold，则判该像元为建筑物上的像元点；

步骤T5、对步骤T3得到的膨胀后的灰度图进行形态学开运算；

本步骤使得离散的像元点被去除，删除了不能包含结构元素的对象区域，且去除了细小的突出部分，使得建筑物的轮廓得到平滑；

步骤T6、在形态学开运算的基础上进行的形态学闭运算；

步骤T7、借助步骤T4中识别出的建筑物上的像元点，判断形态学闭运算后的结果中建筑物是否被识别，如果不能被识别，调整结构元素的种类和/或大小，重复步骤T2到T6，直到识别出所有的建筑物。

4、所述第五步中，首先对第四步处理得到的结果进行二值化处理，然后采用Canny算子提取二值化后图像中建筑物边缘。

本发明的有益效果是：

基于形态学梯度的高分辨InSAR时序分析技术反演城区建筑物沉降的方法对于研究城区建筑物沉降意义重大。本发明实现了基于形态学梯度的高分辨率InSAR技术监测城区建筑物沉降的方法，利用波谱统计特性的相干系数阈值方法，后向散射强度的幅度阈值方法，以及高程阈值方法，结合形态学梯度方法准确地获取了城区建筑物上的点目标，进而可以对建筑物上的点目标进行相关沉降分析，相较于传统方法具有无可比拟的优势。具体有益效果如下：

第一，本发明成功提取了城区的InSAR时序点目标，可进一步应用于建筑物沉降分析和预测建筑物沉降趋势。

第二，本发明采用的数据为COSM-SkyMed卫星影像，分辨率较高，用于建筑物的沉降分析精度较高。

第三，本发明利用地物的波谱统计特性、相干系数和高程的差异，使用基于后向散射强度的幅度阈值方法、基于的相干系数阈值方法和设定一定高程阈值方法实现建筑物目标的粗提取。然后利用形态学梯度的方法提取出建筑物边缘，最终将建筑物上的InSAR时序分析点目标提取出来。

第四，本发明中的各项算法简单有效(主要为InSAR时序分析计算和梯度形态学运算)，精度较高。整个过程无需人工介入，可应用于大面积的建筑物沉降分析。

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明数据处理流程图。

图2是COSMO-SkyMed雷达影像强度图，矩形区域内为本次试验区域。

图3是COSMO-SkyMed雷达影像数据集时间基线和空间基线关系图。

图4是InSAR时序分析结果中建筑物点目标粗提取结果。

图5是形态学梯度建筑物边缘提取结果。

图6是基于形态学梯度建筑物点目标提取结果。

具体实施方式

下面根据附图详细阐述本发明，使本发明的目的和效果变得更加明显。

本发明采用的示例影像数据为12景意大利航天局和意大利国防部共同研发的COSMO-SkyMed高分辨率雷达卫星(SAR)产品，空间分辨率为3m，影像的覆盖区域为江苏省常州市，时间跨度为2011年7月21到2012年7月14日。

图1为基于形态学梯度的高分辨InSAR时序分析技术反演城区建筑物沉降的方法的流程图，该方法的具体实施步骤包括以下步骤：

第一步、高分辨率SAR数据选取。

选取意大利COSMO-SkyMed星载X波段3米分辨率的2011年7月至2012年7月的12景SAR影像为研究数据。图2为SAR影像强度图，矩形框内为选取的建筑物分布示例区域。COSMO-SkyMed的扫描带宽为10km，轨道周期为16天。Cosmo-Skymed雷达影像的主要数据参数如表1所示。

表1

第二步、获取SAR影像平均强度图，并且对高分辨率的SAR数据进行InSAR时序分析解算，获得地理坐标系下的地物目标点的高程改正值和地物目标点的地表形变速率信息。

InSAR时序分析方法采用的小基线集SBAS方法，具体步骤如下：

步骤S1、利用SBAS方法将所获得的SAR数据组成若干个集合。SBAS方法不受长基线导致的几何去相干的影响，使得更多的图像参与计算，增加时间的采样，在本发明中小基线数据集的构成上，采用时间基线T小于180天，空间基线B小于800m，如图3所示。

步骤S2、干涉图的生成。通过步骤S1得到的干涉对需要共轭相乘生成干涉图，而干涉对中的两幅SAR影像由于入射角的差异使得两幅影像不可能完全重合，因而需要对SAR影像配准处理后，本发明中配准值采用0.125个像元。对配准后的主辅影像复共轭相乘得到差分干涉图。本发明采用自适应滤波去除干涉图的噪声相位。使用SRTM DEM去除地形相位。采用最小费用流法MCF(Minimum Cost Flow)进行相位解缠。

步骤S3、候选点目标识别。采用幅度离散指数为1.7，相干系数阈值为0.8共同选取相干点目标候选点(提取永久散射体目标点)。

步骤S4、相干点目标形变结果计算。通过对以上设定的幅度离散指数和相干系数阈值共同得到点目标候选点，建立点目标的二维回归分析模型，迭代求解，最后得到的地物目标点的高程改正值和形变速率值信息。

步骤S5、地理编码。采用地理编码将以上处理得到的雷达坐标系下的测量值转化到地理坐标系下。

第三步、采用相关高程阈值粗略提取建筑物上的点目标，本实验中设定阈值h_threshold为10m。提取的研究区域的建筑物的点目标分布。图4表示SAR影像处理后提取的建筑物粗选点目标的信息。

第四步、使用形态学梯度算法提取SAR影像平均强度图上的建筑物。

该步骤的具体实施如下：

步骤T1、SAR影像中像元梯度值计算。计算建筑物和周围地物过渡处的灰度发生突变的值为梯度值。

步骤T2、对SAR影像进行条件腐蚀。对上述步骤中计算得到的梯度值设定相关阈值，此处设定为SAR影像的灰度均值，以大小为3×3的结构元素进行条件腐蚀。

步骤T3、用同样大小的结构元素作条件膨胀。

步骤T4、计算每一个像元的原始灰度值与膨胀后的灰度值的差，如果大于阈值t_threshold，则判该像元为建筑物上的像元点。

步骤T5、对步骤T3得到的结果进行形态学开运算。

步骤T6、在开运算得到结果的基础上进行闭运算。

步骤T7、判断形态学闭运算后的结果中建筑物是否被识别，如果不能被识别，调整结构元素的种类和/或大小，重复步骤T2到T6，直到识别出所有的建筑物。

第五步、在提取建筑物边缘上，本实验对上述形态学梯度运算后得到的结果做二值化转换，二值化转换的阈值为像元灰度的平均值，然后采用canny算子对二值化的结果做边缘检测。结果如图5所示，可以看出建筑物边缘基本被提取出来。

第六步、将第五步提取的建筑物边缘信息结合InSAR时序分析结果中建筑物粗选点目标信息进行叠置分析计算，去除非建筑点目标，得到最终建筑物点目标的沉降信息。结果如图6所示，可以看出相比较InSAR粗提取的建筑物点目标信息，经过形态学梯度的检测后，最终得出的建筑物上点目标分布与建筑物分布一致，结果较为准确。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。