CN103729853B - 三维gis辅助下的高分遥感影像建筑物倒损检测方法 - Google Patents

Abstract

一种三维GIS辅助下的高分遥感影像建筑物倒损检测方法，包括利用灾前三维GIS或影像RPC参数对灾后遥感影像进行几何纠正；获取灾前三维GIS中各建筑物对象的相应二维矢量数据及高度数据，构成灾前建筑物集合；将灾前二维矢量数据叠加至配准后的灾后遥感影像中，为各建筑物对象构建局部缓冲区，并在缓冲内进行水平集演化分割获得建筑物分割对象；提取灾后建筑物集合多个特征并检测隶属度，结合各个特征证据概率计算建筑物发生倒损的置信度；对灾后遥感影像中未处理区域进行水平集分割，滤掉面积小于预设面积阈值的对象，利用已处理区域训练分类。本发明采用三维GIS、水平集演化、证据理论，从多种特征入手解决建筑物倒损检测问题。

Description

三维GIS辅助下的高分遥感影像建筑物倒损检测方法

技术领域

本发明涉及遥感影像应用技术领域，特别涉及一种三维GIS辅助下的高分遥感影像建筑物倒损检测方法。

背景技术

自然灾害的频繁发生，给国家人民的财产带来极大损失，建筑物作为城市生产生活的重要地物，灾后对其倒损情况快速做出精确评估具有重要意义。遥感具有重访周期短，探测范围大、数据综合性高、经济效益高等特点，是灾害评估的核心手段之一。目前国家减灾部门对建筑物的倒损灾害评估致力于统计模型的建立、定性分析与评估，缺乏定量的特别是精细的定量评估，高空间分辨率(通常空间分辨率小于2米)遥感影像具有分辨率高、地物细节丰富等特点，随着其获取越来越便捷，逐渐成为定量化灾害精细评估的主要手段。

通常自然灾害发生后建筑物的倒损形式包括结构整体倒塌、整体沉降、部分倒塌、下部倒塌、中间层倒塌、墙体裂而不倒等。显而易见，对于建筑物整体沉降，下部倒塌等情况，倒损检测不仅需要二维信息的提取与分析，建筑物高度信息的检测也至关重要。目前利用遥感进行建筑物倒损检测的典型方法包括以下几种：1)利用单张灾后影像的建筑物倒损检测方法，通过对遥感影像解译分类，提取建筑物倒损信息，但该方法难以提取建筑物高度变化信息，即使能利用阴影等方式获得部分建筑物高度，但不具有普适性。2)利用两/多时相遥感影像的建筑物倒损检测方法，包括直接比较法、分析后比较法及统一模型法等。直接比较法即对不同数据源进行直接比较，对象主要包括像素，纹理特征，边缘特征以及各种复杂的变换后特征，如植被指数、主成分变换、独立成分变换、典型相关变化等；分析后比较法即对不同数据源信息提取后比较，对象包括类别、目标对象等；统一模型法即将不同数据源纳入统一的模型进行变化检测，将变化检测的方法和过程作为一个整体，采用统一的平差模型进行迭代求解。以上所述利用高分遥感影像的变化检测方法都取得了不错的效果，但是由于这些方法主要是基于二维数据的变化检测，难以检测建筑物高度变化信息，对检测那些部分倒塌、整体沉降、下部倒塌或中间层倒塌等高度发生变化的建筑物具有先天缺陷。3)利用不同时相LIDAR(LightDetectionandRanging，激光雷达)数据或立体像对的建筑物倒损检测，LIDAR数据及立体像对都包含地物的三维信息，通过DSM(DigitalSurfaceModel，数字表面模型)的提取与比较分析可以很好检测建筑物高度的变化，但是基于LIDAR数据的获取方式及发展状况，灾害发生地区通常很难具备多时相LIDAR数据，立体像对也存在类似问题，即使单时相立体像对可以检测建筑物倒损情况，但是存在立体像对幅面较小，而且需要专业的摄影测量处理软件，获取DSM及建筑物三维需要复杂的工作等问题。因此需要迫切寻找一种数据易获取、检测效率高、检测结果相对精确且顾及到建筑物高度信息的变化检测方法。

发明内容

本发明提出了一种三维GIS(GeographicInformationSystem，地理信息系统)辅助下的高分遥感影像建筑物倒损检测方法，该方法将建筑物高度信息变化情况进行了检测，加强了建筑物变化检测的准确性。

本发明的技术方案为一种三维GIS辅助下的高分遥感影像建筑物倒损检测方法，包括以下步骤：

步骤1，利用灾前三维GIS或影像RPC参数对灾后遥感影像进行几何纠正；

步骤2，获取灾前三维GIS中各建筑物对象的三维GIS数据，包括二维矢量数据及高度数据，构成灾前建筑物集合V_B；

步骤3，记灾前建筑物集合V_B中每个建筑物对象为v_i，将相应二维矢量数据叠加至步骤1几何纠正后的灾后遥感影像中，在灾后遥感影像上为建筑物对象v_i构建局部缓冲区，将v_i所对应的遥感影像区域的轮廓作为初始水平集，在局部缓冲区内部的影像中进行水平集分割，得到建筑物分割对象v_iSeg；所有建筑物分割对象构成灾后建筑物集合V_Seg；

步骤4，针对灾后遥感影像中步骤3所得灾后建筑物集合V_Seg相应部分提取n个特征，并检测各特征的隶属度ρ_j，结合各个特征的预设证据概率R_j计算建筑物发生倒损的置信度，当建筑物对象v_iSeg的置信度P_Judge处于置信区间时，认为建筑物发生倒损，并入集合B_s；P_Judge处于拒绝区间时，得到未发生倒损建筑物对象，并入集合B_n；P_Judge处于不确定区间时，建筑物对象疑似发生倒损，并入集合B_d；

所述计算建筑物发生倒损的置信度采用以下公式，

$P_{Judge}=\frac{\underset{j=1}{\overset{n}{\Σ}}{R}_j\×{\mathrm{\ρ}}_j}{n}$

步骤5，将灾后遥感影像中步骤3所得灾后建筑物集合V_Seg的相应部分作为已处理区域，其他部分作为未处理区域；以集合B_n的区域平均灰度值作为阈值确定初始水平集，对灾后遥感影像中未处理区域进行水平集分割；滤掉分割对象中面积小于预设面积阈值的对象，剩余的对象构成集合B_new；

步骤6，分别从步骤4所得集合B_s、集合B_n形成的区域中随机提取像元光谱特征值作为训练样本训练SVM分类器，利用训练出的SVM分类器对步骤5所得集合B_new中的对象分别进行验证，将验证结果为倒损建筑物类别的对象并入集合B_s，构成新的倒损建筑物集合B′_s；将验证为未倒损建筑物类别的对象并入集合B_n，构成新的未倒损建筑物集合B′_n。

而且，步骤1中，判断输入的灾后遥感影像是否具有RPC参数，有则采用RPC参数对灾后遥感影像进行几何纠正，否则采用三维GIS对灾后遥感影像进行几何纠正。

而且，步骤4中，提取灾后遥感影像中步骤3所得灾后建筑物集合V_Seg相应部分的面积、纹理和高度特征，并检测各特征的隶属度，包括对每个建筑物对象v_i及灾后遥感影像上的相应建筑物分割对象v_iSeg分别进行以下处理，

a)对于某建筑物对象v_i的面积S_iFore、灾后遥感影像上相应的建筑物分割对象v_iSeg的面积为S_iAfter，面积特征提取结果为面积差ΔS_i＝S_iFore-S_iAfter，检测面积特征的隶属度如下，

$R_{i\mathrm{\&Delta;}S}\left(Area\right)=\left\{\begin{array}{cc}0& {\mathrm{\&Delta;S}}_i<Min\mathrm{\&Delta;}S\\ \frac{{\mathrm{\&Delta;S}}_i-Min\mathrm{\&Delta;}S}{Max\mathrm{\&Delta;}S-Min\mathrm{\&Delta;}S}& Min\mathrm{\&Delta;}S\&le;{\mathrm{\&Delta;S}}_i\&le;Max\mathrm{\&Delta;}S\\ 1& {\mathrm{\&Delta;S}}_i>Max\mathrm{\&Delta;}S\end{array}\right.$

其中，MinΔS表示面积差的预设最小阈值，MaxΔS表示面积差的预设最大阈值；

b)对于某建筑物对象v_i，计算灾后遥感影像上相应的建筑物分割对象v_iSeg的灰度共生矩阵中角二阶矩特征值ASM_i作为纹理特征提取结果，检测纹理特征的隶属度如下，

$R_{iASM}\left(ASM\right)=\left\{\begin{array}{cc}1& {\mathrm{ASM}}_i<MinASM\\ \frac{MaxASM-{\mathrm{ASM}}_i}{MaxASM-MinASM}& MinASM\&le;{\mathrm{ASM}}_i\&le;MaxASM\\ 0& {\mathrm{ASM}}_i>MaxASM\end{array}\right.$

其中，MinASM表示预设最小角二阶矩特征值阈值，MaxASM表示预设最大角二阶矩特征值阈值；

c)对于某建筑物对象v_i，建筑物对象v_i灾前的高度值为h_iFore，建筑物对象v_i在灾后遥感影像中的高度值h_iAfter，计算灾前灾后的高度差Δh_i＝h_iFore-h_iAfter作为高度特征提取结果，检测高度特征的隶属度如下，

$R_{i\mathrm{\&Delta;}H}\left(Height\right)=\left\{\begin{array}{cc}0& {\mathrm{\&Delta;H}}_i<Min\mathrm{\&Delta;}H\\ \frac{{\mathrm{\&Delta;H}}_i-Min\mathrm{\&Delta;}H}{Max\mathrm{\&Delta;}H-Min\mathrm{\&Delta;}H}& Min\mathrm{\&Delta;}H\&le;{\mathrm{\&Delta;H}}_i\&le;Max\mathrm{\&Delta;}H\\ 1& {\mathrm{\&Delta;H}}_i>Max\mathrm{\&Delta;}H\end{array}\right.$

其中，MinΔH表示高度差的预设最小值阈值，MaxΔH表示高度差的预设最大值阈值。

而且建筑物对象v_i在灾后遥感影像中的高度值h_iAfter按如下方式提取，

1)由太阳方位角α_sun确定阴影相对于建筑物对象v_i的成像方向，设斜率为k，求出相应建筑物分割对象v_iSeg的最小外接矩形Rect_i，其中一组边平行于阴影成像方向，斜率为k，记为a_i1、a_i2，另一组边垂直于阴影成像方向记为b_i1、b_i2；根据建筑物对象v_i灾前的高度值h_iFore和影像成像时太阳与卫星的角度参数，求出由建筑物对象v_i投影出的阴影长度l_i；将最小外接矩形Rect_i的边a_i1、a_i2沿建筑物成像方向延伸l_i，形成新的矩形Rect′_i，并以预设缓冲区距离d_Buffer2为其构建局部缓冲区Rect′_iBuf，得到待处理影像阴影区域Ω_ishadow；

2)在待处理影像阴影区域Ω_ishadow内进行OSTU分割，然后对所得分割对象sdw_i矢量化后的轮廓曲线进行压缩处理，利用斜率为k的直线以一定的间距对压缩处理后所得多边形进行分块，分别计算每条直线与多边形的两交点间的线段长度，求较长的多条线段的长度平均值作为阴影有效长度sdw_iLen；

3)根据阴影有效长度sdw_iLen和影像成像时太阳与卫星的角度参数，计算得到建筑物对象v_i在灾后遥感影像中的高度值h_iAfter。

经过十多年数字城市的建设，很多城市建立了地理信息系统而且具有该地区的三维模型数据。三维GIS数据通常具有精确的地理坐标和高度信息，可作为建筑物倒损检测中的地物定位及标准参考，而高分遥感影像具有易获取、纹理信息丰富、空间结构清晰等特点，两种数据源融合处理能够优势互补，为倒损建筑物的快速检测提供有力保障。对比现有技术，本发明的特点如下：

1)提出了一种利用多特征融合的方式对建筑物进行变化检测，并运用证据理论对特征提取结果进行分析处理。

2)变化检测由二维扩展到三维空间中，在检测二维平面建筑物倒损情况的基础上，提出了通过提取灾后遥感影像中建筑物阴影来获得灾后建筑物高度，并将其与三维GIS的建筑物高度作比较来检测建筑物对象在三维上的变化。

附图说明

图1是本发明实施例的流程图。

具体实施方式

本发明所提供三维GIS辅助下的高分遥感影像建筑物倒损检测方法是在利用灾前三维GIS建筑物数据对遥感影像几何纠正之后，在三维GIS数据的辅助下提取遥感影像中的建筑物对象，并对不同时相建筑物对象的面积、高度数据作比较，对遥感影像中纹理特征进行提取与检测，利用证据理论对各特征进行组合计算建筑物倒损的置信度，以此判断建筑物倒损的情况。对于新增建筑物利用训练样本的方式来进行验证。本发明技术方案可基于计算机软件技术实现自动运行流程。

以下结合附图和实施例详细说明本发明技术方案，流程图如图1所示，实施例的技术方案流程包括以下步骤：

步骤1，利用灾前三维GIS或影像RPC参数对灾后遥感影像进行几何纠正：对灾后遥感影像进行纠正，可根据影像是否具有RPC参数，分别采用RPC模型参数和三维GIS对灾后遥感影像进行几何纠正，获得具有真实地理坐标的遥感影像。本发明所述灾前三维GIS包括灾前三维GIS建筑物数据。

实施例中，输入的灾后遥感影像为高分辨率光学遥感影像，先判断输入的灾后遥感影像是否具有RPC参数，有则采用RPC参数对灾后遥感影像进行几何纠正。否则，采用三维GIS对灾后遥感影像进行几何纠正，操作时可以分别在三维GIS和灾后遥感影像中选取多对控制点，采用多项式纠正的方式对遥感影像进行几何纠正。利用影像RPC参数求解影像平面坐标为现有技术，本发明不予赘述。

步骤2，获取灾前三维GIS中建筑物对象的三维GIS数据，包括二维矢量数据、建筑物高度数据，并可为每个建筑物对象编号，构成灾前建筑物集合V_B。设从灾前三维GIS中获取的建筑物对象集合为V_B＝{v₁,v₂,v₃,…v_m}，其中v_i为集合中的一个建筑物对象，i＝1,2，…m，m为建筑物对象总数。

实施例中，获取灾前三维GIS中每个建筑物对象的三维数据，包括二维矢量数据、建筑物高度数据，所得二维矢量数据即建筑物对象的相应灾前建筑物矢量数据，并为其建立索引。

步骤3，将要处理的灾前建筑物矢量数据叠加至步骤1几何纠正后的灾后遥感影像中，遍历访问建筑物矢量数据，为待处理建筑物对象建立局部缓冲区，并在缓冲区轮廓以内的影像区域进行分割。

具体实施时，可依次对每一个灾前建筑物矢量数据分别进行处理：将要处理的灾前建筑物矢量数据叠加至步骤1配准后的灾后遥感影像上，依据灾前建筑物集合V_B中建筑物对象与灾后遥感影像的空间关系，在灾后遥感影像上为当前待处理的建筑物对象v_i构建局部缓冲区，将v_i所对应的遥感影像区域的轮廓作为初始水平集，在缓冲区内部的灾后遥感影像中进行局部水平集演化，得到灾后遥感影像的建筑物分割对象v_iSeg，并入灾后建筑物分割集合V_Seg。

步骤3中，根据编号依次对V_B中对象进行处理，每次处理一个建筑物对象的流程可设计为包括如下子步骤，

a)获取待处理的建筑物对象v_i的地理坐标，在灾后遥感影像上相应位置处为v_i构建缓冲区v_iBuf，得到当前待处理影像区域；

b)以当前待处理建筑物对象v_i所对应的遥感影像区域的轮廓作为初始水平集，在待处理影像区域v_iBuf内进行水平集演化分割得到建筑物区域，作为灾后影像的建筑物分割对像v_iSeg。

实施例具体实现方式为，对每一个待处理的建筑物对象v_i，将要处理的灾前建筑物矢量数据叠加到步骤1几何校正后的灾后遥感影像后，执行以下操作：

首先，以预设缓冲区距离d_Buffer1在灾后遥感影像上相应位置处为v_i建立局部缓冲区v_iBuf，缓冲区v_iBuf轮廓以内的影像区域为当前待处理区域Ω。d_Buffer1可由本领域技术人员根据情况自行预先设定。

然后，以当前待处理的建筑物对象v_i叠加到灾后遥感影像后的轮廓边界为初始零水平集，可表示为φ(x,y,t)＝0，下文简记为φ，其中(x,y)为像素点坐标，t为时间参数。在a)中所得当前待处理区域Ω内，进行局部水平集演化分割。曲线φ将待处理区域分为两部分，曲线φ内部φ>0，外部φ<0,这两个区域的平均灰度值分别为u₁、u₂。现有技术基于水平集C-V模型的能量函数为：

$\begin{array}{c}E(\mathrm{\&phi;},c_1,c_2)={\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{\left(u_0\right(x,y)-u_1)}^{2}H\left(\mathrm{\&phi;}\right)dxdy+\\ {\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{\left(u_0\right(x,y)-u_2)}^2(1-H(\mathrm{\&phi;}\left)\right)dxdy+\\ \mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}\right)dxdy+\mathrm{\&nu;}{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\left|\&dtri;\mathrm{\&phi;}\right|dxdy\end{array}---\left(1\right)$

式中，H(φ)是Heaviside函数，其形式为：

δ(φ)为Dirac函数，其形式为u₀(x,y)是待处理影像区域某一点灰度值，为当前点梯度的模，系数λ₁，λ₂＞0，μ，ν≥0为固定参数，一般建议取λ₁＝λ₂＝1，μ＝0，ν＝1，(1)对应的偏微分方程：

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\&CenterDot;\&lsqb;\mathrm{\&mu;}\&dtri;\&CenterDot;\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}-\mathrm{\&nu;}-{\mathrm{\&lambda;}}_1{\&lsqb;u_0(x,y)-c_1\&rsqb;}^2+{\mathrm{\&lambda;}}_2{\&lsqb;u_0(x,y)-c_2\&rsqb;}^2\&rsqb;---\left(2\right)$

其中，参数c₁、c₂根据下式得到，

$\begin{array}{cc}{c}_1=\frac{{\&Integral;}_{\mathrm{\&Omega;}}{u}_0(x,y)H\left(\mathrm{\&phi;}\right)dxdy}{{\&Integral;}_{\mathrm{\&Omega;}}H\left(\mathrm{\&phi;}\right)dxdy}& c_2=\frac{{\&Integral;}_{\mathrm{\&Omega;}}{u}_0(x,y)\&lsqb;1-H\left(\mathrm{\&phi;}\right)\&rsqb;dxdy}{{\&Integral;}_{\mathrm{\&Omega;}}\&lsqb;1-H\left(\mathrm{\&phi;}\right)\&rsqb;dxdy}\end{array}$

根据式(2)，采用C-V水平集演化方法分割提取建筑物区域，作为灾后影像的建筑物分割对像v_iSeg，C-V水平集为现有技术，本发明不予赘述。

步骤4，针对灾后遥感影像中步骤3所得灾后建筑物集合V_Seg相应部分提取n个特征，并检测各特征的隶属度ρ_j，结合各个特征证据概率R_j计算建筑物发生倒损的置信度，当建筑物对象v_iSeg的置信度P_Judge处于置信区间时，认为建筑物发生倒损，并入集合B_s；P_Judge处于拒绝区间时，得到未发生倒损建筑物对象，并入集合B_n；P_Judge处于不确定区间时，建筑物对象疑似发生倒损，并入集合B_d。

具体实施时，本发明技术人员可自行指定特征的数目n和具体相应特征种类。实施例将灾后遥感影像中步骤3所得灾后建筑物集合V_Seg的相应部分作为本步骤的待检测区域，提取灾后遥感影像中待检测区域的面积、纹理、高度特征，并根据不同的隶属函数检测各特征的隶属度，采用分配置信度的方式给予各个特征证据概率ρ_Area、ρ_T、ρ_Height，结合不同特征提取结果的隶属度，计算建筑物发生倒损的置信度。然后判断灾前矢量数据是否处理完毕，如果处理完毕，执行步骤5，否则，返回执行步骤3。

具体实施时，可将每个建筑物对象v_i的二维矢量面积与相应灾后遥感影像上的建筑物分割对象v_iSeg的面积进行比较，计算待检测建筑物面积差，检测面积特征证据隶属度R_iArea；提取灾后遥感影像上的建筑物分割对象v_iSeg区域内纹理特征，根据提取出的角二阶矩(AngularSecondMoment，ASM)能量值检测纹理特征证据隶属度R_iT；提取遥感影像上缓冲区内阴影特征，以此提取建筑物高度并计算其与灾前建筑物高度的差，根据高差检测高度证据隶属度R_iHeight。采用统计专家分配置信度的方式给予各个特征证据概率ρ_Area、ρ_T、ρ_Height，结合不同特征提取结果的隶属度，计算建筑物发生倒损的置信度。

实施例具体实现方式如下：

a)面积特征隶属度检测：基于面积差ΔS，根据隶属函数检测其隶属度，面积差隶属函数如式(3)所示

$R_{\mathrm{\&Delta;}S}\left(Area\right)=\left\{\begin{array}{cc}0& \mathrm{\&Delta;}S<Min\mathrm{\&Delta;}S\\ \frac{\mathrm{\&Delta;}S-Min\mathrm{\&Delta;}S}{Max\mathrm{\&Delta;}S-Min\mathrm{\&Delta;}S}& Min\mathrm{\&Delta;}S\&le;\mathrm{\&Delta;}S\&le;Max\mathrm{\&Delta;}S\\ 1& \mathrm{\&Delta;}S>Max\mathrm{\&Delta;}S\end{array}\right.---\left(3\right)$

式中，R_ΔS(Area)表示面积差特征隶属度，ΔS表示待检测建筑物灾前与灾后面积差，MinΔS表示面积差的预设最小阈值，MaxΔS表示面积差的预设最大阈值，具体实施时本领域技术人员可根据具体建筑物情况设置阈值。

即对于某建筑物对象v_i的面积S_iFore、灾后遥感影像上的相应建筑物分割对象v_iSeg面积为S_iAfter，以及面积差ΔS_i＝S_iFore-S_iAfter，面积特征提取结果为面积差ΔS_i，当ΔS_i小于给定的面积差最小阈值MinΔS时，认为建筑物发生倒损的隶属度为0；当ΔS_i大于给定的面积差最大阈值MaxΔS时，认为建筑物发生倒损的隶属度为1；当MinΔS≤ΔS≤MaxΔS，根据隶属函数确定建筑物发生倒损的隶属度

b)纹理特征隶属度检测：为提取灾后遥感影像上某区域的纹理特征，本发明采用灰度共生矩阵中ASM特征来描述纹理特征，用来检测影像同质性的ASM能量被认为是信息量最丰富的特征，未倒损建筑物ASM能量值较大，倒损建筑物值较小，计算ASM特征值，根据提取出的纹理特征值及隶属函数检测纹理特征证据的隶属度。

实施例中的具体操作步骤包括，将图像灰度值[0,255]准换到灰度区间[1,8]，计算两灰度值m、n在图像中相邻的次数P_d(m,n)，其中m、n在灰度区间内取值，构成灰度共生矩阵P_d。将各元素P_d(m,n)除以各元素之和得到归一化值ASM特征公式如式(4)所示：

$ASM=\underset{m}{\&Sigma;}\underset{n}{\&Sigma;}{\left\{{\hat{P}}_d(m,n)\right\}}^2---\left(4\right)$

式中表示各元素的归一化值。

ASM特征隶属函数如式(5)所示

$R_{ASM}\left(ASM\right)=\left\{\begin{array}{cc}1& ASM<MinASM\\ \frac{MaxASM-ASM}{MaxASM-MinASM}& MinASM\&le;ASM\&le;MaxASM\\ 0& ASM>MaxASM\end{array}\right.---\left(5\right)$

式中R_ASM(ASM)表示纹理特征隶属度，ASM表示灰度共生矩阵中角二阶矩特征值，MinASM表示预设最小ASM特征阈值，MaxASM表示预设最大ASM特征阈值。

即对于某建筑物对象v_i，提取灾后遥感影像上相应的建筑物分割对象v_iSeg区域内遥感影像纹理特征，包括根据式(4)采用灰度共生矩阵中角二阶矩特征值来描述纹理特征，相应记为ASM_i。完整的、未倒损建筑物ASM值较大，纹理杂乱的倒损建筑物ASM值较小，当待检测建筑物ASM_i大于给定的最大阈值MaxASM时，建筑物发生倒损的隶属度为0，当ASM_i小于给定的最小阈值MinASM时，认为建筑物发生倒损的隶属度为1，当MinASM≤ASM≤MaxASM时，根据隶属函数计算建筑物发生倒损的隶属度

c)高度变化特征隶属度检测，通过提取影像中建筑物阴影的方式获取灾后建筑物高度，然后计算待检测建筑物灾前与灾后高度差来检测高度变化特征的隶属度，对于某建筑物对象v_i具体实施如下：

1)由太阳方位角确定阴影相对于待检测建筑物v_i的成像方向，阴影的成像方向与太阳方位角α_sun一致，斜率为k，求出步骤3中执行所得建筑物分割对象v_iSeg的最小外接矩形Rect_i，其中一组边平行于阴影成像方向，斜率为k，记为a_i1、a_i2，另一组边垂直于阴影成像方向记为b_i1、b_i2。

根据灾前的建筑物对象v_i的高度值h_iFore、影像成像时太阳与卫星的角度参数，求出由灾前三维建筑物对象投影出的阴影长度l_i，计算公式如(6)所示：

式中Δα＝|α_sun-α_sat|，即卫星方位角与太阳方位角之间的夹角，θ_sun为太阳高度角，α_sun为太阳方位角，θ_sat为卫星传感器高度角、α_sat为卫星传感器方位角。高度值h_iFore可根据灾前三维GIS中各建筑物对象的相应三维GIS数据得到，即步骤2所获取的灾前三维GIS中建筑物对象的建筑物高度数据。

其中，参数A计算如下，

$A=(1+\frac{{\tan }^2{\mathrm{\&theta;}}_{sun}-{\tan }^2{\mathrm{\&theta;}}_{sat}}{2\&times;({\tan }^2{\mathrm{\&theta;}}_{sun}+{\tan }^2{\mathrm{\&theta;}}_{sat}-2\&times;({\tan }^2{\mathrm{\&theta;}}_{sun}\&times;{\tan }^2{\mathrm{\&theta;}}_{sat}\&times;\cos \left({\mathrm{\&alpha;}}_{sat}-{\mathrm{\&alpha;}}_{sun}\right)\left)\right)})---\left(7\right)$

将最小外接矩形Rect_i的边a_i1、a_i2沿建筑物成像方向延伸l_i，形成新的矩形Rect′_i，并以预设缓冲区距离d_Buffer2为其构建局部缓冲区Rect′_iBuf，得到待处理影像阴影区域Ω_ishadow。d_Buffer2可由本领域技术人员根据情况自行预先设定。

2)阴影有效长度计算：在1)所得Ω_ishadow内进行OSTU分割，然后对分割对象sdw_i矢量化后的轮廓曲线进行压缩处理，利用斜率为k的直线以一定的间距对压缩处理后的多边形进行分块，分别计算每条直线与多边形的两交点间的距离，即每条直线与多边形的两交点所形成线段的长度，求较长的几条线段的长度平均值作为阴影有效长度。

实施例对Ω_ishadow进行OSTU分割，自适应将待处理阴影区域分成阴影与背景两部分，用Douglas-Peucker矢量数据压缩算法对分割出的阴影区域曲线轮廓进行采样简化，形成折线图形，以2个像素为间距对简化后图形进行扫描，扫描线与a_i1、a_i2平行，记录每次斜率为k的扫描线落在图形内的长度，取最大的五个长度值的平均值，作为阴影有效长度sdw_iLen。OSTU分割和Douglas-Peucker矢量数据压缩算法为现有技术，本发明不予赘述。

3)建筑物高度计算：根据2)的sdw_iLen，成像时太阳、卫星角度参数计算建筑物高度，参数之间关系参考式(6)，只需用阴影有效长度sdw_iLen替代投影出的阴影长度l_i，即可根据阴影有效长度sdw_iLen、影像成像时太阳高度角θ_sun、方位角α_sun和卫星传感器的高度角θ_sat、方位角α_sat计算建筑物对象v_i在灾后遥感影像中的高度值h_iAfter。即按以下关系计算：

4)检测提取出的阴影特征的隶属度，高度特征隶属度函数如式(8)所示：

$R_{i\mathrm{\&Delta;}H}\left(Height\right)=\{\begin{array}{cc}0& {\mathrm{\&Delta;h}}_i<Min\mathrm{\&Delta;}H\\ \frac{{\mathrm{\&Delta;H}}_i-Min\mathrm{\&Delta;}H}{Max\mathrm{\&Delta;}H-Min\mathrm{\&Delta;}H}& Min\mathrm{\&Delta;}H\&le;{\mathrm{\&Delta;h}}_i\&le;Max\mathrm{\&Delta;}H\\ 1& {\mathrm{\&Delta;h}}_i>Max\mathrm{\&Delta;}H\end{array}---\left(8\right)$

式中，R_ΔH(Height)表示高度特征的隶属度，Δh_i表示灾前灾后建筑物高度差值，MinΔH表示高度差的预设最小值阈值，MaxΔH表示高度差的预设最大值阈值。

即对于某建筑物对象v_i，计算待检测建筑物对象灾前灾后的高度差Δh_i＝h_iFore-h_iAfter，当Δh_i小于给定的最小建筑物高度差阈值MinΔh时，建筑物发生倒损的隶属度为0，当Δh_i大于给定的最大建筑物差阈值MaxΔh时，建筑物发生倒损的隶属度为1，当MinΔH≤Δh_i≤MaxΔH时，根据隶属函数计算建筑物发生倒损的隶属度

d)根据证据理论原理，建筑物倒损各个特征证据通过分配置信度的方式给予，具体实施时本领域技术人员可根据具体情况预先设置证据概率R_j。将检测出的面积、纹理、阴影特征进行组合来计算建筑物发生倒损的置信度。认为置信度处于置信区间的对象发生倒损，并入对象集B_s，置信度处于拒绝区间的建筑物对象没有发生倒损，并入对象集B_n，处于不确定区间的建筑物有可能发生倒损，提交给人工进行处理。置信区间、拒绝区间、不确定区间可由技术人员根据实际情况预先给定。建筑物倒损置信度公式如公式(9)所示：

$P_{Judge}=\frac{\underset{j=1}{\overset{n}{\&Sigma;}}{R}_j\&times;{\mathrm{\&rho;}}_j}{n}---\left(9\right)$

其中，R_j为证据j的隶属度。ρ_j为证据j的置信度，如果没有检测到相应证据，则ρ_j取0。n为证据总个数，最终计算出建筑物发生倒损的置信度P_Judge。实施例的证据为3个，即n＝3，相应的R_j分别为面积特征证据隶属度R_iArea、纹理特征证据隶属度R_iT、高度证据隶属度R_iHeight，分别如a)、b)、c)采用相应隶属函数得到，各个特征证据概率为ρ_Area、ρ_T、ρ_Height。

步骤5，将灾后遥感影像中步骤3所得灾后建筑物集合V_Seg的相应部分作为已处理区域，其他部分作为未处理区域，未处理区域为新增建筑物可疑区，计算未倒损建筑物集合B_n区域平均灰度值，以该值为阈值确定初始水平集，对灾后遥感影像中未处理区域再次进行C-V水平集分割，剔除分割结果中面积较小的对象，本领域技术人员可自行设置面积阈值。如实施例中预设面积阈值为10m²,剔除面积小于10m²的分割对象。过滤后剩余的待检测的分割对象构成待处理对象集合B_new，其中可能包含新增建筑物。

步骤6，分别从步骤4所得倒损建筑物集合B_s、未倒损建筑物集合B_n形成的区域中随机提取像元光谱特征值作为训练样本，训练SVM分类器，利用训练出的分类器对步骤5中的B_new对象集中的对象分别进行验证，验证结果为倒损建筑物类别(新增建筑物发生倒损)的对象并入对象集B_s，构成新的倒损建筑物集合B′_s；验证为未倒损建筑物的对象(新增建筑物为发生倒损)并入对象集B_n，构成新的未倒损建筑物集合B′_n。具体实现可参见现有SVM分类器技术。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (4)

1.一种三维GIS辅助下的高分遥感影像建筑物倒损检测方法，其特征在于，包括以下步骤：

步骤1，利用灾前三维GIS或影像RPC参数对灾后遥感影像进行几何纠正；

步骤2，获取灾前三维GIS中各建筑物对象的三维GIS数据，包括二维矢量数据及高度数据，构成灾前建筑物集合V_B；

步骤3，记灾前建筑物集合V_B中每个建筑物对象为v_i，将相应二维矢量数据叠加至步骤1几何纠正后的灾后遥感影像中，在灾后遥感影像上为建筑物对象v_i构建局部缓冲区，将v_i所对应的遥感影像区域的轮廓作为初始水平集，在局部缓冲区内部的影像中进行水平集分割，得到建筑物分割对象v_iSeg；所有建筑物分割对象构成灾后建筑物集合V_Seg；

步骤4，针对灾后遥感影像中步骤3所得灾后建筑物集合V_Seg相应部分提取n个特征，并检测各特征的隶属度ρ_j，结合各个特征的预设证据概率R_j计算建筑物发生倒损的置信度，当建筑物对象v_iSeg的置信度P_Judge处于置信区间时，认为建筑物发生倒损，并入集合B_s；P_Judge处于拒绝区间时，得到未发生倒损建筑物对象，并入集合B_n；P_Judge处于不确定区间时，建筑物对象疑似发生倒损，并入集合B_d；

所述计算建筑物发生倒损的置信度采用以下公式，

$P_{Judge}=\frac{\underset{j=1}{\overset{n}{\&Sigma;}}{R}_j\&times;{\mathrm{\&rho;}}_j}{n}$

步骤5，将灾后遥感影像中步骤3所得灾后建筑物集合V_Seg的相应部分作为已处理区域，其他部分作为未处理区域；以集合B_n的区域平均灰度值作为阈值确定初始水平集，对灾后遥感影像中未处理区域进行水平集分割；滤掉分割对象中面积小于预设面积阈值的对象，剩余的对象构成集合B_new；

步骤6，分别从步骤4所得集合B_s、集合B_n形成的区域中随机提取像元光谱特征值作为训练样本训练SVM分类器，利用训练出的SVM分类器对步骤5所得集合B_new中的对象分别进行验证，将验证结果为倒损建筑物类别的对象并入集合B_s，构成新的倒损建筑物集合B′_s；将验证为未倒损建筑物类别的对象并入集合B_n，构成新的未倒损建筑物集合B′_n。

2.根据权利要求1所述的三维GIS辅助下的高分遥感影像建筑物倒损检测方法，其特征在于：步骤1中，判断输入的灾后遥感影像是否具有RPC参数，有则采用RPC参数对灾后遥感影像进行几何纠正，否则采用三维GIS对灾后遥感影像进行几何纠正。

3.根据权利要求1或2所述的三维GIS辅助下的高分遥感影像建筑物倒损检测方法，其特征在于：步骤4中，提取灾后遥感影像中步骤3所得灾后建筑物集合V_Seg相应部分的面积、纹理和高度特征，并检测各特征的隶属度，包括对每个建筑物对象v_i及灾后遥感影像上的相应建筑物分割对象v_iSeg分别进行以下处理，

a)对于某建筑物对象v_i的面积S_iFore、灾后遥感影像上相应的建筑物分割对象v_iSeg的面积为S_iAfter，面积特征提取结果为面积差ΔS_i＝S_iFore-S_iAfter，检测面积特征的隶属度如下，

$R_{i\mathrm{\&Delta;}S}\left(Area\right)=\left\{\begin{array}{cc}0& {\mathrm{\&Delta;S}}_i<Min\mathrm{\&Delta;}S\\ \frac{{\mathrm{\&Delta;S}}_i-Min\mathrm{\&Delta;}S}{Maz\mathrm{\&Delta;}S-Min\mathrm{\&Delta;}S}& Min\mathrm{\&Delta;}S\&le;{\mathrm{\&Delta;S}}_i\&le;Max\mathrm{\&Delta;}S\\ 1& {\mathrm{\&Delta;S}}_i>Max\mathrm{\&Delta;}S\end{array}\right.$

其中，MinΔS表示面积差的预设最小阈值，MaxΔS表示面积差的预设最大阈值；

b)对于某建筑物对象v_i，计算灾后遥感影像上相应的建筑物分割对象v_iSeg的灰度共生矩阵中角二阶矩特征值ASM_i作为纹理特征提取结果，检测纹理特征的隶属度如下，

$R_{iASM}\left(ASM\right)=\left\{\begin{array}{cc}1& {\mathrm{ASM}}_i<MinASM\\ \frac{MaxASM-{\mathrm{ASM}}_i}{MazASM-MinASM}& MinASM\&le;{\mathrm{ASM}}_i\&le;MaxASM\\ 0& {\mathrm{ASM}}_i>MaxASM\end{array}\right.$

其中，MinASM表示预设最小角二阶矩特征值阈值，MaxASM表示预设最大角二阶矩特征值阈值；

c)对于某建筑物对象v_i，建筑物对象v_i灾前的高度值为h_iFore，建筑物对象v_i在灾后遥感影像中的高度值h_iAfter，计算灾前灾后的高度差Δh_i＝h_iFore-h_iAfter作为高度特征提取结果，检测高度特征的隶属度如下，

$R_{i\mathrm{\&Delta;}H}\left(Height\right)=\left\{\begin{array}{cc}0& {\mathrm{\&Delta;h}}_i<Min\mathrm{\&Delta;}H\\ \frac{{\mathrm{\&Delta;h}}_i-Min\mathrm{\&Delta;}H}{Maz\mathrm{\&Delta;}H-Min\mathrm{\&Delta;}H}& Min\mathrm{\&Delta;}H\&le;{\mathrm{\&Delta;h}}_i\&le;Max\mathrm{\&Delta;}H\\ 1& {\mathrm{\&Delta;h}}_i>Max\mathrm{\&Delta;}H\end{array}\right.$

其中，MinΔH表示高度差的预设最小值阈值，MaxΔH表示高度差的预设最大值阈值。

4.根据权利要求3所述的三维GIS辅助下的高分遥感影像建筑物倒损检测方法，其特征在于：建筑物对象v_i在灾后遥感影像中的高度值h_iAfter按如下方式提取，

1)由太阳方位角α_sun确定阴影相对于建筑物对象v_i的成像方向，设斜率为k，求出相应建筑物分割对象v_iSeg的最小外接矩形Rect_i，其中一组边平行于阴影成像方向，斜率为k，记为a_i1、a_i2，另一组边垂直于阴影成像方向记为b_i1、b_i2；根据建筑物对象v_i灾前的高度值h_iFore和影像成像时太阳与卫星的角度参数，求出由建筑物对象v_i投影出的阴影长度l_i；将最小外接矩形Rect_i的边a_i1、a_i2沿建筑物成像方向延伸l_i，形成新的矩形Rect′_i，并以预设缓冲区距离d_Buffer2为其构建局部缓冲区Rect′_iBuf，得到待处理影像阴影区域Ω_ishadow；

2)在待处理影像阴影区域Ω_ishadow内进行OSTU分割，然后对所得分割对象sdw_i矢量化后的轮廓曲线进行压缩处理，利用斜率为k的直线以一定的间距对压缩处理后所得多边形进行分块，分别计算每条直线与多边形的两交点间的线段长度，求较长的多条线段的长度平均值作为阴影有效长度sdw_iLen；

3)根据阴影有效长度sdw_iLen和影像成像时太阳与卫星的角度参数，计算得到建筑物对象v_i在灾后遥感影像中的高度值h_iAfter。