CN113139994A - 基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法 - Google Patents

Abstract

本发明公开一种基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，首先通过合理选择卫星成像参数提高卫星图像的可用性，其次利用阴影在多波段数据中的特点综合计算获取阴影增强图像数据，然后利用图像中典型角点通过角点检测算法获得准确的角点位置信息，利用空间分布进行角点配对，提高阴影长度的测量精度，最后利用阴影成像的几何模型计算建筑物的高度。本发明方法通过优化卫星成像参数提高图像可用性，通过多波段数据融合处理对遥感图像进行增强，并充分利用建筑物和阴影图像中两个边缘方向形成的角点特征实现建筑物高度的提取，解决了阴影长度测量不准确的问题，提高了建筑物高度监测的精度。

Description

基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法

技术领域

本发明属于光学遥感卫星影像数据处理和应用技术领域，具体涉及一种基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法。

背景技术

近年来，卫星遥感技术得到快速发展，商业光学遥感卫星的分辨率已经达到0.3米，市场上优于0.5米分辨率的光学遥感卫星数量已经超过10颗，高空间分辨率卫星遥感影像以直观准确、更新迅速和大区域范围等优点在动态监测、数据采集和信息分析等方面不断拓展新的应用。随着我国经济建设的高速发展，在我国城市和广大农村地区各种建筑物的变化日新月异。建筑物的高度信息对城市规划、城市建设和资源管理具有重要作用，同时各地行政区域范围内在建的建筑物施工进展和主体建设高度等信息，对于分析我国各地的固定资产投资和经济活跃度具有重要参考意义。因此研究一种实用、高效和高性价比的建筑物高度的监测方法十分必要。

目前常用的建筑物高度监测方法可以分为以下几种：①地面测量法，利用GPS、全站仪、激光测距仪等测量工具进行建筑物高度提取，此类方法技术已经非常成熟，适用于城市区域内已经建设完工并投入使用的建筑物，具有较高的测量精度。但是需要在建筑物测量点布设仪器开展外业工作，而且只能分别逐个测量建筑物高度，同时还存在观测周期较长和观测成本较高等问题，因此无法实现区域范围内快速的建筑物高度提取，对于建设施工工地中建筑物的施工进展情况也无法及时监测。②航空测量法，随着无人机技术的发展，基于航空平台结合光学摄影测量技术或激光雷达技术，通过立体像对和点云等处理技术，获得监测区域的DSM，从而得到较高精度的地面建筑物高度。但是同样也存在作业流程繁琐、数据处理复杂，而且效率低下等问题，尤其是很多城市和核心地区还会受到空域管制等限制，监测数据的更新频次受到制约，无法及时对于建筑物情况进行监测。③航天测量法，基于高分辨率卫星遥感数据获取建筑物高度信息的方法可以分为两类，包括立体像对方法，即遥感卫星采用立体测量模式，运用视差原理来获取建筑物高度的立体测量方法，目前已经有“高分七号”等多颗亚米级光学卫星具备立体成像能力，但是该方法对于卫星成像数据质量要求较高，数据处理非常复杂并且数据价格昂贵；基于阴影的提取方法，即利用单幅遥感影像中检测出的建筑物阴影长度，结合成像时建筑物与卫星、太阳等存在的相对空间几何位置关系，可以反演出建筑物的高度信息，该方法在技术上已经比较成熟，也得到较为广泛的应用。

建筑物阴影是由于建筑物本身对太阳光的遮挡，在背光方向形成较暗区域，因此在光学遥感图像上形成局部颜色较黑的区域，从数据上表现为各个光谱波段的DN(DigitalNumber)值会明显低于周围非阴影区域。根据研究，建筑物阴影在遥感图像上的总长度、太阳高度角、太阳方位角与建筑物高度存在一定的几何关系，而且受卫星成像视角的原因，图像中的阴影大小位置还受到卫星高度角和方位角的影响。现有的基于阴影的高分辨率遥感影像建筑物高度估测方法，主要步骤是先利用遥感图像中的光谱特性或空间影像纹理特征对阴影区域进行检测和分类，并根据阴影成像方向统计确定建筑物阴影的长度，最后利用方位角和高度角等信息估测建筑物高度。但是从现有研究成果分析，该方法还主要存在以下问题：(1)阴影检测，阴影的形成受到很多因素影响，并且阴影本身会叠加到地物上，因此仅依靠纹理特征、边缘特征和光谱特征等进行分类识别，会存在较多误差导致阴影区域分类精度不高。(2)阴影长度确定，由于卫星成像轨道、角度、时间以及建筑物的分布和高度等，使得影像中建筑物阴影可能存在遮挡等现象，阴影长度不完全是建筑物阴影投影的总长度；(3)阴影边缘的确定，很多情况下建筑物的阴影边界并不均匀，建筑物与阴影的边界、阴影与地表的边界存在模糊，因此提取的边缘轮廓并不是均匀平滑的曲线。现有大部分方法都是根据阴影方向通过测量阴影两个边缘线的距离确定阴影长度，即使通过统计筛选也容易产生较大的误差。(4)部分建筑不规则或者高度并非一致时候，测量沿阴影方向的两个边缘之间距离，再通过统计阴影距离长度的峰值来计算建筑物高度的方法也存在较大的不确定性。

发明内容

针对上述问题，本发明提出一种基于高分辨率光学卫星遥感影像的建筑物高度的快速监测方法，基于高分辨率光学卫星遥感影像中建筑物及其阴影间的空间几何特性和两个边缘方向的角特征，发挥遥感卫星大幅宽和高分辨率的优势，实现大范围区域内建筑物高度的监测，大型基建工程施工进展情况以及城市建筑群三维立体模型的构建，具备较高的测量精度和较强的可操作性，具有较好的市场价值。

本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，具体步骤如下：

步骤一、选择亚米级高分辨率光学卫星对待监测建筑物进行成像。

步骤二、获取包含待监测区域的原始高分辨率光学卫星图像产品，根据卫星图像产品参数文件获得太阳高度角和方位角，以及卫星的高度角和方位角参数，从而确定待监测建筑物的阴影的成像方向。

步骤三、对原始高分辨率光学卫星图像进行辐射定标、大气校正、全色和多光谱数据融合处理、几何校正以及裁剪，得到高分辨率的待监测区域图像。

步骤四、对高分辨率的待监测区域图像的三波段数据进行增强处理，利用待监测建筑物的阴影在图像中的灰度值较低的特性进行多波段融合的阴影增强处理，得到待监测建筑物的增强阴影图像。

步骤五、根据步骤二确定的待监测建筑物的阴影成像方向，以及步骤四得到的待监测建筑物的增强阴影图像上分别裁切建筑物及其阴影在图像中对应的典型角点区域；对裁切区域的影像角特征进行提取，分别计算获得角点的精确坐标位置，分别得到待监测建筑物及其阴影的角点配对。

步骤六、对步骤五获取角点对配对结果，分别计算欧式距离，再根据卫星参数的几何模型计算建筑物的高度。

本发明的优点在于：

(1)本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法中，提出基于卫星成像参数的建筑物及阴影遥感图像优化方法，通过选择高分辨率光学卫星以及设置合理的卫星成像参数，优化了获取数据的图像质量，降低了数据处理的复杂性，并且有效提高建筑物高度监测的可靠性。

(2)本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法中，提出多波段高分辨率光学遥感图像典型特征增强方法，通过全色和多光谱数据融合处理，获取保持高分辨率和光谱特性的真彩色影像，采用多波段数据的融合处理算法对阴影特性进行增强，显著提升了建筑物图像的边缘特性，为角特征提取奠定基础。

(3)本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法中，提出基于角点欧式距离的建筑物高度计算方法，通过统一的阴影构象几何模型，简化了卫星和太阳不同方位角度情况下的高度计算模型，并对角点进行筛选配对，并且可以监测建筑物不同位置高度，提高测量方法的可操作能力。

附图说明

图1为本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法流程图。

具体实施方式

下面结合附图对本发明进行详细说明。

本发明基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，如图1所示，具体步骤如下：

步骤一、选择亚米级高分辨率光学卫星对待监测建筑物进行成像，根据建筑物情况、卫星成像时间、卫星轨道和成像角度等选择合理的高分辨率光学遥感卫星；

从各卫星数据供应商处挑选分辨率优于0.5米的亚米级高分辨率光学卫星，然后根据待监测建筑物的经纬度位置和监测季节，获得该亚米级高分辨率光学卫星的过境成像时间和侧摆角度，大致确定卫星对建筑物成像时刻的太阳高度角和方位角，优先选择待监测建筑物阴影不存在遮挡以及小侧摆成像的卫星；最后对该卫星进行编程成像确保数据高可用性。

步骤二、获取包含待监测区域的高分辨率光学卫星图像产品，根据卫星图像产品参数文件获得太阳高度角和方位角，以及卫星的高度角和方位角参数，从而确定待监测建筑物的阴影的成像方向。

首先，获取卫星数据商提供的原始高分辨率光学卫星图像产品，以SuperView卫星数据为例，其产品包括一个0.5米分辨率的全色波段图像数据及辅助数据文件、2米分辨率的4个多光谱波段图像数据及辅助数据文件等。

然后，从辅助数据文件中分别读取太阳高度角β、太阳方位角θ、卫星高度角α和卫星方位角

参数。

最后，根据太阳方位角和高度角，以及卫星图像特性，确定待监测建筑物的阴影成像方向。

步骤三、对原始高分辨率光学卫星图像进行辐射定标、大气校正、全色和多光谱数据融合处理、几何校正以及裁剪，得到高分辨率的待监测区域图像。

首先，使用步骤二获取的原始高分辨率光学卫星图像，对图像进行辐射标定和大气校正处理。然后将全色图像和多光谱图像进行波段配准和融合处理，得到分辨率为0.5米的融合后的真彩色影像，其包括了RGB三个波段数据。最后，根据待监测区域范围对融合后的真彩色影像进行裁切，将待监测建筑物的图像提取出来。

步骤四、对高分辨率的待监测区域图像的三波段数据进行增强处理，利用待监测建筑物的阴影在图像中的灰度值较低的特性进行多波段融合的阴影增强处理，得到待监测建筑物的增强阴影图像。

首先，分别计算RGB三个波段数据灰度值的均值和方差值；利用阴影的灰度值在各个多光谱波段中都远远低于正常地物灰度值的特性，计算三个波段组合的增强图像S_en，计算公式如下:

其中，B表示蓝色波段的灰度值，G表示绿色波段的灰度值，R表示蓝色波段的灰度值，μ_B表示蓝色波段的灰度值的均值，μ_G表示绿色波段的灰度值的均值，μ_R表示红色波段的灰度值的均值；σ_B表示蓝色波段的灰度值的方差，σ_G表示绿色波段的灰度值的方差，σ_R表示红色波段的灰度值的方差。

步骤五、根据步骤二确定的待监测建筑物的阴影成像方向，以及步骤四得到的待监测建筑物的增强阴影图像上分别裁切建筑物及其阴影在图像中对应的典型角点区域；对裁切区域的影像角特征进行提取，分别计算获得角点的精确坐标位置，分别得到待监测建筑物及其阴影的角点配对。

步骤501，根据步骤二确定的待监测建筑物的阴影成像方向，在待监测建筑物的增强阴影图像上裁切待监测的建筑物楼顶的角点区域图像I_i(x,y)。其中，(x,y)为像素点坐标，i为选取裁切的角点区域图像数量。建筑物楼顶的角点区域通常是靠近阴影部分的楼角或者楼顶的突出建筑物，裁切大小一般为36*36像素。

步骤502，同上，在待监测建筑物的增强阴影图像上裁切待监测建筑物楼顶的角点区域图像I_i(x,y)对应的阴影区域的角点图像I_i′(x,y)，通常是远离建筑物方向的阴影角点，裁切大小一般为36*36像素。

步骤503，计算图像I(x,y)在水平和垂直方向的一阶偏导函数I_x和I_y，利用高斯函数h(x,y,σ)作为窗口对协方差矩阵

进行平滑滤波，获得像素点的自相关矩阵M，所述

其中，(u,v)表示其它像素点相对于像素点(x,y)的偏移量；ω_u,v为高斯加权窗口函数；此处高斯函数表示为：

其中，σ为正态分布的标准偏差，确定函数的作用范围，σ越大影响的样本越多。

步骤504，计算角点响应函数R，R＝det(M)-k×trace²(M)，其中，det(M)为自相关函数矩阵M的行列式，trace(M)为自相关矩阵M的迹，k是可调参数，一般设置范围为(0.04-0.06)。

步骤505，通过R的值大小来判断像素点是否为角点，对R进行归一化排序，设置阈值T进行判读是否为角点。其中，阈值T是可调参数，本发明中设置T为0.8，只有当R＞T时认为是角点，记录其对应的角点像素坐标(x_i,y_i)_j，j为检测到的角点数量。

步骤506，对阴影区域的角点图像I_i′(x,y)重复步骤503、504和505，得到阴影区域的角点坐标(x_i',y_i')_j，j是检测到的角点数量。

步骤507，对步骤505和步骤506得到的角点进行筛选配对，首先剔除不在阴影和建筑物边缘以及孤立的无效角点，再分别计算(x_i,y_i)_j和(x_i',y_i')_j中剩余角点对连线的角度，根据卫星空间几何模型获取的方位角ω_i，筛选角度一致的为配对角点对(x_i,y_i)_n和(x_i',y_i')_n，其中n是筛选配对后得到的角点对数量，

所述

步骤508，对步骤501和502裁切的i个角点区域图像分别执行503～507步骤，进行角点的提取、筛选和角点配对。

步骤六、对步骤五获取的n对角点配对结果，分别计算欧式距离，再根据卫星参数的几何模型计算建筑物的高度。

步骤601，分别计算配对角点对的欧式距离R_i，所述

其中i＝1,2...,n；

步骤602，根据步骤二的获取的卫星和太阳参数，构建建筑物高度和角点对距离之间的数学模型，计算转换系数K：

步骤603，分别计算建筑物上各角点的高度H_i，所述H_i＝R_i/K。

在实施例中，对两个典型建筑物角点及其阴影角点的距离通过步骤六求得平均为48.66米，根据步骤二获取太阳高度角38.2775°、太阳方位角161.8044°、卫星高度角86.2349°和卫星方位角280.5718°，从而计算出该建筑的高度为37.39米。该建筑位于通州区城市副中心，使用深达威(SNDWAY)测高仪SW-1500进行了实地测量的建筑物顶角高度为37.20米，验证了该方法的测量精度。同时对副中心附近的多栋建筑物进行高度测量，结果如下表，由表1可知：测量结果平均绝对误差0.61米，相对误差1.88％，说明本方法基于角点的建筑物测高方法具有较高的精度。

表1实验区建筑物高度测量结果

序号	地面测量高度	遥感监测高度	绝对误差m	相对误差％
					1	23.25	23.91	0.66	2.82％
2	27.20	27.39	0.19	0.70％
					3	58.54	58.68	0.14	0.25％
4	31.42	32.29	0.87	2.76％
					5	37.85	39.10	1.25	3.31％
6	36.66	37.19	0.53	1.45％
					7	23.25	23.91	0.66	2.82％
平均值			0.61	1.88％

本发明技术方案中除角点数据裁切外，其它数据处理操作可以采用python实现自动化运行处理。

以上所述是依据具体实施方式对本发明做的进一步详细说明，不能认定本发明具体实施仅局限于此。对于本发明所属及相关技术领域的技术人员，在基于本发明技术方法思路前提下，所作的拓展及操作方法、数据的替换，都应属于本发明包含范围。

Claims (4)

1.基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，具体步骤如下：

步骤一、选择亚米级高分辨率光学卫星对待监测建筑物进行成像；

步骤二、获取包含待监测区域的原始高分辨率光学卫星图像产品，根据卫星图像产品参数文件获得太阳高度角和方位角，以及卫星的高度角和方位角参数，从而确定待监测建筑物的阴影的成像方向；

步骤三、对原始高分辨率光学卫星图像进行辐射定标、大气校正、全色和多光谱数据融合处理、几何校正以及裁剪，得到高分辨率的待监测区域图像；

步骤四、对高分辨率的待监测区域图像的三波段数据进行增强处理，利用待监测建筑物的阴影在图像中的灰度值较低的特性进行多波段融合的阴影增强处理，得到待监测建筑物的增强阴影图像；

步骤五、根据步骤二确定的待监测建筑物的阴影成像方向，以及步骤四得到的待监测建筑物的增强阴影图像上分别裁切建筑物及其阴影在图像中对应的典型角点区域；对裁切区域的影像角特征进行提取，分别计算获得角点的精确坐标位置，分别得到待监测建筑物及其阴影的角点配对；

步骤六、对步骤五获取角点对配对结果，分别计算欧式距离，再根据卫星参数的几何模型计算建筑物的高度。

2.如权利要求1所述基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，具体步骤如下：步骤四中，首先分别计算由步骤三得到的高分辨率的待监测区域图像中RGB三个波段数据灰度值的均值和方差值；利用阴影的灰度值在各个多光谱波段中都远远低于正常地物灰度值的特性，计算三个波段组合的增强图像S_en，计算公式如下:

其中，B表示蓝色波段的灰度值，G表示绿色波段的灰度值，R表示蓝色波段的灰度值，μ_B表示蓝色波段的灰度值的均值，μ_G表示绿色波段的灰度值的均值，μ_R表示红色波段的灰度值的均值；σ_B表示蓝色波段的灰度值的方差，σ_G表示绿色波段的灰度值的方差，σ_R表示红色波段的灰度值的方差。

3.如权利要求1所述基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，具体步骤如下：步骤五的具体方法为：

步骤501，根据步骤二确定的待监测建筑物的阴影成像方向，在待监测建筑物的增强阴影图像上裁切待监测的建筑物楼顶的角点区域图像I_i(x,y)；其中，(x,y)为像素点坐标，i为选取裁切的角点区域图像数量；

步骤502，在待监测建筑物的增强阴影图像上裁切待监测建筑物楼顶的角点区域图像I_i(x,y)对应的阴影区域的角点图像I_i′(x,y)；

步骤503，计算图像I(x,y)在水平和垂直方向的一阶偏导函数I_x和I_y，利用高斯函数h(x,y,σ)作为窗口对协方差矩阵

进行平滑滤波，获得像素点的自相关矩阵M，所述

其中，(u,v)表示其它像素点相对于像素点(x,y)的偏移量；ω_u,v为高斯加权窗口函数；此处高斯函数表示为：

其中，σ为正态分布的标准偏差，确定函数的作用范围；

步骤504，计算角点响应函数R，R＝det(M)-k×trace²(M)，其中，det(M)为自相关函数矩阵M的行列式，trace(M)为自相关矩阵M的迹，k是可调参数；

步骤505，通过R的值大小来判断像素点是否为角点，对R进行归一化排序，设置阈值T进行判读是否为角点；其中，阈值T是可调参数，当R＞T时为角点，记录角点对应的角点像素坐标(x_i,y_i)_j，j为检测到的角点数量；

步骤506，对阴影区域的角点图像I_i′(x,y)重复步骤503、504和505，得到阴影区域的角点坐标(x_i',y_i')_j；

步骤507，对步骤505和步骤506得到的角点进行筛选配对，首先剔除不在阴影和建筑物边缘以及孤立的无效角点，再分别计算(x_i,y_i)_j和(x_i',y_i')_j中剩余角点对连线的角度，根据卫星空间几何模型获取的方位角ω_i，筛选角度一致的为配对角点对(x_i,y_i)_n和(x_i',y_i')_n，其中n是筛选配对后得到的角点对数量，

所述

步骤508，对步骤501和502裁切的i个角点区域图像分别执行503～507步骤，进行角点的提取、筛选和角点配对。

4.如权利要求1所述基于角点的高分辨率光学遥感卫星图像建筑物高度监测方法，具体步骤如下：步骤六的具体方法为：

步骤601，分别计算配对角点对的欧式距离R_i，所述

其中i＝1,2...,n；

步骤602，根据步骤二的获取的卫星和太阳参数，构建建筑物高度和角点对距离之间的数学模型，计算转换系数K：

步骤603，分别计算建筑物上各角点的高度H_i，所述H_i＝R_i/K。

Images