CN104268879A - 基于遥感多光谱图像的建筑物实物量损毁评估方法 - Google Patents

Abstract

本发明公开了一种基于遥感多光谱图像的建筑物实物量损毁评估方法。可先对灾前灾后的多时相多光谱遥感图像做预处理，然后采用基于随机游走的变化检测方法获取二值掩膜图像，在利用二值掩膜图像对灾前图像做目标检测以提取建筑物，最后给出建筑物实物量损毁评估结果。本发明可快速、准确、自动地检测建筑物实物量损毁状况，为灾情评估、救灾决策工作提供有力的数据支持。

Description

基于遥感多光谱图像的建筑物实物量损毁评估方法

技术领域

本发明涉及遥感图像处理领域和建筑物损毁评估领域，具体地，涉及一种基于遥感多光谱图像的建筑物实物量损毁评估方法。

背景技术

我国是世界上自然灾害最为严重的国家之一。灾害种类多、分布地域广、发生频率高、危害程度深、造成损失重，给经济社会发展和人民生命财产安全带来严重影响，也对灾害相关部分快速、高效的决策提出了更高的要求。灾害评估是在救灾过程中最先启动、最重要的环节之一，它所提供的灾害评估结果，将为决策部门制定救灾、减灾应对措施提供客观依据，以有针对性地并且在最大程度上减少损失。

及时准确地对灾区建筑物损毁情况进行评估，具有十分重要的现实意义。遥感具有客观性、现实性、直观性和宏观性的优势，特别是卫星遥感，能够及时快速获取较为全面、宏观的灾后第一手资料，能够在缺乏地面调查的情况下，对灾区进行调查与损失的快速评估，为应急救助提供决策。中国从20世纪60年代开始利用遥感技术开展了对生命线工程，如道路、桥梁、电力线等损毁监测应用研究。特别是近年来，利用遥感监测进行灾后评估取得了一定成果，并建立了一系列遥感评估灾情的技术流程。目前文献主要集中在道路损毁评估上，如魏成阶、刘亚岚、陈世荣等等采用卫星遥感监测技术进行的地震和海啸灾后评估，而建筑物损毁评估方面的资料比较少。与道路损毁评估相关的文献包括Wei C J,Liu Y L,Wang S X等，Monitoring andEvaluation of"9.12"Earthquake Damage in Nantou,Taiwan Based on RemoteSensing[A].Remote Sensing in China[C].Earthquake Press 2000；Wang X Q,WeiC J,Miao C G等，The Extraction of Seismic Damage from Remote SensingImages A Case Study of BachuJiash Earthquake with Ms＝68 Occurred on Feb 24,2003[J].Earth Science Frontiers,2003,10(suppl):285-291；Chen S R,Ma H J等，Based on the high resolution remote sensing images of wenchuan earthquakedamaged roads assessment,Journal of remote sensing,2008.11；Liu Y L,ZhangY等，Wenchuan earthquake damaged roads with remote sensing monitoring andevaluation information integration,Journal of remote sensing,2008.11。

在建筑物损毁评估方面，随着IKONOS等高分辨率商业遥感卫星的应用，基于遥感图像处理技术的灾害评估发展迅速，如2003年窦爱霞采用变化检测技术对1998年张北地震和1999年台湾集集地震进行了研究，她采用基于邻域纹理提取办法提取震前、震后纹理图像，对两组纹理图像作相关系数检测变化，进行了评估震害；柳稼航等以区域为基本操作单元，提取以区域为对象的特征信息，并以综合分析的方法对目标区域从整体上进行分析，利用完好建筑与破坏建筑在区域特征和外观表现上的特征差异，自动获取震害建筑物信息。国外也有大量这方面的研究，如Saito和Spence等基于欧洲强震规模(European Macroseismic Scale，简称EMS)用多个不同的级别去评估建筑物的损毁程度，M.Pesaresi等采用超高分辨率遥感图像自动评估海啸后建筑物的损毁程度。相关文献包括窦爱霞，张景发，王晓青，等.SAR图像在震害变化检测中的应用[J].中国遥感应用协会文集[C].中国合肥，北京：宇航出版社，2003；柳稼航，杨建峰，尹京苑.遥感图像区域结构和纹理统计特性相结合的城市震害房屋自动识别－以2001年印度普吉地震和1976年唐山地震为例[J].地震学报,2004，26(6)：623-633；Saito,K.,Spence,R.J.,Going,C.,Markus,M.,2004.Using high-resolution satellite imagesfor post-earthquake building damage assessment:a study following the 26January 2001 Gujarat earthquake.Earthquake Spectra 20(1),145–170；Saito,K.,Spence,R.J.S.,2005.Visual damage assessment using high-resolution satelliteimages following the 2003 Bam,Iran,earthquake.Earthquake Spectra 21(S1),S309–S318；M.Pesaresi,A.Gerhardinger,F.Haag.Rapid damage assessment ofbuilt-up structures using VHR satellite data in tsunami-affected areas,International Journal of Remote Sensing,2007.28:3013-3016。同时也有大量的研究人员利用SAR(合成孔径雷达)数据进行震害建筑物损毁评估。如汶川地震后，刘斌涛等提出了一种结合纹理信息和边缘检测算法的震害信息提取方法，首先提取高分辨率SAR图像的差异性和逆差矩两个纹理特征，然后用边缘检测算子检测边缘信息，最后将这3个分量进行RGB合成，在此基础上进行判读，将建筑物分为完全倒塌、部分倒塌和低密度建设区。这种方法利用了影像纹理特征，但后期的目视解译工作量仍然很大。张景发等采用ERS影像分析了张北地震建筑物破坏情况，通过相关性处理、平均灰度差异性处理以及灰度方差差异性处理，获取了定量的建筑物震害信息。相关文献包括：刘斌涛，陶和平，范建容等.高分辨率SAR数据在5.12汶川地震灾害监测与评估中的应用.山地学报2008，26(3)267-271；.张景发，谢礼立，陶夏新.建筑物震害遥感图像的变化检测与震害评估.自然灾害学报，2002，11(2)59-64。

发明内容

本发明的目的是提供一种方法，该方法可快速准确地进行灾后建筑物实物量损毁评估。

为了实现上述目的，本发明提供一种基于遥感多光谱图像的建筑物实物量损毁评估方法，该方法包括：获取同一区域在第一时相(例如灾前)的第一多光谱图像和在第二时相(例如灾后)的第二多光谱图像，两个所述多光谱图像均由同样的B个波段组成，B＞1；对所述第一多光谱图像进行预处理得到第三图像，对所述第二多光谱图像进行预处理得到第四图像，所述第三图像和所述第四图像中的像素一一对应；生成所述第三图像和所述第四图像的差分图像，并对所述差分图像进行主成分分析以获取掩膜图像，所述掩膜图像用于区分所述第四图像相比所述第三图像的变化区域和非变化区域；以及在所述第三图像中对应于所述变化区域的区域内提取建筑物实物量。本发明还提供了基于随机游走算法对像素强度值位于中间区域的像素进行判断的方法。

通过上述技术方案，可对比感兴趣区域的多时相遥感图像，获取用于区分感兴趣区域中的变化部分和非变化部分的掩膜图像，然后提取变化区域中的图像特征以获取被损毁的建筑物实物量。本发明可基于主成分分析和随机游走算法进行快速准确地变化检测，并且将变化检测和建筑物检测相结合以实现全自动的建筑物实物量损毁评估。因为从低分辨率图像中难以进行建筑物提取，所以用于提取建筑物的遥感图像需要具有较高的分辨率。发明人应用四川雅安、凯里等高分辨率影像数据对本发明进行验证，均可获得良好的评估效果。本发明可为救灾、减灾工作提供有力的数据支持。

本发明的其它特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1是根据本发明的基于遥感多光谱图像对建筑物实物量损毁进行评估的流程图；

图2是根据本发明的优选实施方式的基于主成分分析和随机游走算法的变化检测的流程图；

图3(a)示出了图像的4连通图模型示意图；

图3(b)示出了图像的8连通图模型示意图；以及

图4是根据本发明的一个优选实施方式的建筑物提取流程图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

图1是根据本发明的基于遥感多光谱图像对建筑物实物量损毁进行评估的流程图。

步骤S11，可获取同一区域的多时相遥感多光谱图像。可采用如资源三号卫星等在第一时相(例如，地震前)在第二时相(例如，地震后)所拍摄的同一区域的卫星遥感图像。可将分别在第一时相和第二时相获取的图像称为第一多光谱图像和第二多光谱图像。为了进行后续的变化检测，第一多光谱图像和第二多光谱图像均由同样的B个波段组成，B＞1。

步骤S12，对多时相的多光谱数据分别进行预处理。该预处理可包括辐射校正、图像融合和图像配准等。图像校正是为了对由大气辐射等引起的图像失真进行修复。图像融合指将多光谱图像和全色图像进行融合，在保持多光谱信息的同时尽可能提高图像的空间分辨率，是否进行图像融合取决于对图像分辨率的要求和是否有相应的全色图像等因素。图像配准指对多时相图像进行配准，以使不同时相的图像中的像素彼此对应。可将第一多光谱图像和第二多光谱图像经过预处理后得到的图像分别称为第三图像和第四图像，第三图像和第四图像都是感兴趣区域的图像，并且第三图像中的像素和第四图像中的像素一一对应。

步骤S13中，对第三图像和第四图像的差分图像进行变化检测。首选，可计算第四图像和第三图像中的对应像素的差向量，生成多光谱差分图像；然后，对差分图像中的所有像素进行主成分分析(Principal ComponentsAnalysis，简称PCA)，以判断每个像素为变化状态还是非变化状态；接下来，可将所有被标记为变化状态的像素用一个值(例如，“0”)来表示，将所有被标记为非变化状态的像素用另一个值(例如，“1”)来表示，可输出感兴趣区域的掩膜图像。在感兴趣区域中，像素被标记为变化状态的区域可被称为变化区域，像素被标记为非变化状态的区域可被称为非变化区域。后文中将参考图2来详细说明根据本发明的一个优选实施方式的变化检测方法。

图1中的步骤S14中，从灾前图像(即此处的第三图像)中对应于掩膜图像的变化区域的区域中提取的建筑物。可先对第三图像中对应的变化区域进行图像分割，以得到多个区块；然后计算每个区块的特征；以及根据该特征判断每个区块是否为建筑物。后文中将参考图4来详细说明根据本发明的一个优选实施方式的建筑物提取方法。

步骤S15中，可根据被提取的建筑物，评估建筑物实物量的损毁状况。可计算被损毁的建筑物的面积和数量、被损毁的建筑物的面积占感兴趣区域的面积的百分比以及被损毁的建筑物数量占感兴趣区域中的全部建筑物数量的百分比等信息。

图2是根据本发明的优选实施方式的基于主成分分析和随机游走算法的变化检测的流程图。

在步骤S21中，可生成多光谱差分图像。第三图像中的任意像素可被表示光谱向量第四图像中的任意像素可被表示光谱向量f为波段序号，B为多光谱图像所包括的波段数，i为像素序号，N为每个图像中的像素数量，t1和t2表示图像采集时间。可通过差运算得到两幅图像中的对应像素的变化矢量从而生成感兴趣区域的差分图像。

在步骤S22中，可对用于表示差分图像的所有变化矢量进行PCA变换，用变换得到的PCA变化矢量来表示差分图像中的像素v_i(1≤i≤N)，c_b,i为像素v_i的主成分b的绝对值(本发明中也称为强度值)，b为B个主成分中的主成分的序号，B为主成分数量，i为像素序号，N为差分图像所包括的像素数量。PCA变换后得到的主成分的数量和多光谱差分图像所包括的波段数量相同。

在步骤S23中，可根据主成分的统计特性构建数学模型。针对B个主成分中的每个主成分b，设差值图像中的所有像素可按照主成分b的强度值被分为三类：变化(用ω_c,b表示)、非变化(用ω_n,b表示)、无标记(用ω_nl,b表示)。“变化”意味着主成分b的强度值较高，“非变化”意味着主成分b的强度值较低，“无标记”意味着主成分b的强度值在中间区域导致难以判断。可用下面的高斯混合模型(GMM模型)H(c_b,i)来表示依据主成分b的三类状态划分的三类像素的统计特性：

$H\left(c_{b,i}\right)=w_{n,b}N(c_{b,i},{\mathrm{\μ}}_{n,b},{\mathrm{\σ}}_{n,b}^2)+w_{\mathrm{nl},b}N(c_{b,i},{\mathrm{\μ}}_{\mathrm{nl},b},{\mathrm{\σ}}_{\mathrm{nl},b}^2)+w_{c,b}N(c_{b,i},{\mathrm{\μ}}_{c,b},{\mathrm{\σ}}_{c,b}^2),$    式(1)其中，表示主成分b为非变化状态的高斯模型，表示主成分b为无标记状态的高斯模型，表示主成分b为变化状态的高斯模型，c_b,i表示像素i的主成分b的强度值，w_n,b表示主成分b为非变化状态的高斯模型的权重，μ_n,b表示主成分b为非变化状态的高斯模型的均值，表示主成分b为非变化状态的高斯模型的方差，w_nl,b表示主成分b为无标记状态的高斯模型的权重，μ_nl,b表示主成分b为无标记状态的高斯模型的均值，表示主成分b为无标记状态的高斯模型的方差，w_c,b表示主成分b为变化状态的高斯模型的权重，μ_c,b表示主成分b为变化状态的高斯模型的均值，表示主成分b为变化状态的高斯模型的方差，可知w_n,b+w_nl,b+w_c,b＝1，可采用最大期望算法(EM算法)估计上述高斯混合模型的参数{w_n,b,μ_n,b,σ_n,b}、{w_nl,b,μ_nl,b,σ_nl,b}、{w_c,b,μ_c,b,σ_c,b}。

在步骤S24中，可进行自适应主成分选取。可根据差分图像中每个主成分的变化信息的显著性，从B个主成分中选取若干个主成分，以提高建筑物实物量损毁评估的准确性。可基于上述GMM模型选取主成分。具体步骤如下：

(1)可先通过下式计算变化指数：

$F_b=\frac{\frac{{({\mathrm{\μ}}_{c,b}-{\mathrm{\μ}}_{n,b})}^2}{{\mathrm{\σ}}_{n,b}^2}+\frac{{({\mathrm{\μ}}_{c,b}-{\mathrm{\μ}}_{\mathrm{nl},b})}^2}{{\mathrm{\σ}}_{\mathrm{nl},b}^2}}{\frac{{({\mathrm{\μ}}_{\mathrm{nl},b}-{\mathrm{\μ}}_{n,b})}^2}{{\mathrm{\σ}}_{n,b}^2}}$    式(2)

(2)按照下式对变化指数进行归一化：

$f_b=\frac{F_b}{\underset{j=1}{\overset{B}{\mathrm{\Σ}}}{F}_j}$    式(3)

(3)可按照如下规则选取主成分：选择预设门限t(0＜t＜1)，优选地，t不应小于如果则仅最大的f_b所对应的主成分被选取，主成分选取结束；否则抛弃最小的f_b所对应的主成分，并判断剩余主成分是否满足Σf_b≤t，如果满足，则未被抛弃的主成分为被选取的主成分，主成分选取结束，如果不满足，则继续抛弃剩余主成分中最小的f_b所对应的主成分，直至满足Σf_b≤t，未被抛弃的主成分为被选取的主成分。设共选取了S个主成分，各个主成分的序号用s表示，1≤S≤B。

在步骤S25中，可针对S个主成分中的每个主成分，对差分图像中的像素进行预标记。可基于式(1)所描述的GMM模型进行预标记。具体步骤如下：

(1)针对S个主成分中的每个主成分s，基于上述GMM模型计算差分图像中的每个像素的主成分s为变化状态(ω_c,s)、非变化状态(ω_n,s)和无标记状态(ω_nl,s)的概率，并将每个像素的主成分s标记为计算得到的三个概率值中最大的概率值所对应的状态。可按下式标记像素v_i的主成分s的状态ω_i,s：

${\mathrm{\ω}}_{i,s}=\underset{{\mathrm{\ω}}_{j,s}\∈({\mathrm{\ω}}_{n,s},{\mathrm{\ω}}_{c,s},{\mathrm{\ω}}_{\mathrm{nl},s})}{\arg \max }P\left({\mathrm{\ω}}_{j,s}\right)P\left(c_{i,s}\right|{\mathrm{\ω}}_{j,s})$    式(4)

其中，ω_n,s表示主成分s为非变化状态，ω_c,s表示主成分s为非变化状态，ω_nl,s表示主成分s为无标记状态，P(ω_j,s)为先验概率，P(c_i,s|ω_j,s)为将像素v_i的主成分s的强度值c_i,s代入对应的高斯模型后得到的概率值，例如，P(c_i,s|ω_n,s)即为将c_i,s代入后得到的概率值。设经过上述标记后，差分图像中主成分s被标记为ω_c,s的所有像素的位置的集合为主成分s被标记为ω_n,s的所有像素的位置的集合为

(2)消除差分图像噪声，并对滤波后的图像进行标记。PCA变化虽然去除了差分图像中的一部分噪声，但是主成分中的剩余噪声仍然会影响变化检测的精度。可通过例如高斯栈(Gaussian Stack)滤波的方法有效消除噪声。可将用上述差分图像可作为高斯栈的第0层，用Cs表示由差分图像的所有像素的主成分s的强度值所构成的矩阵，则高斯栈的第n层可表示为：

$C_s^n=G(w,\mathrm{\&sigma;})*C_s^{n-1},1<n\&le;l$    式(5)

其中，G(w,σ)为高斯滤波器，w为该高斯滤波器的大小，σ为高斯标准差。可根据经验选取w、σ。可指定高斯栈的层数为l，l可在2至4间选取。高斯栈顶层可被记为其所表示的图像可被称为第五图像。可知第五图像和差分图像中的像素一一对应。对于S个主成分中的每个主成分s，可以以与步骤23类似的方法建立用于表示第五图像中的主成分s的统计特性的GMM模型，并按照式(4)所示的方式对第五图像中的每个像素进行标记。设经过标记后，第五图像中主成分s被标记为ω_c,s的所有像素的位置的集合为主成分s被标记为ω_n,s的所有像素的位置的集合为

(3)综合步骤S25中(1)和(2)的标记结果，重新标记差分图像中的像素的主成分状态。针对每个主成分s，位置上的像素的主成分s的状态被标记为ω_c,s，位置上的像素的主成分s的状态被标记为ω_n,s，其他位置上的像素被标记为ω_nl,s，换言之，针对每个主成分s，针对差分图像中的每个像素v_i，只有像素v_i的主成分s被标记为ω_c,s，并且其在第五图像中的对应像素的主成分s也被标记为ω_c,s，才将像素v_i的主成分s标记为ω_c,s；只有像素v_i的主成分s被标记为ω_n,s，并且其在第五图像中的对应像素的主成分s也被标记为ω_n,s，像素v_i的主成分s才被标记为ω_n,s；否则，该像素的s主成分被标记为ω_nl,s。

在步骤S26中，可基于随机游走算法将差分图像中主成分的无标记状态判定为变化状态或不确定状态。

(1)针对S个主成分中的每个主成分s，计算差分图像的相邻像素权重矩阵W_s，W_s中第i行第j列的元素为：

   式(6)

其中，v_i表示差分图像中的第i个像素，1≤i≤N，c_i,s表示像素v_i的主成分s的强度值，式(6)中的“相邻”指4连通图模型(如图3(a)所示)中的相邻或8连通图模型(如图3(b)所示)中的相邻。图3(a)中，像素5共有四个相邻元素，分别为像素2、4、6和8；图3(b)中，像素5共有8个相邻元素，分别为像素1、2、3、4、6、7、8和9。

(2)针对S个主成分中的每个主成分s，计算差分图像的Laplacian(拉普拉斯)矩阵L_s，L_s中的第i行第j列的元素为：

   式(7)

其中， $d_{s,i}=\underset{1\&le;j\&le;N}{\mathrm{\&Sigma;}}{w}_{s,i,j}.$

(3)针对S个主成分中的每个主成分s，构建Dirichlet(狄利克雷)模型：

$D\left[u_s\right]=\frac{1}{2}{u_s}^T{L}_s{u}_s$    式(8)

其中，u_s(v_i)表示像素v_i的主成分s为变化状态的概率， $u_s=\left[\begin{array}{c}{u}_s\left(v_1\right)\\ u_s\left(v_2\right)\\ ...\\ u_s\left(v_N\right)\end{array}\right].$

(4)求解上述Dirichlet模型。下面将具体介绍本优选实施方式中Dirichlet模型的求解过程：

对 $u_s=\left[\begin{array}{c}{u}_s\left(v_1\right)\\ u_s\left(v_2\right)\\ ...\\ u_s\left(v_N\right)\end{array}\right]$ 中的元素重新排序，将主成分s已被标记为ω_c,s和ω_n,s的像素所对应的概率排在前面，将主成分s已被标记为ω_nl,s的像素所对应的概率排在后面，重排后得到 $u_s=\left[\begin{array}{c}{u}_{l,s}\\ u_{n,s}\end{array}\right],$ 其中u_l,s是由主成分s被标记为ω_c,s和ω_n,s的像素所对应的概率组成的向量，u_n,s是由主成分s被标记为ω_nl,s的像素所对应的概率组成的向量；

将 $D\left[u_s\right]=\frac{1}{2}{u_s}^T{L}_s{u}_s$ 分解为：

$D[u_{l,s},u_{n,s}]=\frac{1}{2}\left[\begin{array}{cc}{u_{l,s}}^T& {u_{n,s}}^T\end{array}\right]\left[\begin{array}{cc}{L}_{l,s}& R_s\\ {R_s}^T& L_{n,s}\end{array}\right]\left[\begin{array}{c}{u}_{l,s}\\ u_{n,s}\end{array}\right]$    式(9)

$=u_{l,s}^T{L}_{l,s}{u}_{l,s}+2u_{n,s}^T{R_s}^T{u}_{l,s}+u_{n,s}^T{L}_{n,s}{u}_{n,s}$

其中，L_s被相应地分解为L_l,s、L_n,s、R_s；

然后求解得到：

u_n,s＝-L_n,s ^-1R_s ^Tu_l,s   式(10)

从而得到差分图像中主成分s被标记为无标记状态的所有像素的主成分为变化状态的概率。

(5)根据计算得到的u_n,s对主成分被标记为ω_nl,s的像素的主成分s进行判定。预设概率可为0.5。对于任意在步骤S25中主成分s被标记为ω_nl,s的像素v_i，将像素v_i对应的概率u_s(v_i)和预设概率进行比较。如果u_s(v_i)＞0.5，则判定像素v_i的主成分s为变化状态，并将其标记为ω_c,s，否则判定像素v_i的主成分s为非变化状态，并将其标记为ω_n,s。

在步骤S27中，可综合S个主成分的变化检测结果，生成掩膜图像。对于差分图像中的每个像素v_i，如果像素v_i的S个主成分中任意一个主成分被标记为变化状态，则将像素v_i标记为变化状态ω_c；否则，将像素v_i标记为非变化状态ω_n，以得到掩膜图像。该掩膜图像为二值掩膜图像，其中，对应于ω_c的区域可被称为变化区域，对应于ω_n的区域可被称为非变化区域。

图4是根据本发明的一个优选实施方式的建筑物提取流程图。

在步骤S41中，可从第三图像(例如，灾前图像)中对应于掩膜图像的变化区域的区域中去除植被部分和水体部分。具体步骤为：

采用式(11)和式(12)计算第三图像中对应于掩膜图像中的变化区域的区域中的各个像素的归一化植被指数(NDVI)和归一化水体指数(NDWI)，其中，ρ_nir表示近红外反射辐射通量；ρ_green和ρ_red分别表示绿色光和红色光的反射辐射通量；

$\mathrm{NDVI}=\frac{{\mathrm{\&rho;}}_{\mathrm{nir}}-{\mathrm{\&rho;}}_{\mathrm{red}}}{{\mathrm{\&rho;}}_{\mathrm{nir}}+{\mathrm{\&rho;}}_{\mathrm{red}}}$    式(11)

$\mathrm{NDWI}=\frac{{\mathrm{\&rho;}}_{\mathrm{green}}-{\mathrm{\&rho;}}_{\mathrm{nir}}}{{\mathrm{\&rho;}}_{\mathrm{green}}+{\mathrm{\&rho;}}_{\mathrm{red}}}$    式(12)

可将NDVI和NDWI的判断阈值都设为0，则NDVI＞0的区域可被判断为植被区域，NDVI≤0的区域可被判断为非植被区域；NDWI＞0的区域可被判断为水体区域，NDWI≤0的区域可被判断为非水体区域。然后去除对应于变化区域的区域中的植被部分和水体部分。

在步骤S42中，可对第三图像中对应于掩膜图像的变化区域的区域进行平滑滤波，以去除背景噪声。本实施方式中选择双边滤波算法进行平滑滤波。经双边滤波后，图像的前景边缘基本保持不变，而背景变得较为平滑。

根据本发明的另一实施方式，也可先对第三图像进行平滑滤波，再去除植被部分和/或水体部分。

在步骤S43中，可在经过上述步骤S41和S42处理后的第三图像的变化区域中，进行图像分割。本实施方式中，采用Mean Shift算法进行图像分割，以得到多个区块。

还可在经分割后的变化区域中进行连通域标记，以将多个区块分为不同类别。

在步骤S44中，可计算各个区块的特征。本实施方式中，所采用的特征可包括几何特征、颜色特征和纹理特征。具体介绍如下：

(1)几何特征。本实施方式可采用形状指数、矩形指数、分散度表示每个区块的几何特征，三种指数的计算方法如下：

$S=\frac{4\sqrt{A}}{P}$    式(13)

$R=\frac{A}{A_{\min \mathrm{irect}}}$    式(14)

$D=\frac{P^2}{A}$    式(15)

其中S表示形状指数，A表示区块面积，P表示区块周长；R表示矩形指数，A_minirect表示区块的最小外接矩形的面积，区块越接近矩形，R值越接近于1；D表示分散度指数，是度量区块紧凑度的指标。

此外，由于建筑物的面积大小具有一定的范围，此处还可以通过设置阈值去除面积过大(如一些农田、树林等背景区块)或过小(如汽车)的区块，提高提取效率。本实施方式中可将建筑物的最小面积设为20平方米，最大面积设为200平方米。

(2)颜色特征。遥感卫星探测器检测的是地物发射或辐射特定波长电磁波的强度，在遥感影像中是以灰度值形式显示的。为了更好地区分建筑物与周围地物的颜色特征，可对图像的颜色空间进行变换。基于此，本实施方式采用Luv颜色空间。

Luv颜色空间是对CIE XYZ空间进行简单变换得到的，具有视觉统一性。其中，用L*表示物体的亮度，用u*和v*表示色度。对于一般图像，u*和v*的取值范围为-100至+100，L*的取值范围为0至100。对于一幅RGB图像，可以先将其转换至XYZ空间，再将其转换至Luv颜色空间。

从RGB空间转到XYZ空间可以通过下式计算(转换之前需要将RGB值归一化到0到1之间)：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}0.412453& 0.357580& 0.180423\\ 0.212671& 0.715160& 0.072169\\ 0.019334& 0.119193& 0.950227\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]$    式(16)

从XYZ空间转到Luv空间可以通过下式计算：

$L^{*}=\left\{\begin{array}{cc}\frac{{\left(\frac{29}{3}\right)}^{3}Y}{Y_n},& \frac{Y}{Y_n}\&le;{\left(\frac{6}{29}\right)}^3\\ 116{\left(\frac{Y}{Y_n}\right)}^{\frac{1}{3}}-16,& \frac{Y}{Y_n}\&rang;{\left(\frac{6}{29}\right)}^3\end{array}\right.$

u^*＝13L^*·(u'-u'_n)

                                        式(17)

v^*＝13L^*·(v'-v'_n)

其中，u'_n＝0.19793943和v'_n＝0.46831096是色度图上一个白色点的色度坐标，Y_n是CIE XYZ系统中参考白光的Y值，通常取1。u'和v'可以通过下式计算：

$\begin{array}{c}{u}_n^{\&prime;}=\frac{4X}{X+15Y+3Z}=\frac{4x}{-2x+12y+3}\\ v_n^{\&prime;}=\frac{9X}{X+15Y+3Z}=\frac{9x}{-2x+12y+3}\end{array}$    式(18)

在Luv色彩空间中，可以直接通过计算两点坐标(L^*,u^*,v^*)的欧几里得距离得到该两点之间的色差。

转换得到的Luv图像中每个像素的L、u、v值已被归一化到0至255的范围内。可将这一范围分为若干等长的区段，本实施方式中将0-255分成26个区段，对L、u、v分量，分别统计各分割区块中所有像素落在每个区段的数量，从而得到各分割区块的L、u、v颜色分布统计直方图。

(3)纹理特征。局部二值模式(LBP)是一种对光照具有不变性的纹理描述因子。LBP通过记录一个像素与周围其他像素点的对比信息，描述该像素周围邻域的纹理特征信息。LBP特征的计算过程如下：选中图像中的一个像素，将该像素的灰度值设为阈值；比较该像素相邻的8个像素的灰度值和阈值，若某个相邻像素的灰度值大于所设阈值，则将该像素标记为1，否则标记为0。将这8个像素的标记值按顺时针排列得到一个8位的二进制序列，该二进制序列对应的十进制数就是被选中的像素的LBP值。可计算出各区块中每个像素的LBP值，并统计其直方图，即可得出各区块的纹理特征信息

本优选实施方式中采用LBP特征的一个变种，称为Uniform LBP(简称ULBP)。在一个二进制模式中，将二进制序列首尾相连，如果该循环二进制数从0到1或从1到0最多有两次跳变(例如00000000(没有跳变)、11000011(两次跳变))，则该二进制模式是一个等价模式。对于旋转不变的LBP，可以将等价模式根据二进制数中1的数目分为9类。对于3×3邻域，等价模式的二进制序列可以占到全部序列的90％，因此可以把剩下的非等价模式的序列归为一类，这样得出的ULBP二进制模式就被归为10类。

同之前颜色特征的表示类似，计算各分割区块中所有像素的LBP值，判断该像素属于上述10类中的哪一类，然后统计每类中的像素数量，就可以得到各分割区块的ULBP类别分布统计直方图。再对其进行归一化处理，即将各类别的统计数值除以该分割区块像素数，使得10个类别的数值之和为1。最终得到的就是该区块纹理的一个10维特征向量。

在步骤S45中，检测各个区块的特征以判断其中的哪些区块为建筑物。可采用支持向量机(SVM)、adaboost等分类器对各个区块进行检测。本实施方式中可采用SVM分类器。可先对SVM进行样本训练，并将训练结果设置为建筑物和非建筑物两类；然后将计算得到的各个区块的特征输入SVM分类器，以判断每个区块是否为建筑物。

然后，即可根据检测结果进行建筑物实物量损毁评估。

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (10)

1.一种基于遥感多光谱图像的建筑物实物量损毁评估方法，该建筑物实物量损毁评估方法包括：

获取同一区域在第一时相的第一多光谱图像和在第二时相的第二多光谱图像，两个所述多光谱图像均由同样的B个波段组成，B＞1；

对所述第一多光谱图像进行预处理得到第三图像，对所述第二多光谱图像进行预处理得到第四图像，所述第三图像和所述第四图像中的像素一一对应；

生成所述第三图像和所述第四图像的差分图像，并对所述差分图像进行主成分分析以获取掩膜图像，所述掩膜图像用于区分所述第四图像相比所述第三图像的变化区域和非变化区域；以及

在所述第三图像中对应于所述变化区域的区域内提取建筑物实物量。

2.根据权利要求1所述的建筑物实物量损毁评估方法，其中，对所述差分图像进行主成分分析以获取掩膜图像包括：

获取所述差分图像的一个或多个主成分；

将所述差分图像中的每个像素的所述一个或多个主成分中的每个主成分s的状态标记为变化状态、非变化状态和无标记状态中的一者，s用于表示所述一个或多个主成分中的主成分的序号；

针对所述一个或多个主成分中的每个主成分s，将所述差分图像中主成分s的状态被标记为无标记状态的所有像素的主成分s的状态判定为变化状态或者非变化状态；以及

对于所述差分图像中的每个像素，如果该像素的任意一个所述主成分的状态为变化状态，则将该像素标记为变化状态，否则将该像素标记为非变化状态，以获取所述掩膜图像。

3.根据权利要求2所述的建筑物实物量损毁评估方法，其中，获取所述差分图像的一个或多个主成分包括：

对所述差分图像进行主成分分析变换，得到用于表示所述差分图像的B个主成分；以及

根据所述差分图像中所述B个主成分中的每个主成分的变化信息的显著性从所述B个主成分中选取所述一个或多个主成分，包括：

(1)针对所述B个主成分中的每个主成分b，设所述差分图像中的所有像素按照主成分b的强度值被分为三类：变化、非变化、以及无标记，并建立如下形式的高斯混合模型H(c_b,i)表示三类像素的统计特性：

$H\left(c_{b,i}\right)=w_{n,b}N(c_{b,i},{\mathrm{\&mu;}}_{n,b},{\mathrm{\&sigma;}}_{n,b}^2)+w_{\mathrm{nl},b}N(c_{b,i},{\mathrm{\&mu;}}_{\mathrm{nl},b},{\mathrm{\&sigma;}}_{\mathrm{nl},b}^2)+w_{c,b}N(c_{b,i},{\mathrm{\&mu;}}_{c,b},{\mathrm{\&sigma;}}_{c,b}^2),$

其中，表示主成分b为非变化状态的高斯模型，表示主成分b为无标记状态的高斯模型，表示主成分b为变化状态的高斯模型，b用于表示所述B个主成分中的主成分的序号，c_b,i表示像素vi的主成分b的强度值，w_n,b表示主成分b为非变化状态的高斯模型的权重，μ_n,b表示主成分b为非变化状态的高斯模型的均值，表示主成分b为非变化状态的高斯模型的方差，w_nl,b表示主成分b为无标记状态的高斯模型的权重，μ_nl,b表示主成分b为无标记状态的高斯模型的均值，表示主成分b为无标记状态的高斯模型的方差，w_c,b表示主成分b为变化状态的高斯模型的权重，μ_c,b表示主成分b为变化状态的高斯模型的均值，表示主成分b为变化状态的高斯模型的方差；

(2)针对所述B个主成分中的每个主成分b，通过下式计算变化指数：

$F_b=\frac{\frac{{({\mathrm{\&mu;}}_{c,b}-{\mathrm{\&mu;}}_{n,b})}^2}{{\mathrm{\&sigma;}}_{n,b}^2}+\frac{{({\mathrm{\&mu;}}_{c,b}-{\mathrm{\&mu;}}_{\mathrm{nl},b})}^2}{{\mathrm{\&sigma;}}_{\mathrm{nl},b}^2}}{\frac{{({\mathrm{\&mu;}}_{\mathrm{nl},b}-{\mathrm{\&mu;}}_{n,b})}^2}{{\mathrm{\&sigma;}}_{n,b}^2}};$

并按照下式归一化所述变化指数：

$f_b=\frac{F_b}{\underset{j=1}{\overset{B}{\mathrm{\&Sigma;}}}{F}_j}$

；以及

(3)按照如下规则选取所述一个或多个主成分：如果 则仅最大的f_b所对应的主成分被选取，主成分选取结束；否则抛弃最小的f_b所对应的主成分，并判断剩余主成分是否满足Σf_b≤所述预设门限，如果满足，则未被抛弃的主成分为被选取的主成分，主成分选取结束，如果不满足，则继续抛弃剩余主成分中最小的f_b所对应的主成分，直至满足Σf_b≤所述预设门限，未被抛弃的主成分为被选取的主成分，其中，0＜所述预设门限＜1。

4.根据权利要求2所述的建筑物实物量损毁评估方法，其中将所述差分图像中的每个像素的所述一个或多个主成分中的每个主成分s的状态标记为变化状态、非变化状态和无标记状态中的一者包括：

针对所述一个或多个主成分中的每个主成分s，设所述差分图像中的所有像素按照主成分s的强度值被分为三类：变化、非变化、以及无标记，并建立数学模型表示主成分s的三类像素的统计特性；以及

针对所述一个或多个主成分中的每个主成分s，基于主成分s的所述数学模型，分别计算所述差分图像中每个像素的主成分s为变化状态、非变化状态和无标记状态的概率值，并将每个像素的主成分s标记为计算得到的三个概率值中最大的概率值所对应的状态。

5.根据权利要求4所述的建筑物实物量损毁评估方法，还包括：

对所述差分图像进行滤波，得到第五图像，用所述一个或多个主成分表示所述第五图像，所述第五图像和所述差分图像中的像素一一对应；

针对所述一个或多个主成分中的每个主成分s，设所述第五图像中的所有像素按照主成分s的强度值被分为三类：变化、非变化、以及无标记，并建立数学模型表示主成分s的三类像素的统计特性；

针对所述一个或多个主成分中的每个主成分s，基于主成分s的所述数学模型，分别计算所述第五图像中每个像素的主成分s为变化状态、非变化状态和无标记状态的概率，并将该像素的主成分s标记为计算得到的三个概率值中最大的概率值所对应的状态；以及

重新标记所述差分图像中的像素的主成分的状态，针对所述一个或多个主成分中的每个主成分s，针对所述差分图像中的每个像素，如果该像素的主成分s被标记为变化状态，并且所述第五图像中的对应像素的主成分s也被标记为变化状态，则将该像素的主成分s标记为变化状态；如果该像素的主成分s被标记为非变化状态，并且所述第五图像中的对应像素的主成分s也被标记为非变化状态，则将该像素的主成分s标记为非变化状态；否则将该像素的主成分s标记为无标记状态。

6.根据权利要求4或5所述的建筑物实物量损毁评估方法，其中，所述数学模型都为高斯混合模型。

7.根据权利要求2所述的建筑物实物量损毁评估方法，其中，将所述差分图像中主成分s的状态被标记为无标记状态的所有像素的主成分s的状态判定为变化状态或者非变化状态是基于随机游走算法。

8.根据权利要求7所述的建筑物实物量损毁评估方法，其中，基于随机游走算法将已被标记的无标记状态判定为变化状态或非变化状态包括：

针对所述一个或多个主成分中的每个主成分s，计算所述差分图像的相邻像素权重矩阵W_s，W_s中第i行第j列的元素为：

像素v_i和像素v_j相邻

其中，c_s,i表示像素v_i的主成分s的强度值；

针对所述一个或多个主成分中的每个主成分s，计算所述差分图像的Laplacian矩阵L_s，L_s中的第i行第j列的元素为：

其中， $d_{s,i}=\underset{1\&le;j\&le;N}{\mathrm{\&Sigma;}}{w}_{s,i,j};$

针对所述一个或多个主成分中的每个主成分s，构建Dirichlet模型：

$D\left[u_s\right]=\frac{1}{2}{u_s}^T{L}_s{u}_s$

其中，u_s(v_i)表示像素v_i的主成分s为变化状态的概率， $u_s=\left[\begin{array}{c}{u}_s\left(v_1\right)\\ u_s\left(v_2\right)\\ ...\\ u_s\left(v_N\right)\end{array}\right];$ 以及

求解所述Dirichlet模型，得到 $u_s=\left[\begin{array}{c}{u}_s\left(v_1\right)\\ u_s\left(v_2\right)\\ ...\\ u_s\left(v_N\right)\end{array}\right],$ 对于其中任意已被标记为无标记状态的像素v_i，如果u_s(v_i)＞预设概率，则判定该像素v_i的主成分s为变化状态；否则，判定该像素v_i的主成分s为非变化状态。

9.根据权利要求8所述的建筑物实物量损毁评估方法，其中，求解所述Dirichlet模型包括：

对 $u_s=\left[\begin{array}{c}{u}_s\left(v_1\right)\\ u_s\left(v_2\right)\\ ...\\ u_s\left(v_N\right)\end{array}\right]$ 中的元素重新排序，将主成分s已被标记为变化状态和非变化状态的像素所对应的概率排在前面，将主成分s已被标记为无标记状态的像素所对应的概率排在后面，即 $u_s=\left[\begin{array}{c}{u}_{l,s}\\ u_{n,s}\end{array}\right],$ 其中u_l,s是由主成分s已被标记为变化状态和非变化状态的像素所对应的概率组成的向量，u_n,s是由主成分s已被标记为无标记状态的像素所对应的概率组成的向量；

将 $D\left[u_s\right]=\frac{1}{2}{u_s}^T{L}_s{u}_s$ 分解为：

$D[u_{l,s},u_{n,s}]=\frac{1}{2}\left[\begin{array}{cc}{u_{l,s}}^T& {u_{n,s}}^T\end{array}\right]\left[\begin{array}{cc}{L}_{l,s}& R_s\\ {R_s}^T& L_{n,s}\end{array}\right]\left[\begin{array}{c}{u}_{l,s}\\ u_{n,s}\end{array}\right]$

其中，L_s被相应地分解为L_l,s、L_n,s、R_s；以及

求解得到：

u_n,s＝-L_n,s ^-1R_s ^Tu_l,s。

10.根据权利要求1所述的建筑物实物量损毁评估方法，其中在所述第三图像中对应于所述变化区域的区域内提取建筑物实物量包括：

对所述第三图像中对应于所述变化区域的区域进行图像分割，以得到多个区块；

计算所述多个区块中每一者的特征；以及

对所述多个区块中每一者的所述特征进行检测，以判断所述多个区块中的每个区块是否为建筑物。