一种基于立体影像和激光测高数据的建筑物倒塌评估方法

1.本发明属于图像处理的技术领域，具体涉及一种基于立体影像和激光测高数据的建筑物倒塌评估方法。


背景技术：

2.传统的灾害评估能够获取受损房屋的位置、面积等二维信息，但是难以评估出房屋的三维信息的变化。立体卫星具有幅度宽、响应周期短等特点，能够计算出房屋三维信息的变化，十分契合地震灾区评估工作要求的快响应、广范围。使用hrsi进行地震灾害评估需要通过地面控制点提升hrsi的垂直及平面定位精度,具体方法为将地面控制点与hrsi进行配准并构建控制点误差方程，hrsi影像间连接点提取并构建连接点误差方程，求解误差方程得到对立体影像的补偿系数。卫星立体影像的定位精度严重依赖于立体影像本身的地面分辨率以及地面控制点的布设，然而由于震前震后的环境差异较大，地面控制点的快速获取十分困难，导致震前震后立体影像之间的平面、垂直定位精度存在不一致性，同时房屋受损评估时面临的房屋角点震前震后受损导致获取困难等问题，进而导致房屋评估结果存在不准确性。


技术实现要素：

3.针对使用hrsi进行受损房屋的评估会遇到诸如平面定位不一致、高程定位精度不足、灾后地面控制点布设困难、受损房屋的像方难以选择等问题，本发明提出了一种基于高分七号(gf-7)震前震后立体影像及激光测高数据的建筑物倒塌三维精细化评估方法，首先针对震前震后立体影像平面定位不一致、垂直定位精度不足，使用gf-7和icesat1激光测高数据作为高程控制、震前震后四重连接点作为约束进行区域网平差；其次本研究针对震前震后立体影像的辐射差异以及房屋角点受损前后的巨大差异性，提出了一种由粗到精的亚像元房屋角点提取方法，先由单应性矩阵、核线约束进行粗配准，再由最小二乘模板匹配进行精配准，获取震前震后立体影像上房屋角点的像方；最后依据震前震后立体影像计算房屋角点高度变化以及ems分成的三种类别进行建筑物倒塌三维精细化评估。
4.本发明可通过以下技术方案实现：
5.一种基于立体影像和激光测高数据的建筑物倒塌评估方法，包括以下步骤：
6.步骤一、采用由粗到精匹配方法对震前、震后的多张立体影像进行匹配并结合仿射模型构建连接点误差方程，同时，对内部激光测高数据即足印影像与立体影像进行匹配以实现激光点的提取并结合仿射模型构建内部激光点误差方程，以及对外部激光测高数据与立体影像进行匹配并结合仿射模型构建外部激光点误差方程，共同构建基于有理函数成像模型的平差模型；
7.步骤二、对平差模型进行求解，完成区域网平差并获取每张立体影像对应的补偿系数，从而获取每张立体影像对应的平差后有理函数成像模型；
8.步骤三、采用由粗到精定位方法实现建筑物角点的自动量测，获取各个角点对应
的像方坐标；
9.步骤四、将步骤三中的像方坐标分别输入震前、震后的立体影像的平差后有理函数成像模型，计算出建筑物的高度变化，实现建筑物倒塌程度的评估。
10.进一步，所述步骤一中构建连接点误差方程的方法包括以下步骤：
11.步骤ⅰ、选取多张立体影像中的一张为主影像，其余为副影像，使用sift算法对主影像进行特征点提取；
12.步骤ⅱ、计算主影像上每个特征点对应的经度及纬度，以此为目标特征点，结合主影像的最大及最小高程值，输入其他副影像的有理函数成像模型，获取包含目标特征点的影像块，再以主影像上的目标特征点为模板，采用最小二乘匹配算法，计算出目标特征点在其他副影像上的精细坐标；
13.步骤ⅲ、以主影像的目标特征点和其他副影像对应的匹配点作为连接点，构建连接点误差方程。
14.进一步，所述步骤ⅰ中使用sift算法对主影像进行特征点提取时，先对主影像进行分块，再使用sift对每一图像块进行特征点提取，然后使用anms算法均匀化对每个影像块中的特征点进行均匀化，保留响应值大于20％的特征点。
15.进一步，所述步骤一中构建内部激光点误差方程的方法包括以下步骤：
16.步骤
ⅰ
、使用sift对足印影像进行特征点提取；
17.步骤
ⅱ
、计算足印影像上每个特征点对应的经度及纬度，以此为目标特征点，结合足印影像的最大及最小高程值，输入各个立体影像的有理函数成像模型，获取包含目标特征点的影像块，再以足印影像的目标特征点为模板，采用最小二乘匹配算法，计算出目标特征点在各个立体影像上的精细坐标；
18.步骤
ⅲ
、由足印影像和各个立体影像之间的匹配点对计算出由足印影像到各个立体影像的单应性矩阵，输入激光点在足印影像上的像方坐标到该矩阵，获得激光点在立体影像上的像方坐标，构建内部激光点误差方程。
19.进一步，将外部激光测高数据投影至各个立体影像上，再利用根据权利要求2获得的连接点构建对应的单应性矩阵对激光测高数据的粗配准，然后对所有立体影像分别利用最小二乘模板匹配进行激光测高数据的精配准，以外部激光测高数据和各个立体影像对应的匹配点，构建外部激光点误差方程的方法。
20.进一步，所述步骤一中构建基于有理函数成像模型的平差模型如下：
[0021][0022]
其中，第一类方程为连接点误差方程，第二类为虚拟控制点误差方程，第三类为内部和外部共同组成的激光点误差方程，v
tp
，v
vcp
，v
las
分别代表连接点、虚拟控制点以及激光点的误差项，a
tp
，a
vcp
，a
las
分别代表连接点、虚拟控制点以及激光点的对仿射变换补偿模型中的待定求解系数的偏导项，x
aff
为待定求解的补偿系数，b
tp
，b
vcp
，b
las
分别代表连接点、虚拟控制点以及激光点的对点的物方的偏导项，x
tp
，x
vcp
，x
las
分别代表连接点、虚拟控制点以及激光点的物方项，l
tp
，l
vcp
，l
las
分别代表连接点、虚拟控制点以及激光点的残差项，p
tp
，p
vcp
，p
las
分别代表连接点、虚拟控制点以及激光点的权重项。
[0023]
进一步，所述步骤三中实现建筑物角点的自动量测如下：首先在主影像上通过目视解译选取建筑物角点的精细像方坐标，结合根据权利要求2获得的连接点构建对应的单应性矩阵，实现建筑物角点在其他副影像上的粗定位，然后以主影像上建筑物角点的精细像方坐标为模板，采用最小二乘匹配算法，保留相关系数大于0.8的点，自动计算出建筑物角点在其他副影像上的精细像方坐标。
[0024]
本发明有益的技术效果在于：
[0025]
使用的是激光测高数据并统一了震前震后立体影像的平面定位精度，避免使用耗时耗力且不利于地震灾害评估工作要求的时效性地面控制点，仅考虑垂直方向的定位精度，且完成了激光测高数据的自动配准。有控情况下，相较于震前立体影像单独定位情况下，统一震前、震后定位精度后，高程精度得到明显提升，提升为52％，垂直精度rmse从2m提升至0.98m，影像内精度rmse从0.63提升至0.33；同时，无控情况下，震前的定位精度在统一震前、震后，定位精度提升为82％，垂直精度rmse从10.28m提升至2.3m，影像内精度rmse从0.57变为0.65。
[0026]
对于建筑物角点测量，以目视解译为基础，使用单应性矩阵及最小二乘模板匹配的由粗到精的自动测量方法，能够完成对震前震后多个建筑物角点的自动测量，测量效果良好，大大节省了人力，避免不同个体量测精度参差不齐等问题。对灾区选取四处受损严重的区域进行评估，正确率达到了a区域正确率93.6％，b区域47％，c区域97.3％，d区域100％，满足灾区房屋受损情况的评估的精度要求。
附图说明
[0027]
图1为本发明的总体流程示意图；
[0028]
图2为本发明的克服不同时相辐射差异的同名点提取方法流程示意图；
[0029]
图3为本发明的足印影像和立体影像的匹配流程示意图；
[0030]
图4为本发明的房屋角点自动获取的流程示意图。
[0031]
图5为本发明的针对单片区域的评估结果示意图。
[0032]
图6为本发明的使用目视解译选取了a、b、c、d四处产生变化的房屋区域的示意图。
具体实施方式
[0033]
下面结合附图及较佳实施例详细说明本发明的具体实施方式。
[0034]
如图1所示，本发明提供了一种基于立体影像和激光测高数据的建筑物倒塌评估方法，基于sift和最小二乘匹配可以克服不同时相数据即震前震后多张立体影像的辐射差异以提取同名点，结合外部激光测高辅助数据通过区域网平差消除震前、震后几何定位不一致性并提升三维空间垂直定位精度且消除相对平面定位误差后，分别对震前、震后立体影像前方交会实现对房屋角点高度变化的评估，再根据灾区房屋的建筑标准来评估灾区每一栋建筑的受损程度，具体如下：
[0035]
步骤一、采用由粗到精匹配方法对震前、震后的多张立体影像进行匹配并构建连接点误差方程，如震前、震后的各两张立体影像，同时，对内部激光测高数据即足印影像与立体影像进行匹配以实现激光点的提取并构建内部激光点误差方程，以及对外部激光测高数据与立体影像进行匹配并构建外部激光点误差方程。
[0036]
震前与震后的立体影像数据是不同时相的数据，采集于不同时间，存在基准不统一的问题。在不同时相情况下，卫星的轨道定位、姿态、星下点定位等定位参数都不同，导致不同时相立体影像的成像范围、水平、高程定位精度存在很大的差异性。本发明基于震前、震后投影光线的数量，分别针对连接点、激光点、虚拟控制点构建了不同的误差方程，形成更严格的约束，以期求得更精准的定位。
[0037]
1、震前震后立体影像的同名点提取，以构建连接点及其连接点误差方程
[0038]
相关研究表明，连接点的分布对立体影像之间的相对定位误差的消除有着直接影响，因此，我们以震前震后的多张立体影像如各两张为研究数据，选取其中一张为主影像，其余为副影像，如图2所示，首先对主影像进行分块，使用sift算法对每一影像块进行特征点提取，这样可以获得在整张影像上分布较为均匀的特征点，再使用anms算法对单个影像块中的特征点进行均匀化，保留响应值大于20％的特征点；
[0039]
然后，在立体影像之间采用由粗到精的策略来获取已有特征点在其他副影像上的同名点，即对主影像的每个特征点，以此作为目标特征点，通过主影像的地面分辨率及其行列号，可以粗略计算出其经度及纬度，结合主影像的最大及最小高程值(可由rpcs计算得出)，输入其他副影像的有理函数成像模型，就可以在各个副影像上可以计算出包含目标特征点的影像块，完成粗定位；
[0040]
最后，以主影像上的特征点为模板，采用最小二乘匹配算法，可以计算出目标特征点在其他副影像上的精细坐标，从而完成同名点的提取。
[0041]
由主影像的每个特征点和各个副影像与之匹配的特征点组成连接点，结合仿射变换模型构建连接点误差方程。
[0042]
2、内部激光测高数据及外部激光测高数据的配准
[0043]
高分七号搭载了全波形激光测高仪及其足印相机，用于获取高精度的激光控制点即内部激光测高数据。足印影像的特点在于分辨率低(3.2m)，如何将足印影像和立体影像进行配准是能否成功提取激光点的关键，如图3所示，本发明采用寻找足印影像和立体影像之间的转换关系的方法来进行足印影像和立体影像的配准。
[0044]
首先，使用sift算法对足印影像进行特征点提取，以此作为目标特征点，再采用和连接点提取时一致的粗定位方法，计算出各个立体影像上对应包含目标特征点的影像块，
[0045]
然后，对粗定位出的影像块进行降采样至与足印影像分辨率一致，以在足印影像上的特征点为模板，采用最小二乘匹配算法，得到其在降采样之后的影像块上的精细坐标。
[0046]
对于gf-7本身的激光点的筛选而言，选取三级产品中ecp_flag为1和2的部分作为质量较高的激光点。由足印影像和立体影像之间的匹配点对可以计算出由足印影像到立体影像的单应性矩阵，输入激光点在足印影像上的像方到单应性矩阵，可以得到激光点在立体影像上的像方，再结合仿射变换模型构建内部激光点误差方程。
[0047]
在实际的建筑物受损评估中，待评估区域可能存在激光点分布不均匀的问题。本发明还使用外部激光测高数据与立体影像进行匹配，具体的步骤为：
[0048]
对于外部激光测高数据的筛选，选取坡度小于2的作为质量较高的激光点。然后将外部激光测高数据投影至单张立体影像上得到像方坐标i1，根据上述方法中同名点的提取结果构建单应性矩阵，使用该单应性矩阵进行激光测高数据的粗匹配(将i1输入至单应性矩阵得到激光测高数据在其他立体影像上的粗略像方i2,i
3,
i4)，最后对所有立体影像分别
利用最小二乘模板匹配进行激光测高数据的精细化匹配，并以匹配点对结合仿射变换模型构建外部激光点误差方程。
[0049]
步骤二、以上述连接点误差方程、内部和外部激光点误差方程为基础，构建基于立体影像的有理函数成像模型的平差模型，对该平差模型进行求解，完成区域网平差并获取每张立体影像对应的补偿系数，从而获取每张立体影像对应的平差后有理函数成像模型。
[0050]
立体影像对/组之间本身存在相对定位误差，往往基于有理函数成像模型和已提取的同名点对/组构建误差模型，对误差模型进行求解，最终消除立体影像之间的相对定位误差。
[0051]
有理函数成像模型是一种对严格成像模型的数学拟合形式，为了消除卫星推扫过程、不精确的姿态以及轨道定位带来的误差，研究往往采用仿射变换模型去对立体卫星影像的行和列进行补偿以修正这些误差。基于有理函数成像模型的平差模型构建如下：
[0052]
v＝ax1 bx
2-l
ꢀꢀꢀ
(1)
[0053]
其中，v为误差余项；a为对仿射变换补偿模型中的待定求解系数的偏导，x1为与之对应的待定系数补偿量；b为对点的物方的偏导，x2为与之对应的物方补偿量；l为加入补偿模型后与实际的像方残差；
[0054]
假设提上文中选取的立体影像分别为建筑物震前的后视立体影像、前视立体影像与震后的后视立体影像、前视立体影像，其同名点分别记为与震后的后视立体影像、前视立体影像，其同名点分别记为补偿系数分别记为则此时的误差模型由原来的二重增至四重，误差方程中对应的系数矩阵a、b及残差矩阵l变为：
[0055][0056]
其中，表示有理函数成像模型补偿形式对仿射补偿系数的偏导，表示有理函数成像模型补偿形式对同名点物方的偏导，x
rfm
表示由立体影像的有理函数成像模型联合物方计算得出的像方。
[0057]
上述公式(1)是平差模型的一般形式，在实际研究中，往往需要针对不同的点作出相应的调整，如唐新明提出的平面、高程控制分离，针对激光点测高点，由于其极高的测距精度但同时平面精度较差，往往将高程值视为精确值，水平方向(经度、纬度)视为存在误差的值，即上述平差模型中不需要对高程进行偏导，而只需要对水平方向进行偏导。
[0058]
本发明的平差模型主要包含以下三种类型：
[0059][0060]
其中，第一类方程为连接点误差方程，第二类为虚拟控制点误差方程，第三类为激光点误差方程。v
tp
，v
vcp
，v
las
分别代表连接点、虚拟控制点以及激光点的误差项，a
tp
，a
vcp
，alas
分别代表连接点、虚拟控制点以及激光点的对仿射变换补偿模型中的待定求解系数的偏导项，x
aff
为待定求解的补偿系数，b
tp
，b
vcp
，b
las
分别代表连接点、虚拟控制点以及激光点的对点的物方的偏导项，x
tp
，x
vcp
，x
las
分别代表连接点、虚拟控制点以及激光点的物方项，l
tp
，l
vcp
，l
las
分别代表连接点、虚拟控制点以及激光点的残差项，p
tp
，p
vcp
，p
las
分别代表连接点、虚拟控制点以及激光点的权重项。针对不同类型的点，x中待求解的未知量也不同，例如，在连接点误差方程中，经度、纬度、高程三个未知量都需要求解；但是在激光点误差中，由于激光点较高的垂直精度但水平精度较低，经度、纬度都需要求解，高程视为固定值。此外，由于本发明不需要高精度的绝对水平控制，只使用连接点和激光点会导致法方程病态，所以需要加入了虚拟控制点约束了平面自由度，将单张立体影像的平均高程(由有理函数定位参数获得)作为虚拟控制点的高程并利用立体影像左上角的物方以及立体影像的地面分辨率获得虚拟控制点的经度以及纬度共同确定虚拟控制点的物方，将立体影像进行均匀分块如九块，每块取一个虚拟控制点，利用公式(1)构建虚拟控制点的误差方程，来加速了方程求解的收敛速度。针对不同类型的点，x中待求解的未知量也不同，例如，在连接点误差方程中，经度、纬度、高程三个未知量都需要求解；但是在激光点误差方程中，由于激光点较高的垂直精度但水平精度较低，经度、纬度都需要求解，高程视为固定值；在虚拟控制点误差方程中，将物方视为固定值，由于其是假设的平面约束，权重设为连接点权重的1/100。
[0061]
步骤三、采用由粗到精定位方法实现建筑物角点的自动量测，获取各个角点对应的像方坐标。
[0062]
震前与震后的建筑物如房屋之间存在较大的差异性，如变形、掀瓦和倒塌等，此外，由于拍摄时间不同，数据之间的辐射差异比较大。选取建筑物变化处的像方，采用消除立体影像之间相对定位误差后的定位模型，得到建筑物的物方，再将得到的物方反投影至另一对立体影像上，然后通过目视解译，就可以得到亚像素精度的建筑物像方。
[0063]
立体影像的成像范围较大，一般需要通过目视解译去选择出精细的建筑物像方，这种方式需要在多张影像上进行比对，耗时耗力，如图4所示，因此本发明提出一种由粗到精的自动定位方法，用于房屋像方的精细量测。
[0064]
首先，在主影像上通过目视解译选取房屋屋顶的角点的精细像方坐标，结合上文提取的连接点所构建的单应性矩阵即核线约束，可以实现房屋角点在其他副影像上的粗定位，再以主影像上的房屋角点的精细像方坐标为模板，结合粗定位的影像块，输入最小二乘模板匹配算法，保留相关系数大于0.8的点，就可以自动计算出房屋角点在其他副影像上的精细像方坐标。
[0065]
值得指出的是，若震前、震后匹配成功了，说明当前房屋在震前、震后没有遭到破坏；若是震前、震后匹配没有匹配成功，则说明震前、震后房屋变化较大，无法成功匹配，此时可以直接使用连接点构建的单应性矩阵即核线约束的变换结果作为倒塌房屋的位置，即获得房屋的包括经纬度信息的像方坐标。
[0066]
步骤四、将步骤三中的像方坐标输入立体影像的平差后有理函数成像模型，计算出建筑物的高度变化，实现建筑物倒塌程度的评估。
[0067]
european macroseismic scale 1998(ems98)将对房屋受损程度的评估划分的五个等级，分别为轻微受损，中度受损，严重受损，十分严重受损和倒塌。本研究结合ems98并挖掘遥感影像的最大评估效果，将第二至第四等级合并为一个等级，即轻微受损、中度受
损、倒塌三个等级，如下表所示。
[0068][0069]
(摘自ems98)
[0070]
使用立体影像计算出的高度的基准为wgs84，并不是房屋的真实高度，本研究以此为基础，结合灾区房屋的高度变化的比例，推算出房屋的受损程度，有受损程度评估公式(6)，dj表示单栋房屋的四个角点的受损程度，其中h
pre
,h
post
分别为震前、震后房屋的量测高度,h
ground
为房屋附近的地表高度，式(7)表示倒塌房屋楼层数的具体计算方法，其中c表示当地房屋的建筑标准高度。
[0071]
为了充分验证本发明提出的联合激光测高数据与震前、震后的立体影像的灾区房屋三维变化检测与房屋倒塌评估方法，实验选取了拍摄于墨西哥的瓦哈卡州的高分七号立体影像。该区域于2020年6月23日发生了7.4级的地震，造成了一定的房屋破坏和人员伤亡，本发明基于拍摄于灾前和灾后的两组高分七号立体影像(4张立体影像)对上文提出的联合激光测高数据与震前、震后的立体影像的处理方法以及对房屋受损程度的评估方法进行验证，具体结果如图5所示。
[0072]
为了验证本发明提出的震前震后联合平差的方法，设计了四种影像组合，分别为：
①
震前单独定位
②
震后单独定位
③
震前震后立体影像联合平差后使用震前影像单独定位
④
震前震后立体影像联合平差后使用震后影像单独定位。如表1所示，有控情况下即使用激光测高数据进行评估的情况下，相较于震前立体影像单独定位情况下，统一震前、震后定位精度后，高程精度得到明显提升，提升为52％，垂直精度rmse从2m提升至0.98m；同时，无控情况下，震前的定位精度在统一震前、震后，定位精度提升为82％，垂直精度rmse从10.28m提升至2.3m。
[0073]
表1
[0074]
[0075][0076]
本发明结合灾区报道情况及立体影像范围内房屋受损目视解译情况，选取了a、b、c、d四处产生变化的房屋区域，分布如图6所示，并利用前文提出的房屋受损程度评估的方法进行评估，评估结果如表2所示，其中a区域正确率达到了93.6％，b区域正确率达到了47％，c区域正确率达到了97.3％，d区域正确率达到了100％，满足灾区房屋受损情况的评估的精度要求。
[0077]
表2给出了受损/变化房屋的结果。表2给出了检测到的房屋变化的编号，由震前、震后立体影像联合平差后，由震前、震后立体影像定位得到房屋角点的高度及其变化，实际高度，实际层数，评估变坏层数以及评估的受损程度(完好，部分倒塌，完全倒塌，新建)。
[0078]
表2中，a、b处震前、震后四个角点均出现较大的高度变化，并且减少的高差几乎相同，所以被检测为完全倒塌；c、d处震后相较于震前，四个角点均出现了明显的高度增量，用于辅助对照倒塌房屋的检测，且也被检测为新建房屋，说明本发明所提出的方法能够检测到房屋明显的倒塌或新建的高度变化。
[0079]
表2
[0080][0081]
以上所述，以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。