一种基于高分辨率遥感数据的地震损毁评估方法和装置与流程

本发明涉及快速评估技术领域，特别是指一种基于高分辨率遥感数据的地震损毁评估方法和装置。

背景技术：

地震灾害受灾范围快速评估是指在地震灾害刚刚发生后的(通常为48小时以内)，在灾区地面调查结果有限或者不足的情况下，使用灾区获取的实时遥感影像数据快速获取本次地震灾害中受影响区域范围及受影响程度。在近年来的地震灾害后评估实际研究中，已经表明，城市范围的地震损毁情况可以通过卫星影像或者高分辨率的遥感影像分析得到。

相关研究已经证实基于变化监测和高分辨率遥感影像的方法可以有效评估地震灾害中的损毁情况，例如Keiko Saito,等针对通过遥感影像目视解译结果与地面数据对比验证了高分辨率卫星影像对基于区域的损毁评估的可靠性。Charles K.Huyck等使用QuickBird数据和Neighborhood Edge Dissimilarities特征方法，用于评估bam地震再次高分辨率数据用于快速评估的可行性。

但是在现有技术中，基于遥感的损毁评估图只能通过监测灾害前后的遥感影像变化进行绘制，而基于卫星的损毁评估图只能通过卫星数据。也就是说，由于卫星遥感影像数据不同传感器生成的影像空间分辨率、配准偏差、传感器入射角度影像的不同，会影响变化检测方法结果的准确性,现有方法只能灾害前后的两影像数据需要同种传感器获得。

技术实现要素：

有鉴于此，本发明的目的在于提出一种基于高分辨率遥感数据的地震损毁评估方法和装置，解决现有对地震损毁评估只能使用卫星影像或者航空影像的问题。

基于上述目的本发明提供基于高分辨率遥感数据的地震损毁评估方法，包括步骤：

根据监测边缘特征变化，获得地震前后遥感影像的梯度图像；

分别获得所述地震前后梯度图像的相似度特征图像，提取所述相似度图像的相似度特征；

根据评估倒毁率绘制当地地震灾害快速评估图。

在本发明的一些实施例中，所述根据监测边缘特征变化，获得地震前后遥感影像的梯度图像，包括：

获取灾害地区地震前后遥感影像，并进行地理配准；

根据地理配准后的所述地震前后遥感影像，分别计算地震前后遥感影像的梯度图像。

在本发明的一些实施例中，通过如下公式，分别计算地震前后遥感影像的梯度图像：

$G(x,y)=\sqrt{{\left(f\right(x,y)-f(x+1,y+1\left)\right)}^2+{\left(f\right(x+1,y)-f(x,y+1\left)\right)}^2}$

$G(x^{\′},y^{\′})=\frac{G(x,y)-G_{\min }(x,y)}{G_{max}(x,y)-G_{min}(x,y)}$

其中，f(x,y)为源图像的单通道，G(x',y')为梯度图像。

在本发明的一些实施例中，分别获得所述地震前后梯度图像的相似度特征图像，提取所述相似度图像的相似度特征，包括：

将地震前后遥感影像的梯度图像网格化；

分别计算网格化后的地震前后梯度图像的相似度特征图像；

根据所述地震前后梯度图像的相似度图像的相似度特征，获取完全变化和严重变化网格数量。

在本发明的一些实施例中，根据如下公式，计算网格化后的地震前后梯度图像的相似度特征图像：

$S\left(W\right)=\frac{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)*\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}{\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)}^2}*\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}^2}}$

其中，G_A(x，y)表示T1时间(灾前)网格W的梯度特征图像；G_B(x，y)表示T2时间灾后网格W的梯度特征图像；

还有，G_A-mean(x，y)＝(∑∑G_A(x，y))/H，表示目标A区域中像元的梯度平均值。G_B-mean(x，y)＝(∑∑G_B(x，y))/H，表示目标B区域中像元的梯度平均值；H则为网格中像元的总个数，k和l分别为目标区域的宽度和高度。

在另一方面，本发明还提供了一种基于高分辨率遥感数据的地震损毁评估装置，包括：

梯度图像获取单元，用于根据监测边缘特征变化，获得地震前后遥感影像的梯度图像；

相似度特征获取单元，用于分别获得所述地震前后梯度图像的相似度图像，提取所述相似度图像的相似度特征；

评估单元，用于根据评估倒毁率绘制当地地震灾害快速评估图。

在本发明的一些实施例中，梯度图像获取单元，还用于：

获取灾害地区地震前后遥感影像，并进行地理配准；

根据地理配准后的所述地震前后遥感影像，分别计算地震前后遥感影像的梯度图像。

在本发明的一些实施例中，通过如下公式，分别计算地震前后遥感影像的梯度图像：

$G(x,y)=\sqrt{{\left(f\right(x,y)-f(x+1,y+1\left)\right)}^2+{\left(f\right(x+1,y)-f(x,y+1\left)\right)}^2}$

$G(x^{\′},y^{\′})=\frac{G(x,y)-G_{\min }(x,y)}{G_{max}(x,y)-G_{min}(x,y)}$

其中，f(x,y)为源图像的单通道，G(x',y')为梯度图像。

在本发明的一些实施例中，所述相似度特征获取单元，还用于：

将地震前后遥感影像的梯度图像网格化；

分别计算网格化后的地震前后梯度图像的相似度特征图像；

根据所述地震前后梯度图像的相似度图像的相似度特征，获取完全变化和严重变化网格数量。

在本发明的一些实施例中，根据如下公式，计算网格化后的地震前后梯度图像的相似度特征图像：

$S\left(W\right)=\frac{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)*\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}{\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)}^2}*\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}^2}}$

其中，G_A(x，y)表示T1时间(灾前)网格W的梯度特征图像；G_B(x，y)表示T2时间灾后网格W的梯度特征图像；

还有，G_A-mean(x，y)＝(∑∑G_A(x，y))/H，表示目标A区域中像元的梯度平均值。G_B-mean(x，y)＝(∑∑G_B(x，y))/H，表示目标B区域中像元的梯度平均值；H则为网格中像元的总个数，k和l分别为目标区域的宽度和高度。

从上面所述可以看出，本发明提供的一种基于高分辨率遥感数据的地震损毁评估方法和装置，通过基于边缘相似度的变化监测分析和倒毁率评估，实现从像素到对象再到评估区逐级提取损毁信息的评估，达到了对建成区灾害损失快速评估，及倒损率的初步判定。还有，本发明能够同时应用于卫星影像与航空影像，且可以兼容用于变化监测的灾害前后影像为异源传感器。

附图说明

图1为本发明第一实施例中基于高分辨率遥感数据的地震损毁评估方法的流程示意图；

图2为本发明可参考实施例中基于高分辨率遥感数据的地震损毁评估方法的流程示意图；

图3为本发明实施例中基于高分辨率遥感数据的地震损毁评估装置的结构示意图；

图4为本发明实施例中地震灾害前后遥感影像数据图；

图5为本发明实施例中图4基础上的梯度图像；

图6为本发明实施例中图5基础上的相似度特征图像；

图7为本发明实施例中基于图6的变化检测方法分类结果图；

图8为本发明实施例中图7基础上的倒毁率评估图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

需要说明的是，本发明实施例中所有使用“第一”和“第二”的表述均是为了区分两个相同名称非相同的实体或者非相同的参量，可见“第一”“第二”仅为了表述的方便，不应理解为对本发明实施例的限定，后续实施例对此不再一一说明。

参阅图1所示，为本发明实施例中基于高分辨率遥感数据的地震损毁评估方法的流程示意图，所述基于高分辨率遥感数据的地震损毁评估方法包括：

步骤101，根据监测边缘特征变化，获得地震前后遥感影像的梯度图像。

作为实施例，可以根据地震前后遥感影像，分别计算地震前后遥感影像的梯度图像。较佳地，通过如下公式，分别计算地震前后遥感影像的梯度图像：

$G(x,y)=\sqrt{{\left(f\right(x,y)-f(x+1,y+1\left)\right)}^2+{\left(f\right(x+1,y)-f(x,y+1\left)\right)}^2}$

$G(x^{\′},y^{\′})=\frac{G(x,y)-G_{\min }(x,y)}{G_{max}(x,y)-G_{min}(x,y)}$

其中，f(x,y)为源图像的单通道，G(x',y')为梯度图像。

优选地，在计算地震前后遥感影像的梯度图像之前，需要对获得的灾害地区地震前后遥感影像进行地理配准。其中，可以利用Arcgis软件将获取的地震前后遥感影像进行地理配准。

另外，还可以对地震前后遥感影像的梯度图像进行网格化。较佳地，可以根据当地评估目标(建筑物)的平均长度L，以及获取的遥感影像空间分辨率R，选择网格大小。优选地，采用香浓采样原理，选取网格大小。在该实施例中，接近S＝0.5*L/R，当然根据实际情况可以进行调整。

步骤102，分别获得所述地震前后梯度图像的相似度图像，提取所述相似度图像的相似度特征。

在该实施例中，计算网格化后的地震前后梯度图像的相似度图像，通过如下公式：

$S\left(W\right)=\frac{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)*\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}{\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)}^2}*\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}^2}}$

其中，G_A(x，y)表示T1时间(灾前)网格W的梯度特征图像。G_B(x，y)表示T2时间(灾后)网格W的梯度特征图像。还有，G_A-mean(x，y)＝(∑∑G_A(x，y))/H，表示目标A区域中像元的梯度平均值。G_B-mean(x，y)＝(∑∑G_B(x，y))/H，表示目标B区域中像元的梯度平均值。H则为网格中像元的总个数，k和l分别为目标区域的宽度和高度。

步骤103，根据评估倒毁率绘制当地地震灾害快速评估图。

在本发明的一个可参考的实施例中，如图2所示，所述的基于高分辨率遥感数据的地震损毁评估方法包括：

步骤201，获取灾害地区地震前后遥感影像，并进行地理配准。

其中，可以利用Arcgis软件将获取的地震前后遥感影像进行地理配准。在实施例中，地理配准是将前后的两景影像按照相同投影方式，投影到同一坐标系中，该过程是进行地理计算的基础。

步骤202，根据地理配准后的所述地震前后遥感影像，分别计算地震前后遥感影像的梯度图像。

优选地，通过如下公式，分别计算地震前后遥感影像的梯度图像：

$G(x,y)=\sqrt{{\left(f\right(x,y)-f(x+1,y+1\left)\right)}^2+{\left(f\right(x+1,y)-f(x,y+1\left)\right)}^2}$

$G(x^{\′},y^{\′})=\frac{G(x,y)-G_{\min }(x,y)}{G_{max}(x,y)-G_{min}(x,y)}$

其中，f(x,y)为源图像的单通道，G(x',y')为梯度图像。

步骤203，将地震前后遥感影像的梯度图像网格化。

较佳地，可以根据当地评估目标(建筑物)的平均长度L，以及获取的遥感影像空间分辨率R，选择网格大小。优选地，采用香浓采样原理，选取网格大小。在该实施例中，接近S＝0.5*L/R，当然根据实际情况可以进行调整。例如：当地建筑物平均长度为20米，遥感影像空间分辨率为0.5米，则网格尺度选择为20。

步骤204，分别计算网格化后的地震前后梯度图像的相似度特征图像。

优选地，根据如下公式，计算网格化后的地震前后梯度图像的相似度特征图像：

$S\left(W\right)=\frac{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)*\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}{\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_A\right(x,y)-G_{A-mean}(x,y\left)\right)}^2}*\sqrt{{\mathrm{\Σ}}_{x=1}^k{\mathrm{\Σ}}_{y=1}^l{\left(G_B\right(x,y)-G_{B-mean}(x,y\left)\right)}^2}}$

其中，G_A(x，y)表示T1时间(灾前)网格W的梯度特征图像。G_B(x，y)表示T2时间(灾后)网格W的梯度特征图像。

还有，G_A-mean(x，y)＝(∑∑G_A(x，y))/H，表示目标A区域中像元的梯度平均值。G_B-mean(x，y)＝(∑∑G_B(x，y))/H，表示目标B区域中像元的梯度平均值。H则为网格中像元的总个数，k和l分别为目标区域的宽度和高度。

步骤205，根据所述地震前后梯度图像的相似度图像的相似度特征，获取完全变化和严重变化网格数量。具体的实施过程包括：

通过ENVI软件选取适量的训练样本，并使用最大似然分类法对相似度图像进行分类。其中，在选取样本时，可以选取30至50个训练样本，另外参考建筑物按照EMS-1998标准，将Grades 3,4划分为完全变化类，将Grades1,2划分为严重变化类，Grade5划分为轻度变化类。

步骤206，根据评估倒毁率绘制当地地震灾害快速评估图。

优选地，利用Arcgis软件根据行政区划矢量数据及公式：

D_i＝A_i/B_i

计算每个行政区的倒毁率评估结果，上式中，A_i为行政区i对应的完全变化及严重变化网格数量，B_i为行政区i对应的完全变化及严重变化网格数量，D_i为行政i对应的估计倒毁率。也就是说，使用ArcGis软件，将行政区划中每个行政区按照Di值赋予不同的颜色，制成当地地震灾害快速评估图。

在本发明的另一方面，还提供了一种基于高分辨率遥感数据的地震损毁评估装置，如图3所示，所述基于高分辨率遥感数据的地震损毁评估装置包括依次连接的梯度图像获取单元301、相似度特征获取单元302以及评估单元303。其中，梯度图像获取单元301根据监测边缘特征变化，获得地震前后遥感影像的梯度图像。相似度特征获取单元302分别获得所述地震前后梯度图像的相似度图像，提取所述相似度图像的相似度特征。评估单元303根据评估倒毁率绘制当地地震灾害快速评估图。

需要说明的是，在本发明所述的基于高分辨率遥感数据的地震损毁评估装置的具体实施内容，在上面所述基于基于高分辨率遥感数据的地震损毁评估方法中已经详细说明了，故在此重复内容不再说明。

本发明根据上面所述的方法和装置，对2010年4月14日青海省玉树藏族自治州玉树县发生Ms7.1级地震进行了快速评估。图4显示的是根据本发明的一个具体实施例的地震灾害前后遥感影像数据图，其中左侧为地震前获取的ikonos影像，右侧震后获取的QuickBird影像；图5为图4基础上分别计算的梯度图像；图6是图5基础上的相似度特征图像；图7是基于图6的变化检测方法分类结果。图8是图7基础上的倒毁率评估图。

其中，如图4所示，获取当地地震前Ikonos遥感影像数据和地震后QuickBird遥感影像数据，并通过Arcgis软件进行空间配准。如图5所示，分别计算灾害前后影像的梯度图像，依据当地建筑物平均长度为20米，影像空间分辨率为0.6米，确定网格大小为17。还有，如图8所示，根据评估倒毁率绘制当地地震灾害快速评估图，按受损程度可以使用不同颜色展示。

综上所述，本发明提供的一种基于基于高分辨率遥感数据的地震损毁评估方法和装置，可在地震发生后，快速评估由地震引发的建城区建筑物倒毁情况，对选择启动应急预案的应急响应级别具有较高参考价值，为决策者指定灾民的营救和安置计划提供依据；进一步地，本发明通过先提取图像特征，后进行相似度对比的方式，进而可以兼容用于变化监测的灾害前后影像为异源传感器；从而，本发明具有广泛、重大的推广意义；最后，整个所述基于高分辨率遥感数据的地震损毁评估方法和装置紧凑，易于控制。

所属领域的普通技术人员应当理解：以上任何实施例的讨论仅为示例性的，并非旨在暗示本公开的范围(包括权利要求)被限于这些例子；在本发明的思路下，以上实施例或者不同实施例中的技术特征之间也可以进行组合，步骤可以以任意顺序实现，并存在如上所述的本发明的不同方面的许多其它变化，为了简明它们没有在细节中提供。

另外，为简化说明和讨论，并且为了不会使本发明难以理解，在所提供的附图中可以示出或可以不示出与集成电路(IC)芯片和其它部件的公知的电源/接地连接。此外，可以以框图的形式示出装置，以便避免使本发明难以理解，并且这也考虑了以下事实，即关于这些框图装置的实施方式的细节是高度取决于将要实施本发明的平台的(即，这些细节应当完全处于本领域技术人员的理解范围内)。在阐述了具体细节(例如，电路)以描述本发明的示例性实施例的情况下，对本领域技术人员来说显而易见的是，可以在没有这些具体细节的情况下或者这些具体细节有变化的情况下实施本发明。因此，这些描述应被认为是说明性的而不是限制性的。

尽管已经结合了本发明的具体实施例对本发明进行了描述，但是根据前面的描述，这些实施例的很多替换、修改和变型对本领域普通技术人员来说将是显而易见的。例如，其它存储器架构(例如，动态RAM(DRAM))可以使用所讨论的实施例。

本发明的实施例旨在涵盖落入所附权利要求的宽泛范围之内的所有这样的替换、修改和变型。因此，凡在本发明的精神和原则之内，所做的任何省略、修改、等同替换、改进等，均应包含在本发明的保护范围之内。