CN101894382A - 一种集成LiDAR点云的卫星立体影像阴影计算方法 - Google Patents

Abstract

本发明一种集成LiDAR点云的卫星立体影像阴影计算方法，首先对原始LiDAR点云数据进行滤波，得到包含建筑物和树木等非地面点的LiDAR点云，然后使用带距离控制的点包裹算法得到地物点集的外包多边形，通过判断点与多边形的位置关系得到地物轮廓多边形，对地物轮廓多边形进行滤波和形状简化得到建筑物外部轮廓多边形；由影像拍摄瞬间太阳高度角和方位角信息可以求解建筑物阴影的物方坐标，根据卫星立体影像像方RFM几何校正结果确定物方坐标与像方坐标之间的转换关系，从而得到建筑物阴影在影像像方空间的坐标和位置。本发明可以高效、快速地获取高分辨率卫星影像的阴影区域与坐标范围，具有稳健性，得到的阴影位置坐标精确。

Description

一种集成LiDAR点云的卫星立体影像阴影计算方法

技术领域

本发明属于LiDAR点云和卫星立体遥感数据城市空间数据采集与更新应用领域，涉及一种集成LiDAR点云的卫星立体影像阴影计算方法，特别是基于机载LiDAR点云和高分辨率卫星立体影像的城市建筑物阴影计算方法。

背景技术

20世纪以来地球科学进步的一个突出标志是对地观测技术的飞速发展，主要包括机载对地观测技术和星载对地观测技术。近年来的机载和星载对地观测技术分别以机载激光扫描(LiDAR，Light Detection And Ranging)的出现和高分辨率卫星影像的广泛应用为主要代表，于是对应的是人们从遥感数据中获取的信息越来越丰富，速度越来越快，自动化程度也越来越高。在高分辨率卫星影像的各种研究中，核心的问题是影像的目标识别与信息解译。在城市地区，由于建筑物特别是高层建筑物的存在，进行影像分析面临的突出问题是建筑物存在阴影，阴影会遮挡到其他建筑物、道路等地物，给影像地物信息识别解译造成很大的困难。

要减弱或消除影像上建筑物阴影造成的影响，首先需要确定阴影在影像上的位置，再对其进行相应的处理。对于遥感影像上阴影区域的检测，现在一般有两种方式，一种是基于影像颜色空间的信息，这种方式不能很好地区分与阴影颜色相近的地物如水体等；另外一种是利用影像的几何信息，如阴影边缘检测等，这种方式很多时候也会造成不同地物的误判。

发明内容

为了解决现有依靠单一卫星影像数据进行阴影检测与计算方法的不足，本发明目的在于提供一种集成LiDAR点云的卫星立体影像阴影计算方法，结合机载LiDAR点云数据和高分辨率卫星立体影像，可以稳健地得到高分辨率卫星立体影像上的建筑物阴影位置和精确坐标。

为达到以上目的，本发明所采用的解决方案是：

一种集成LiDAR点云的卫星立体影像阴影计算方法，包括以下步骤：

(1)计算机载LiDAR点云数据建筑物外部轮廓多边形；

首先对原始的机载LiDAR点云数据进行滤波，去除地面点和不属于任何地物的噪声点；然后使用带距离控制的卷包裹算法得到地物点集的外包多边形，判断点与多边形的位置关系，得到地物轮廓多边形；进行多边形滤波，得到建筑物多边形，再进行建筑物轮廓多边形形状简化，得到最终的建筑物外部轮廓多边形。

(2)确定建筑物阴影在物方空间坐标位置；

根据卫星立体影像(同轨或异轨)在拍摄瞬间的太阳高度角信息和太阳方位角信息，列立建筑物棱线高度和对应的阴影的方向和长度的关系式，求解建筑物多边形中每个点在地面阴影的坐标，得到建筑物阴影在物方空间的坐标。

(3)卫星立体影像几何定位校正，确定物方坐标与像方坐标之间的转换关系；

根据卫星影像提供的有理多项式系数RPC列立卫星立体影像的有理函数模型RFM，并通过加入控制点的方式进行像方空间RFM几何校正，选择像方RFM几何校正模型(如平移模型、平移缩放模型、仿射变换模型、二阶多项式模型)，在影像区域内加入6-8个平面分布均匀的控制点，提高几何定位的精度，建立卫星立体影像物方坐标和像方坐标之间的精确转换关系。

(4)确定建筑物阴影在影像像方空间的坐标和位置；

根据步骤(2)所得到的建筑物阴影在物方空间的坐标，以及步骤(3)所得到的卫星立体影像物方坐标和像方坐标之间的精确转换关系，分别计算建筑物阴影区域每个点在影像像方空间的坐标，进而得到建筑物阴影在影像像方的整体坐标以及位置。

所述去除地面点和噪声点，具体方法为：去除地面点是以一定的高程值为阈值，去除小于该阈值的LiDAR点，保留大于该阈值的LiDAR点；去除噪声点是计算LiDAR点云中某点与其8个最邻近点的平均距离，该距离大于一定阈值则去除。

所述带距离控制的卷包裹算法，具体方法为：a).在去除地面点和噪声点的点集中找到X坐标最小及最大、Y坐标最小及最大的点，取其中之一作为出发点，并记该顶点为P1；b).从点P1出发，做一条垂直线段L，长度为D，D为距离控制阈值；c).以P1为固定点，逆时钟方向旋转L，直到接触到点集中的某点P2，则P2必是地物外部轮廓多边形上的另一个顶点；d).以P2为固定点继续旋转L，直到接触到点集中的某点P3，则P3也必是地物外部轮廓多边形中的一个顶点；e).重复步骤d)直到L又重新回到P1为止。

当以某点为固定点旋转L，接触到的点有多个时，则将这些点按照与固定点之间的距离由小到大的排列成一组，则这一组点都是地物外部轮廓多边形的顶点，然后将与固定点之间距离最大的点作为新的固定点，继续旋转L。

得到地物点集的外包多边形后，去除外包多边形的顶点及外包多边形内部的点，只保留外包多边形外部的数据点，对保留的数据点再次应用带距离控制卷包裹算法，得到另一个地物的外部轮廓多边形，重复前述方法直至提取所有的地物轮廓多边形。

采用面积阈值进行多边形滤波；采用最小二乘原理对多边形形状简化。

由于采用了上述方案，本发明具有以下特点：本发明可以集成机载LiDAR点云，高效、快速地获取高分辨率卫星影像的阴影区域与坐标范围，不受卫星影像上地物颜色与波谱等特征的限制，具有稳健性，得到的阴影位置坐标精确，可以促进城市地区卫星影像目标识别与信息解译的研究。

附图说明

图1为本发明方法的流程图。

图2为建筑物阴影与太阳高度角和太阳方位角之间的关系。

具体实施方式

以下结合附图所示实施例对本发明作进一步的说明。

参见附图1，本发明集成机载LiDAR点云数据与卫星立体影像进行影像阴影特别是建筑物阴影的检测与坐标计算，首先对原始LiDAR点云数据进行滤波，得到包含建筑物和树木等非地面点的LiDAR点云，然后使用带距离控制的点包裹算法得到地物点集的外包多边形，通过判断点与多边形的位置关系得到地物轮廓多边形，对地物轮廓多边形进行滤波和形状简化得到建筑物外部轮廓多边形；由影像拍摄瞬间太阳高度角和方位角信息可以求解建筑物阴影的物方坐标，根据卫星立体影像像方RFM几何校正结果确定物方坐标与像方坐标之间的转换关系，从而得到建筑物阴影在影像像方空间的坐标和位置。

为了更好地理解本发明的技术方案，以下结合附图和实施例作进一步的详细描述。

(1)计算机载LiDAR点云数据建筑物外部轮廓多边形；

对于原始的机载LiDAR点云，由于包含所有的地面和地物扫描点，首先要进行点云数据的滤波，除去点云中的地面点和噪声点，去除地面点一般以一定的高程值为阈值，小于此值的LiDAR点去除，大于阈值的点保留；对于噪声点，不属于任何地物对象，可以计算LiDAR点云中某点与其k个最邻近点的平均距离，若此距离大于一定阈值，则去除，实施例中取k＝8，阈值为4m，从92729个点的实验数据中去除噪声点149个。

将原始的LiDAR点云除去地面点和噪声点后，将其投影到某一高程为H的平面上，则这些投影点在平面上表现为一系列离散的点族，其中每一个点族代表了某一地物，用带距离控制的卷包裹算法分别求取这些点族的外包，得到的就是这些点族所表示的地物在高程为H处的外部轮廓多边形。

带距离控制的卷包裹算法，其具体过程如下：

a).在点集中找到X坐标最小及最大、Y坐标最小及最大的点，这四个点肯定在多边形上，取其中之一作为出发点，并记该顶点为P1；b).从点P1出发，做一条垂直线段L，长度为D，D为距离控制阈值；c).以P1为固定点，逆时钟方向旋转L，直到接触到点集中的某点P2，则P2必是地物外部轮廓多边形上的另一个顶点。以P2为固定点继续旋转L，直到接触到点集中的某点P3，则P3也必是地物外部轮廓多边形中的一个顶点……一直这样下去，直到L又重新回到P1为止。需要说明的是，当以某点为固定点旋转L，接触到的点有多个时，则将这些点按照与固定点之间的距离由小到大的排列成一组，则这一组点都是地物外部轮廓多边形的顶点。然后将与固定点之间距离最大的点作为新的固定点，继续旋转L……

使用卷包裹算法时，距离控制是十分重要的，距离太小会将一个建筑物分割得支离破碎，距离太大则会忽略建筑物的一些细节，或者使得相互靠近的两栋建筑或者是一栋建筑和周围的树木融合到一起，实施例中使用D为3-4m的距离可以得到比较好的效果。

投影到高程为H的水平面上的点族运用带距离控制的卷包裹算法求得一个地物外包轮廓多边形后，必须把外包多边形的所有顶点和位于外包多边形内部的所有点都去除掉，只留下位于外包多边形外部的数据点，对这些数据点再次应用带距离控制的卷包裹算法，获取另一个地物的外部轮廓多边形，依此不停地循环，直到把所有的地物轮廓多边形都提取出来。

得到地物轮廓多边形中，包含建筑、树木，电线杆等，采用面积阈值将表示树木、电线杆等其他地物的多边形去除掉，只留下代表建筑物的多边形，实施例中采用面积阈值为8m²可以得到比较好的过滤效果。

在得到建筑物轮廓多边形以后，其顶点数目一般比较多，需要对多边形进行简化，采用最小二乘原理对建筑物多边形的边线进行简化，减少顶点的个数，以方便后续的计算过程。

(2)确定建筑物阴影在物方空间坐标位置；

得到某一高程H的外部轮廓多边形以后，通过改变H的值，可以获得建筑物在不同高度的外部轮廓多边形。对于一栋高层建筑物，从某一较小高程平面H＝h开始，以较小的高程差(1m)，不断获取建筑的外部轮廓多边形，直到h＞＝建筑物高程为止，可以得到建筑物体在空间的外部轮廓多边形变化情况。

根据卫星立体影像在拍摄瞬间的太阳高度角和太阳方位角信息，列立建筑物棱线高度和对应的阴影的方向和长度的关系式。如附图2，假设直线AO为建筑物一条竖直的棱线，点A的坐标为(X₀，Y₀，H)，点A′为点A在太阳光下的阴影，点O为点A在地平面上的竖直投影，其坐标为(X₀，Y₀，h)，太阳高度角为α，太阳方位角为β，则AO在太阳光下的阴影长度为：

$s=\frac{(H-h)}{\tan \mathrm{\α}}---\left(1\right)$

点A的阴影A′的坐标为：

$\left\{\begin{array}{c}{X}_A=X_0+\mathrm{\ΔX}=X_0+s\cos (180-\mathrm{\β})=X_0-\frac{(H-h)}{\tan \mathrm{\α}}\cos \mathrm{\β}\\ Y_A=Y_0+\mathrm{\ΔY}=Y_0-\sin (180-\mathrm{\β})=Y_0-\frac{(H-h)}{\tan \mathrm{\α}}\sin \mathrm{\β}\\ Z_A=h\end{array}\right.---\left(2\right)$

在得到建筑物多边形中每个点在地面上阴影的物方坐标以后，即得到整个建筑物的阴影在物方空间的坐标。

(3)确定卫星影像物方坐标与像方坐标之间的转换关系；

卫星影像物方坐标与像方坐标之间的转换关系通过对卫星立体影像进行基于像方有理函数模型RFM的几何校正实现。有理函数模型RFM是目前高分辨率卫星影像几何定位的通用传感器模型，由于其形式简单、定位精度高等优点而被广泛采用，RFM模型的定义式如下：

$\begin{array}{c}{r}_n=\frac{P_1(X_n,Y_n,Z_n)}{P_2(X_n,Y_n,Z_n)}\\ c_n=\frac{P_3(X_n,Y_n,Z_n)}{P_4(X_n,Y_n,Z_n)}\end{array}---\left(3\right)$

其中r_n，c_n是分别是标准化的影像行列坐标，单位为像素；X_n，Y_n，Z_n为地面坐标系像方坐标的标准化坐标。坐标的标准化过程使用以下的等式计算：

$\begin{array}{ccccc}{r}_n=\frac{r-r_0}{r_s}& c_n=\frac{c-c_0}{c_s}& X_n=\frac{X-X_0}{X_s}& Y_n=\frac{{Y-Y}_0}{Y_s}& Z_n=\frac{Z-Z_0}{Z_s}\end{array}---\left(4\right)$

其中r₀，c₀是两个影像坐标的平移参数，r_s，c_s为缩放参数。同样，X₀，Y₀，Z₀为三个地面坐标的平移参数，X_s，Y_s，Z_s为缩放参数。

多项式中每一项地面坐标的最大幂通常不超过3，每一项各个坐标分量的幂的总和也不超过3。P_i(i＝1，2，3，4)是物方地面坐标(X_n，Y_n，Z_n)的三阶多项式函数，称为有理多项式系数RPC。

由于RFM中的RPC系数一般仅仅采用卫星轨道参数计算获得，一般只能获取20m左右的精度，针对此问题，通常利用控制点对原始的RFM模型进行几何校正，常用的有物方空间校正和像方空间校正。像方RFM几何校正模型可以在公式(3)的基础上加入变换参数得到：

$\begin{array}{c}\frac{(r+A_0+A_{1}r+A_{2}c+A_3\mathrm{rc}+A_4{r}^2+A_5{c}^2)-r_0}{r_s}=\frac{P_1(X_n,Y_n,Z_n)}{P_2(X_n,Y_n,Z_n)}\\ \frac{(c+B_0+B_{1}c+B_{2}r+B_3\mathrm{rc}+B_4{r}^2+B_5{c}^2)-c_0}{c_s}=\frac{P_3(X_n,Y_n,Z_n)}{P_4(X_n,Y_n,Z_n)}\end{array}---\left(5\right)$

式(5)中的变换参数有以下四种选择：

1)A₀，B₀，表示像方坐标平移，代表平移模型；

2)A₀，A₁，B₀，B₁，表示像方坐标平移和缩放，为平移缩放模型；

3)A₀，A₁，A₂，B₀，B₁，B₂，表示像方仿射变换模型；

4)A₀，A₁...A₅，B₀，B₁...B₅，表示像方二阶多项式模型。

实施例中的卫星立体影像为QuickBird高分辨率卫星立体影像，通过加入GPS测量点作为控制点，采用最小二乘平差求解像方RFM几何校正模型，可以得到卫星影像的物方坐标与像方坐标之间的转换关系，通过分析控制点个数、分布以及不同的像方RFM几何校正模型，得到最后的转换关系。

(4)确定建筑物阴影在影像像方空间的坐标和位置；

根据机载LiDAR数据得到的建筑物阴影在物方空间的坐标，以及卫星立体影像物方坐标和像方坐标之间的精确转换关系，分别计算建筑物阴影区域每个点在影像像方空间的坐标，进而得到建筑物阴影在影像像方的整体坐标以及位置。

上述的对实施例的描述是为便于该技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对这些实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于这里的实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (10)

1.一种集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：其具体包括以下步骤：

1)计算机载LiDAR点云数据建筑物外部轮廓多边形；

2)确定建筑物阴影在物方空间坐标位置；

3)卫星立体影像几何定位校正，确定物方坐标与像方坐标之间的转换关系；

4)确定建筑物阴影在影像像方空间的坐标和位置。

2.根据权利要求1所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述计算机载LiDAR点云数据建筑物外部轮廓多边形，具体方法为：首先对原始的机载LiDAR点云数据进行滤波，去除地面点和噪声点；然后使用带距离控制的卷包裹算法得到地物点集的外包多边形，判断点与多边形的位置关系，得到地物轮廓多边形；进行多边形滤波，得到建筑物多边形，再进行建筑物轮廓多边形形状简化，得到最终的建筑物外部轮廓多边形。

3.根据权利要求1所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述确定建筑物阴影在物方空间坐标位置，具体方法为：根据卫星立体影像同轨或异轨在拍摄瞬间的太阳高度角信息和太阳方位角信息，列立建筑物棱线高度和对应的阴影的方向和长度的关系式，求解建筑物多边形中每个点在地面阴影的坐标，得到建筑物阴影在物方空间的坐标。

4.根据权利要求1所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述确定物方坐标与像方坐标之间的转换关系，具体方法为：根据卫星影像提供的有理多项式系数RPC列立卫星立体影像的有理函数模型RFM，并通过加入控制点的方式进行像方空间RFM几何校正，选择像方RFM几何校正模型，在影像区域内加入6-8个平面分布均匀的控制点，建立卫星立体影像物方坐标和像方坐标之间的精确转换关系。

5.根据权利要求1所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述确定建筑物阴影在影像像方空间的坐标和位置，具体方法为：根据建筑物阴影在物方空间的坐标，以及卫星立体影像物方坐标和像方坐标之间的精确转换关系，分别计算建筑物阴影区域每个点在影像像方空间的坐标，进而得到建筑物阴影在影像像方的整体坐标以及位置。

6.根据权利要求2所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述去除地面点和噪声点，具体方法为：去除地面点是以一定的高程值为阈值，去除小于该阈值的LiDAR点，保留大于该阈值的LiDAR点；去除噪声点是计算LiDAR点云中某点与其8个最邻近点的平均距离，该距离大于一定阈值则去除。

7.根据权利要求2所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：所述带距离控制的卷包裹算法，具体方法为：a).在去除地面点和噪声点的点集中找到X坐标最小及最大、Y坐标最小及最大的点，取其中之一作为出发点，并记该顶点为P1；b).从点P1出发，做一条垂直线段L，长度为D，D为距离控制阈值；c).以P1为固定点，逆时钟方向旋转L，直到接触到点集中的某点P2，则P2必是地物外部轮廓多边形上的另一个顶点；d).以P2为固定点继续旋转L，直到接触到点集中的某点P3，则P3也必是地物外部轮廓多边形中的一个顶点；e).重复步骤d)直到L又重新回到P1为止。

8.根据权利要求7所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：当以某点为固定点旋转L，接触到的点有多个时，则将这些点按照与固定点之间的距离由小到大的排列成一组，则这一组点都是地物外部轮廓多边形的顶点，然后将与固定点之间距离最大的点作为新的固定点，继续旋转L。

9.根据权利要求2所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：得到地物点集的外包多边形后，去除外包多边形的顶点及外包多边形内部的点，只保留外包多边形外部的数据点，对保留的数据点再次应用带距离控制卷包裹算法，得到另一个地物的外部轮廓多边形，重复前述方法直至提取所有的地物轮廓多边形。

10.根据权利要求2所述的集成LiDAR点云的卫星立体影像阴影计算方法，其特征在于：采用面积阈值进行多边形滤波；采用最小二乘原理对多边形形状简化。

Images