0 引言

相较传统地面调查的灾情获取手段，遥感技术更为便捷、宏观，利用全极化SAR可以全天时、全天候工作，还可以记录4种不同极化组合方式下的后向散射信息，包含更丰富的对地观测内容，有助于人们更好地理解地物特征，快速实现建筑物震害灾情的获取，在地震应急中发挥着越来越重要的作用(金鼎坚等，2012；崔丽萍等，2016a)。

目前，各国学者大多利用震前、震后多幅全极化SAR影像变化检测来提取震害信息。Yamaguchi等(2010)对多幅震前、震后的全极化SAR影像做极化分解，对影像散射信息的变化进行统计分析，从而有效地提取了震害信息；Esra等(2011)基于特征量KL距离的变化，实现极化SAR影像的变化检测，得到了较好的检测效果；张红等(2009)提出了针对特定目标区域的一种基于极化SAR影像的变化检测方法；吴樊等(2009)提出极化似然比检验算法，通过Wishart分布的地物目标进行分类后，进一步实现变化检测。

由于震前全极化SAR数据缺乏，不能满足地震应急的实时性需求，因此对基于震后单时相全极化SAR影像数据的震害评估方法进行研究更具有必要性。

基于震后单时相全极化SAR影像数据的震害提取方法主要分为基于统计特征提取法和基于散射机理提取法两类。基于统计特征的震害提取方法是基于统计模型进行影像的极化散射矢量和极化协方差矩阵的分布计算，实现影像中地物目标的提取(龚丽霞等，2013)。Lee等(1994)基于Wishart分布定义各个类别之间的距离，利用最大似然法对目标数据进行非监督提取。基于散射机理的提取方法是根据散射机理与目标地物震害的对应特征，对目标地物的散射过程进行分解，实现震害的提取。Huynen(1978)使用Muller矩阵对不同的地物目标散射进行分解，提出奇次散射、偶次散射以及体散射3种散射机理，其中建筑物的散射特性归结为偶次散射；Cloude等(1997)对极化特征值进行分解，将目标特征的散射过程转换为相干矩阵，并对3种相干矩阵加权求和，从而表征不同的目标；Freeman等(1998)基于协方差矩阵将相干目标分解为表面散射、二次散射和体散射3种散射成分；Yamaguchi等(2010)针对建筑物的散射特性，在Freeman等(1998)分解的基础上，建立了适用于建筑物极化特征提取的模型，将螺旋散射成分加入其中，形成四成分散射模型。

目前，研究极化SAR提取建筑物震害的定量化指标较少。本文实验数据为Radarsat-2全极化影像，利用Pauli-Wishart监督分类方法对极化分解提取的建筑物震害信息进行分类识别，提取出未倒塌建筑物和倒塌建筑，依据提取分类结果，以街区为单元进行评估，并确定评估的精度。

1 极化特征提取与分类方法 1.1 Pauli分解

Pauli分解的核心是针对不同地物具有不同极化这一基本特征，通过相干矩阵来表征目标的极化特征，从而实现对目标地物的提取。地物目标的散射特性可以由HH、VV、HV、VH四种极化散射之和表示，在数学上可由4个2×2的基本散射矩阵(即Pauli基)表示，即

$ {\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{a}}} = \frac{1}{{\sqrt 2 }}\left[ {\begin{array}{*{20}{l}} 1&0\\ 0&1 \end{array}} \right] $

(1)

$ {\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{b}}} = \frac{1}{{\sqrt 2 }}\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&{ - 1} \end{array}} \right] $

(2)

$ {\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{c}}} = \frac{1}{{\sqrt 2 }}\left[ {\begin{array}{*{20}{l}} 0&1\\ 1&0 \end{array}} \right] $

(3)

$ {\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{d}}} = \frac{1}{{\sqrt 2 }}\left[ {\begin{array}{*{20}{c}} 0&{ - i}\\ i&0 \end{array}} \right] $

(4)

其中，S_a表示单次或奇次散射，代表球体、平坦表面和三面角反射器等的散射机制，在本文中，主要对应于损坏建筑物所产生的散射现象；S_b和S_c表示二次或偶次散射，分别代表 0°和45°方向角，主要对应未倒塌建筑物；S_d表示的散射机制是现实中不存在的不对称分量，通常不予考虑(薛腾飞，2017)。

极化散射效应则可以采用如下复矩阵表示

$ \begin{array}{l} \mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} {{S_{{\rm{HH}}}}}&{{S_{{\rm{HV}}}}}\\ {{S_{{\rm{VH}}}}}&{{S_{{\rm{VV}}}}} \end{array}} \right] = a{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{a}}} + b{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{b}}} + c{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{c}}} + d{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{d}}} \end{array} $

(5)

式中，a、b、c、d分别代表各分量的复系数，元素S_ij表示发射为i极化、接收为j极化(i，j=H或V)的回波信号分量。

根据互易性定理，即S_HV=S_VH，极化散射矩阵S矢量化可得到散射矢量D，即

$ D=[a b c]^{T}=\frac{1}{\sqrt{2}}\left[S_{\mathrm{HH}}+S_{\mathrm{NV}} S_{\mathrm{HH}}-S_{\mathrm{VV}} 2 S_{\mathrm{HV}}\right]^{T} $

(6)

根据$\mathit{\boldsymbol{T}} = \left\langle {D\; {D^H}} \right\rangle $，由极化散射矢量D转换得到极化相干矩阵T，其中，T₁₁、T₂₂、T₃₃分别对应极化相干矩阵T主对角线上的值，可表示为

$ T_{11}=\frac{1}{2}\left|S_{\mathrm{HH}}+S_{\mathrm{VV}}\right|^{2} $

(7)

$ T_{22}=\frac{1}{2}\left|S_{\mathrm{HH}}-S_{\mathrm{VV}}\right|^{2} $

(8)

$ T_{33}=2\left|S_{\mathrm{HV}}\right|^{2} $

(9)

根据极化相干矩阵，可得到极化Pauli分解总功率，即

$ \operatorname{SPAN}=\left|S_{\mathrm{HH}}\right|^{2}+2\left|S_{\mathrm{HV}}\right|^{2}+\left|S_{\mathrm{VV}}\right|^{2} $

(10)

1.2 Wishart监督分类

在极化SAR图像监督分类中，可以假设样本数据和待分类数据服从同一分布，因此，可以采用Wishart分类算法进行分类(翟玮，2016)。对N个视数的独立极化SAR数据，构造矩阵M (李广，2018)，即

$ \boldsymbol{M}=\sum\limits_{i=1}^{N} \delta(i) \delta(i)^{H} $

(11)

其中，δ表示复散射矢量，H表示共轭转置。

矩阵M服从复Wishart分布，对应的概率密度函数为

$ P(\boldsymbol{M})=\frac{|\boldsymbol{M}|^{N-Q} \exp \left(-N {Tr}\left(\left(E^{-1}(M) M\right)\right)\right)}{K(N, Q) E^{N}(M)} $

(12)

其中，Tr为矩阵M主对角矩阵元素的和，Q为极化系统的参数，本文Radarsat-2为单站互易极化，Q取值为3；K(N，Q)则按下式取值

$ K(N, Q)=\pi^{\frac{1}{2}\left(Q^{2}-Q\right)} \Gamma(N) \cdots \Gamma(N-Q+1) $

(13)

其中，Γ为Gamma函数。

采用已知倒塌和未倒塌建筑物样本，通过上述方法确定M的概率密度函数后，再通过最大似然法计算M到各类别中心的距离，通过类别间距离的收敛判断，迭代得到倒塌和未倒塌建筑物的分类结果。

1.3 遥感震害指数

在倒塌和未倒塌建筑物分类基础上，为进一步评估震害程度，采用下式计算每个街区的平均遥感震害指数

$ \overline{D_{i}^{R S}}=\frac{\sum\limits_{j} d_{i j}^{R S} s_{i j}}{\sum\limits_{j} s_{i j}} $

(14)

其中，d_ij^RS表示第i类遥感解译的建筑物破坏等级为j对应的震害指数，对于遥感震害评估，建筑物震害等级一般划分为2类(倒塌取值为1，未倒塌取值为0)；s_ij表示这一空间单元内第i类建筑物遥感破坏等级为j的建筑物区域的面积或像素数。

2 分解提取建筑物震害流程

按照上述极化SAR分解、震害分类和震害程度评估方法，确定建筑物震害提取流程，如图 1所示。

(1) 预处理：首先将原始Radarsat-2全极化数据转换成极化相干矩阵和极化协方差矩阵，然后对影像进行多视、滤波等处理，减少斑点噪声的影响，最后对影像进行地理编码，并与光学影像做配准。

(2) Pauli分解：对SAR图像进行极化分解，生成表面散射(Pauli_Odd)、0°角偶次散射(Pauli_Vol)和45°角偶次散射(Pauli_Dbl)3种散射机理的影像，将分解得到的极化特征影像波段叠加，融合生成假彩色影像图。

(4) Wishart分类：基于波段叠加影像图选取样本，计算每个类别的平均协方差矩阵、概率密度分布函数、类间距，使用Wishart分类器距离公式进行最大似然监督分类。

(5) 震害程度评估与分类精度验证：根据提取的建筑物震害级别，以街区为统计单元，按不同类别的像素数及其震害级别，计算遥感震害指数；以地震现场震害实地调查或遥感目视判读结果评估的震害程度为参照(真实分类)，验证本方法震害程度评估结果的精度。

3 实验及结果分析 3.1 实验数据

选择2010年4月14日青海玉树7.1级地震发生后一天(4月15日)获取的Radarsat-2全极化影像作为实验数据，如图 2所示。该影像距离向分辨率为13.295m，经过多视处理后，方位向分辨率与距离向分辨率保持一致。实验区域位于玉树县城区，该城区建筑物大约10000栋，以单层土木、砖木结构平顶为主(占比75%)，其次为多层砖混、框架结构房屋(占比24%)，以及少量厂房(郭华东等，2010；王晓青等，2013)，典型结构类型现场拍摄照片见图 3。

3.2 实验结果与分析

结古镇主城区Radarsat-2影像的Pauli分解三通道分量及假彩色合成图，如图 4所示。根据合成结果，选取多数倒塌建筑物区域与未倒塌建筑物区域样本，如图 5所示。与光学影像对比可以发现，倒塌严重的街区对应着分解合成图中发生奇次散射的区域，视觉展示上为蓝、绿色调为主合成的颜色(偏青色)，后续选取样本时，该区域标记为多数倒塌建筑物区域；光学影像中建筑物完好的区域，对应分解合成图中偶次散射的区域，视觉展示上为红、绿色调为主合成的颜色，后续选取样本时，该区域标记为多数未倒塌建筑物区域。

选取完样本后，使用Wishart分类器进行迭代分类，得到分类结果，并按照街区尺度进行分割，如图 6所示，分类结果图中红色区域为多数倒塌建筑物区域，绿色区域为多数未倒塌建筑物区域。

3.3 震害评估

以玉树主城区227个街区为统计单元，根据Wishart分类得到的倒塌和未倒塌建筑物结果，按照式(13)计算每个街区的平均遥感震害指数。

需要说明的是，本文分解基于散射和二面角机理，对于建筑物结构类型的抗震能力无法计算，故无法得到各结构类型遥感震害指数之间的准确关系。因此，本文未对建筑物进行结构类型分类，将综合遥感震害指数简化为平均遥感震害指数，其中，2个震害等级的建筑物区域平均震害指数分别取为0(倒塌)和1(未倒塌)。

参考崔丽萍等(2016b)的遥感震害等级划分方法，根据计算确定的各街区遥感震害指数，以0.55为阈值，将街区震害程度分成2类，多数未倒塌建筑物街区(类别0)取值为0~0.55，多数倒塌建筑物街区(类别1)取值为0.55~1，得到基于SAR影像的建筑物震害分类图，见图 7(a)。各街区不同震害程度的面积及遥感震害指数统计结果，见表 1。

在玉树地震震后应急期间，中国地震局地震预测研究所组织开展了基于高分光学遥感影像的建筑物震害(倒塌、局部倒塌、未倒塌3类)目视解译，共计9889栋建筑物(王晓青等，2013)。采用式(14)计算街区震害遥感指数，并以0.45为阈值分成多数倒塌和多数未倒塌2类(图 7(b))，该结果作为本文基于全极化SAR影像分解分类提取的震害分类评估结果的比对参考。

结果表明，基于SAR影像的震害分类结果与光学影像目视解译图像结果基本吻合，其中破坏较严重的区域位于研究区西部及南部。以光学影像目视解译的建筑物震害分类结果作为真实结果，统计不同建筑物震害程度的街区数，对SAR影像震害分类结果精度进行检验，其混淆矩阵如表 2所示。

精度分析结果表明，本文分类方法的准确率达到81%，一致性系数kappa值为61%。从精度验证结果可以看出，本文结果与光学影像目视解译结果具有较高的一致性，提取与整体分类效果较好。

4 结论

本文根据倒塌建筑物和未倒塌建筑物的全极化SAR散射机理，建立基于建筑物震害特征Pauli分解方法和基于Wishart监督分类的建筑物倒塌程度分类方法，并采用遥感震害指数评估街区单元的震害程度。选择玉树县城区2010年4月14日玉树7.1级地震震后极化SAR影像，进行了SAR图像去噪和几何校正、通过Pauli分解得到散射特征图像，利用极化分解监督分类方法，提取了玉树县城区的建筑物震害信息，并进行了以街区为单元的遥感震害指数评估以及建筑物多数倒塌和多数未倒塌街区分类。将本文分类结果与基于震后高分光学影像建筑物震害提取的建筑物多数倒塌和多数未倒塌街区分类结果进行对比验证，结果表明，本文所建立的极化SAR建筑物震害提取结果总体精度可达81%，Kappa值达到61%，精度评价指标可以看出本文的方法在一定程度上实现了建筑物震害信息的提取。

本文的方法实现半监督全自动化提取，相对于监督分类的耗时耗力和人工解译中易造成的错选和漏选等问题，效率更高，对于地震应急具有重要意义；此外，通过低分辨率单幅的全极化SAR影像便可实现较好震害提取效果。综上所述，本方法有较好的实际使用价值。今后，随着更多较高分辨率全极化SAR数据的方便获取，可以针对单体建筑物进行更为精确的提取。