CN112347913B - 基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法及系统 - Google Patents

Abstract

本发明公开了基于快速Huynen‑Euler分解的受灾建筑物损坏等级估计方法及系统，所述方法包括：根据全极化合成孔径雷达的极化数据，分别获得目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]；根据目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]，分别计算目标灾前Huynen参数和目标灾后的Huynen参数；根据目标灾前Huynen参数和目标灾后的Huynen参数，利用快速Huynen‑Euler分解，分别计算目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2；根据目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2，计算目标跳跃角变化率参数Δv_n；根据跳跃角变化率参数Δν_n，估计受灾建筑物损坏等级DL。该方法解决了自然灾害突发时对城区受灾程度进行评估的过程耗时较大的问题，实现对受灾地区受灾程度快速准确的估计。

Description

基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法 及系统

技术领域

本发明涉及全极化合成孔径雷达数据的极化分解领域，特别涉及基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法及系统。

背景技术

随着微波技术的发展，能够有效获取目标的全极化散射特性的极化合成孔径雷达作为微波遥感的重要手段之一，被广泛应用于目标检测与识别、地物分类、地表参数反演等方面。近年来，由于自然灾害的频发及其程度的加剧，灾害监测和建筑区域损坏评估作为极化合成孔径雷达的新兴应用引起了广泛关注，因为其能够全天时、全天候、不受云雨气象条件的影响、经由遥感的方式迅速的获取大范围的有关受灾区域的数据。经由比较并分析灾前与灾后该区域的散射特性的变化，可以将其准确识别为损坏区域、水淹区域、未受损区域，等等。模型分解提取的二次散射功率是目前最常用的进行受灾区域状态判断的参数，因为其被认为是由建筑竖直的墙体和地表构成的二面角结构所产生的。故已有算法通过比较二次散射功率的变化来判断建筑区域受损程度，如文献1[S.-W.Chen,X.-S.Wang,andM.Sato,“Urban damage level mapping based on scattering mechanisminvestigation using fully polarimetric SAR data for the 3.11east japanearthquake,”IEEE Trans.Geosci.Remote Sens.,vol.54,no.12,pp.6919–6929,2016.]。然而，采用模型分解得到二次散射功率的过程相对耗时，因为在对这一参数进行求解的过程中需要经过去取向、同极化功率相对大小判断、根据判断结果选择求解方法、消除负功率影响等多个步骤，而灾害应用对参数提取的时效性要求很高——为了保证能够迅速的检测到灾害、对城区受灾程度进行快速评估以制定高效及时的救援措施，挽救人民群众生命安全和财产。

Huynen-Euler分解是Huynen提出的对目标极化散射矩阵进行的一种分解算法(文献2[J.R.Huynen,“Phenomenological theory of radar targets,”Ph.D.dissertation,Dept.Elect.Eng.,Math.Comput.Sci.,Tech.Univ.Delft,Delft,The Netherlands,1970.])，分解得到的五个欧拉参数能够有效反应目标不同方面的散射特性。其中跳跃角参数v与信号的反射次数相关，其本质上反应了与二次散射功率相同的信息。尽管理论上可行，Huynen提出的跳跃角参数在实际应用中存在角度模糊问题。采用快速Huynen-Euler分解可以通过解析的方式提取改进的跳跃角参数v_n，不仅避免了旧跳跃角参数角度模糊的问题，而且实现了速度远超于模型分解的对建筑区域的快速识别，从而达到对受灾地区快速检测和受灾程度迅速估计的目的。

发明内容

本发明的目的在于解决自然灾害突发时对城区受灾程度进行评估的过程耗时较大的问题。

为了实现上述目的，本发明提供了一种基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法，所述方法包括：

根据全极化合成孔径雷达的极化数据，分别获得目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]；

根据目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]，分别计算目标灾前Huynen参数和目标灾后的Huynen参数；

根据目标灾前Huynen参数和目标灾后的Huynen参数，利用快速Huynen-Euler分解，分别计算目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2；

根据目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2，计算目标跳跃角变化率参数Δv_n；

根据跳跃角变化率参数Δv_n，估计受灾建筑物损坏等级DL。

作为上述方法的一种改进，所述目标灾前的相干矩阵[T₁]为；

其中，T1₁₁、T1₁₂、T1₁₃、T1₂₁、T1₂₂、T1₂₃、T1₃₁、T1₃₂和T1₃₃为目标灾前的相干矩阵[T₁]的元素；

目标灾后的相干矩阵[T₂]为：

其中，T2₁₁、T2₁₂、T2₁₃、T2₂₁、T2₂₂、T2₂₃、T2₃₁、T2₃₂和T2₃₃为目标灾后的相干矩阵[T₂]。

作为上述方法的一种改进，所述根据目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]，分别计算目标灾前Huynen参数和目标灾后的Huynen参数；具体包括：

将目标灾前和灾后的Huynen参数分别记为A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁、G₁和H₁，其计算公式为：

C₁＝Re{T1₁₂}

D₁＝-Im{T1₁₂}

E₁＝Re{T1₂₃}

F₁＝Im{T1₂₃}

G₁＝Im{T1₁₃}

H₁＝Re{T1₁₃}

将目标灾后的Huynen参数分别记为A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂和H₂，其计算公式为：

C₂＝Re{T2₁₂}

D₂＝-Im{T2₁₂}

E₂＝Re{T2₂₃}

F₂＝Im{T2₂₃}

G₂＝Im{T2₁₃}

H₂＝Re{T2₁₃}

其中，Re{·}表示取实部，Im{·}表示取虚部。

作为上述方法的一种改进，所述根据目标灾前Huynen参数和目标灾后的Huynen参数，利用快速Huynen-Euler分解，分别计算目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2；具体包括：

计算目标灾前的跳跃角参数v_n1的公式为：

计算目标灾后的跳跃角参数v_n2的公式为：

其中，v_n1∈[0，45°]，v_n2∈[0，45°]。

作为上述方法的一种改进，所述根据目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2，计算目标跳跃角变化率参数Δv_n；具体包括：

计算跳跃角变化率参数Δv_n的公式为：

其中，Δv_n∈[0,1]。

作为上述方法的一种改进，所述根据跳跃角变化率参数Δv_n，估计受灾建筑物损坏等级DL，具体包括：

受灾建筑物损坏等级DL为：

其中，损坏等级DL∈[0,1]。

本发明还提出了一种基于快速Huynen-Euler分解的受灾建筑物损坏等级估计系统，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的方法。

本发明的优势在于：

本发明提出的方法和系统解决了自然灾害突发时对城区受灾程度进行评估的过程耗时较大的问题，既避免了旧跳跃角参数的角度模糊问题，通过快速分解方法实现对建筑区域的快速识别，解决了自然灾害突发时对城区受灾程度进行评估的过程耗时较大的问题，能够得到对受灾地区受灾程度快速准确的估计。

附图说明

图1为本发明的基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法的流程图；

图2(a)为受灾地区日本石卷市的灾前的光学图像；

图2(b)为受灾地区日本石卷市的灾后的光学图像；

图3为建筑物损坏等级估计结果的示意图；

图4为利用模型分解方法得到的建筑物损坏等级估计结果。

具体实施方式

下面结合附图和具体实施例对本发明进行详细的说明。

Huynen-Euler分解提取的跳跃角参数v与信号的反射次数相关，理论上能够对城区进行识别，但在实际应用中其存在角度模糊的问题。采用快速Huynen-Euler分解可以通过解析的方式提取改进的跳跃角参数v_n，既能够有效避免旧跳跃角参数角度模糊的问题，而且求解速度极快。本发明通过提出改进的跳跃角参数v_n并提出求解v_n的快速Huynen-Euler分解方法，以解决自然灾害突发时对城区受灾程度进行评估的过程耗时较大的问题，从而得到一种既能避免旧跳跃角参数的角度模糊，又能对建筑区域快速识别，进而实现对受灾地区受灾程度快速准确估计的方法。

如图1所示，本发明的实施例1提出了的基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法，包括以下步骤：

步骤1)根据灾前和灾后的全极化合成孔径雷达的极化数据，分别获得目标灾前和灾后的相干矩阵[T₁]和[T₂]，分别计算灾前和灾后的Huynen参数A₀，B₀，B，C，D，E，F，G和H，并分别记为A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁，G₁，H₁(灾前)和A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂，H₂(灾后)；

步骤2)根据步骤1)得到的灾前和灾后的Huynen参数A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁，G₁，H₁(灾前)和A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂，H₂(灾后)，通过快速Huynen-Euler算法，分别计算目标灾前和灾后的改进的跳跃角v_n1和v_n2；

步骤3)根据步骤2)得到的目标灾前和灾后的改进的跳跃角v_n1和v_n2，计算目标的跳跃角变化率Δv_n；

步骤4)根据步骤3)得到的目标的跳跃角变化率Δv_n，估计受灾建筑物损坏等级DL；

下面通过实例进一步描述。

在一个实例中，受灾地区日本石卷市的灾前与灾后的光学图像如图2(a)和图2(b)所示，该图像由Google Earth分别在2010年6月25日和2011年4月6日获取。在2011年3月11日，发生了9.0级的日本东北地方太平洋近海地震，并引发最大爬高40.1米的海啸，是日本史上最大的天然灾难。地震造成至少15893人死亡、2553人失踪、伤者6152人，遭受破坏的房屋1292417栋，为日本二战后伤亡最惨重的自然灾害。由图2可以看出，灾前灾后的沿海城区发生了明显变化，远离海边的地区则没有明显的变化。这是因为沿海城区受灾极其严重，多数房屋倒塌且残骸被随之而来的海啸冲走，而内陆地区受这两种灾害的影响则较小。该实施例中采用的全极化合成孔径雷达数据是由JAXA的星载ALOS PALSAR系统分别在2009年4月2日和2011年4月8日获取。

在步骤1)中，根据灾前和灾后的全极化合成孔径雷达的极化数据，分别获得目标灾前和灾后的相干矩阵[T₁]和[T₂]，分别计算灾前和灾后的Huynen参数A₀，B₀，B，C，D，E，F，G和H，并分别记为A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁，G₁，H₁(灾前)和A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂，H₂(灾后)；其中，

目标的相干矩阵[T]表示如下：

将目标灾前灾后的相干矩阵分别记为[T₁]和[T₂]。

Huynen参数A₀，B₀，B，C，D，E，F，G和H的计算公式如下：

其中，Huynen参数A₀，B₀，B，C，D，E，F，G和H都是实数。将目标的相干矩阵[T]采用实数参数表示具有清晰明了的优点，此外还能够方便的应用快速惠能欧拉分解。利用Huynen参数进行目标散射信息的表示具有简洁直观的优点。上述两种表示是等价的。将目标灾前和灾后的Huynen参数分别记为A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁，G₁，H₁(灾前)和A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂，H₂(灾后)。

在步骤2)中，根据步骤1)得到的灾前和灾后的Huynen参数A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁，G₁，H₁(灾前)和A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂，H₂(灾后)，通过快速Huynen-Euler算法，分别计算目标灾前和灾后的改进的跳跃角v_n1和v_n2；其中，

快速Huynen-Euler算法是本发明提出的一种能够快速实现Huynen-Euler分解并提取相应的Huynen-Euler参数的算法。Huynen-Euler分解是Huynen基于包含目标的全部极化信息的散射矩阵[S]提出的一种分解方法，表示如下：

其中，

各参数都具备特定的与目标散射特性相关的含义：m²表示目标的最大雷达散射截面积；是目标的取向角，与目标绕雷达视线取向有关；τ是螺旋角，与目标对称性有关(人造目标τ＝0，自然地物τ＝45°)；γ是极化度角，与目标极化敏感度有关(线性目标γ＝0，如偶极子；球或平面目标γ＝45°)；v是跳跃角，与多次散射有关(单次散射v＝0，二次散射v＝45°)。由此可以看出v能够有效识别城区等二次散射目标，与二次散射功率参数具有相似的性质。Huynen对各个参数的取值范围给出了描述：

已有的算法通过对[S]进行伪特征分解或对进行特征分解来求解五个Huynen-Euler参数(其中上标干表示取共轭转置)，但这一过程较为耗时。本发明通过将散射矩阵[S]转化为与其等价的相干矩阵[T]，并进一步利用Huynen参数对目标全部的散射信息进行表示，推导出五个Huynen-Euler参数的解析解形式，实现快速Huynen-Euler分解并得到五个Huynen-Euler参数。其中，将散射矩阵[S]转化为与其等价的相干矩阵[T]利用如下变换：

然后再根据式(2)便可以获得目标的Huynen参数。根据式(2)、(3)和(6)，我们可以将Huynen参数用五个Huynen-Euler参数表示，如下：

利用式(7)可以推导出Huynen-Euler参数的解析解如下：

需要说明的是，原Huynen-Euler分解提取的跳跃角参数v对于建筑物的识别存在角度模糊问题，例如，对于两个散射矩阵为如下所示对角阵的目标[S₁]和[S₂]

原Huynen-Euler分解提取的跳跃角参数v₁＝-45°，v₂＝45°，即将[S₁]和[S₂]识别为不同的目标，而实际上可以推导出

[S₁]＝e^j180°[S₂] (14)即二者实际上只相差一个绝对相位。考虑到跳跃角参数v反映的是与目标二次散射有关的信息，而标准二次散射最重要的特点是同极化通道的相位相反，因此本发明认为v的取值只需要反应同极化通道的相对相位即可，即，采用|v|对目标的二次散射信息进行描述。综上，本发明对于旧的跳跃角参数进行了改进，通过令

v_n＝|v| (15)

本发明采用v_n对目标的二次散射信息进行描述，且有

因此，本发明根据式(8)-(11)和(16)，实现对目标的快速Huynen-Euler分解，并利用式(16)分别计算得到目标灾前和灾后的改进的跳跃角参数v_n1和v_n2。

在步骤3)中，根据步骤2)得到的目标灾前和灾后的改进的跳跃角v_n1和v_n2，计算目标的跳跃角变化率Δv_n；其中，

Δv_n反映了目标跳跃角v_n在受灾前后的变化情况，能够反映二面散射减少的相对量，从而反映出受灾城区的受损程度。

在步骤4)中，根据步骤3)得到的目标的跳跃角变化率Δv_n，估计受灾建筑物损坏等级DL；其中，

受灾建筑物损坏等级DL的估计方法如下：

其中，Δv_n＜0.2时候我们认为目标受灾程度为0，这是因为尽管理论上v_n的变化只与二次散射相关，但实际上由于获取数据具有一定的时间间隔，目标的散射特性极有可能发生变化，噪声等因素的影响也会导致v_n一定程度的变化，因此为了避免随机因素和时间因素的影响，本发明采取了这一措施。同理，我们将Δv_n≤1也作为基本的约束条件，以保证获得的建筑物受灾程度结果有意义。

在本实施例中，利用快速Huynen-Euler分解估计得到的受灾建筑物损坏等级如图3所示，可以看出，其能够准确表示严重受损的沿海城区的与非沿海区域的无损特点，而且不受海洋、山区和植被区域的干扰，即这些非城区得到的受灾程度估计结果都近乎为0。已有的利用模型分解得到的二次散射功率进行城区受损程度估计的算法的估计结果如图4所示，可以看出，其结果与本发明的结果很相近，二者的相干性高达0.7471，但很显然，已有算法不能避免海洋、山区和植被区域的干扰，其在这些地区都有较高的DL值。

除此之外，本发明的快速Huynen-Euler分解另一巨大优势在于求解的高效性。通过对本组实施例数据的计算时间记录(本实施例数据为复数矩阵1248×18432)，可以得到：本发明提出的算法得到计算DL值所需参数v_n的时间是0.2458秒，而已有算法通过模型分解得到计算DL值所需参数二次散射功率的时间是4.9178，差距不可谓之不大。因而，本发明——快速Huynen-Euler分解——的求解的高效性对于灾害检测、受灾城区受损程度估计、救援措施的迅速制定与执行有着重要意义。

本发明的实施例2提出了一种基于快速Huynen-Euler分解的受灾建筑物损坏等级估计系统，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现实施例1的方法。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (4)

1.一种基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法，所述方法包括：

根据全极化合成孔径雷达的极化数据，分别获得目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]；

根据目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]，分别计算目标灾前Huynen参数和目标灾后的Huynen参数；

根据目标灾前Huynen参数和目标灾后的Huynen参数，利用快速Huynen-Euler分解，分别计算目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2；

根据目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2，计算目标跳跃角变化率参数Δν_n；

根据跳跃角变化率参数Δv_n，估计受灾建筑物损坏等级DL；

所述目标灾前的相干矩阵[T₁]为；

其中，T1₁₁、T1₁₂、T1₁₃、T1₂₁、T1₂₂、T1₂₃、T1₃₁、T1₃₂和T1₃₃为目标灾前的相干矩阵[T₁]的元素；

目标灾后的相干矩阵[T₂]为：

其中，T2₁₁、T2₁₂、T2₁₃、T2₂₁、T2₂₂、T2₂₃、T2₃₁、T2₃₂和T2₃₃为目标灾后的相干矩阵[T₂]；

所述根据目标灾前的相干矩阵[T₁]和目标灾后的相干矩阵[T₂]，分别计算目标灾前Huynen参数和目标灾后的Huynen参数；具体包括：

将目标灾前和灾后的Huynen参数分别记为A₀₁，B₀₁，B₁，C₁，D₁，E₁，F₁、G₁和H₁，其计算公式为：

C₁＝Re{T1₁₂}

D₁＝-Im{T1₁₂}

E₁＝Re{T1₂₃}

F₁＝Im{T1₂₃}

G₁＝Im{T1₁₃}

H₁＝Re{T1₁₃}

将目标灾后的Huynen参数分别记为A₀₂，B₀₂，B₂，C₂，D₂，E₂，F₂，G₂和H₂，其计算公式为：

C₂＝Re{T2₁₂}

D₂＝-Im{T2₁₂}

E₂＝Re{T2₂₃}

F₂＝Im{T2₂₃}

G₂＝Im{T2₁₃}

H₂＝Re{T2₁₃}

其中，Re{·}表示取实部，Im{·}表示取虚部；

所述根据目标灾前Huynen参数和目标灾后的Huynen参数，利用快速Huynen-Euler分解，分别计算目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2；具体包括：

计算目标灾前的跳跃角参数v_n1的公式为：

计算目标灾后的跳跃角参数ν_n2的公式为：

其中，ν_n1∈[0,45°]，ν_n2∈[0,45°]。

2.根据权利要求1所述的基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法，其特征在于，所述根据目标灾前的跳跃角参数v_n1和目标灾后的跳跃角参数v_n2，计算目标跳跃角变化率参数Δv_n；具体包括：

计算跳跃角变化率参数Δv_n的公式为：

其中，Δv_n∈[0,1]。

3.根据权利要求2所述的基于快速Huynen-Euler分解的受灾建筑物损坏等级估计方法，其特征在于，所述根据跳跃角变化率参数Δv_n，估计受灾建筑物损坏等级DL，具体包括：

受灾建筑物损坏等级DL为：

其中，损坏等级DL∈[0,1]。

4.一种基于快速Huynen-Euler分解的受灾建筑物损坏等级估计系统，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-3任一项所述的方法。