基于高分辨率全极化SAR图像的取向角校正方法

孙 翔 宋红军 王 宇② 李 宁

①(中国科学院电子学研究所 北京 100190)

②(中国科学院大学 北京 100049)

③(河南大学 开封 475004)

1 引言

极化合成孔径雷达(Polarimetric Synthetic Aperture Radar, PolSAR)是一种获取丰富地物散射信息的手段，在军事和民用领域有着广泛的应用和研究价值[1]。极化目标分解是极化SAR图像解译的一个重要分支，是目标特征参数反演、目标识别和图像分类的理论基础[2]。极化散射矩阵将目标散射的能量特性、相位特性以及极化特性统一起来，相对完整地描述了雷达目标的电磁特性[3]。目标分解理论最早由Huynen提出，它利用极化散射矩阵揭示散射体的物理机理，促进极化信息的充分利用[4]。而后Krogager[5]、vanZly[6]、Cloude[7]和Pottier[8]等人做了大量的研究，为极化分解理论奠定了基础。

然而，地形目标(散射体)的散射依赖于散射取向、形状、介电特性、散射机制等。复杂地形表面的散射目标往往是随机取向的，会引起随机起伏的回波。随机取向和随机分布的散射目标难以分类。不同散射机制的不同取向的散射粒子可能会产生类似的散射；反之，相同散射机制的散射体则在随机的取向角下会造成不同的散射，导致混乱的分解结果与分类结果[9]。去取向概念的引入是为了减少随机波动取向的影响，将目标自身的物理特点突出地表现出来[10]。去取向的概念类似于斜坡补偿的概念，不同的是去取向适用于各种目标，包括目标和目标所在的背景。

2011年，为了将大片植被区域与城市建筑区域区分开来，Yamaguchi提出了使用根据取向角大小旋转相干矩阵后，再进行分解的方法[11]。该方法首次将去取向角引入了分解算法中。然而，城区取向角早在2000年就受到了科学家们的关注，那时就已经对方位向斜坡变化的补偿方法进行了研究[12]。他们提出了两种补偿方法，并利用反射对称的概念对估计算法进行了统一的分析[13]。尤其针对建筑物地区的取向角变化利用了后向散射模型进行了重点的研究[14]。偶次散射模型与奇次散射模型被推广使用于交叉极化项与非对角项中[15]。

为了研究不同特点的建筑物区域，2013年陈思伟等人提出了主导极化方向角(Dominant Polarization Orientation Angle,DPOA)的概念[16]。DPOA定义为在某一区域内所有像素点POA分布直方图的顶点所对应的POA值。在研究了多个不同DPOA值的建筑物区域去取向结果后，他们发现DPOA的大小与去取向角后建筑物地物分解的结果有较大的关系。经过去取向操作后基于模型的目标分解方法在|DPOA|＜22.5°时较为有效，体散射分量过估的问题得到有效解决，二面角散射分量得到了大幅额度提升，由11%提升到44%。然而，当|DPOA|＞22.5°时，分解结果中体散射分量只有小幅度的减少[16]。即使经过了常规去取向处理，在纯粹的建筑区域中体散射分量仍占主导地位[17]。

分解算法针对取向角问题改进至今，虽然误分解问题得到了一定的改善，但是由于基于模型的分解算法过于强调体散射分量，代表建筑物的偶次散射机制依然无法在城区占主导地位。为了解决这一问题，本文提出了基于高分辨率全极化SAR图像的取向角校正方法，使用了四川都江堰地区机载全极化图像进行算法验证。

2 极化取向角对极化分解结果的影响

极化方向角(Polarization Orientation Angle,POA)即极化取向角，定义为极化椭圆长轴与水平轴之间的夹角，如图1所示。

在经历极化基变换后，目标的散射相关信息不会发生变化。假设将目标散射矩阵变换极化基，变换后的极化基是原本的极化基沿视线旋转Ψ后得到的，那么新的极化基下的散射矩阵与旧极化基下的散射矩阵之间的关系为：

本文使用的数据是中科院电子所使用机载X波段全极化合成孔径雷达于2009年获得的四川省都江堰地区高分辨率数据。图像的斜距分辨率和方位分辨率均为0.5 m。入射角范围在20°到70°之间变化。图2(a)展示了该地区数据的Pauli分解结果。

图中用白色的边框标出了A, B, C 3个区域以便进行定量分析，其中A, C都是城市建筑物区域，B区域是平原低矮植被区域，存在少量独立的房屋建筑。从Pauli分解图可以看出，有大片的建筑物区域呈现绿色，反映的是体散射机制的性质，只有在C区域周边存在一部分较为整块的紫红色区域，反映了建筑物应显示的二面角散射特性。为了解决这一问题，常用的解决方式是去除使得T33最小的取向角后再重新分解，去除的取向角大小如式(7)所示。

由于计算取向角是使用了相干矩阵，所以去除取向角后进行分解时直接采用了基于相干矩阵的非相干目标分解方法，Yamaguchi分解方法。去取向后的Yamaguchi分解结果如图2(b)所示。

可以看出经过取向角旋转后的Yamaguchi分解图像依旧呈现大片绿色区域，即呈现了体散射特性，这与该区域的真实地貌不相符。出现该问题的原因在上节已经介绍。在这片区域中，大多数的建筑物与飞行方向之间的夹角较大，导致了极化旋转角度较大。大部分建筑物区域|DPOA|＞22.5°，使得传统的去取向角方法不再适用。为了解决这一问题，必须利用图像的其他特性将建筑物的取向角进行校正。

3 高分辨率图像POA校正算法

本文针对传统的去取向角方法无法校正由于取向角导致的分解体散射过估问题进行了研究，提出了新的适用于高分辨率图像的POA校正算法。算法分为两个部分：第1个部分利用高分辨率图像在城区取向角跳变的现象对需要进行POA校正的区域进行大致提取；第2部分针对提取出的区域进行迭代逼近的方法找出新的使得T33最小的POA。完整的算法流程图如图3所示。

下面将对算法的两个部分进行详细地介绍：

第1部分：大致提取需要进行取向角重校正区域，将无法通过常用POA计算方法去除取向角影响的区域选取出来。这个部分利用了POA的随机性进行研究。POA随机性的差异是对城市和其他地区进行分类的适当指标。基于POA的跳变现象，提出一个区分城市地区和平原地区的算法。算法流程如下：

步骤1 将POA分为5个部分。分界线分别为24°, 15°, 3°, –3°, –15°和–24° (图4(a)所示)；

步骤2 按照传统方法计算每个像素点的POA大小。将每个像素点根据POA的大小参照5类进行编号；

步骤3 假设一个像素为参考像素(Reference pixel)A，将其上下左右4个像素的编号与A的编号进行比较。此时定义一个新的跳变参数(Outburst Parameter, OP)。如果所有4个相邻像素的编号都与A相邻或相同，则参考像素A的跳变参数记为0；否则，记为1。如此计算图像中每个像素点的跳变参数(图4(b)所示)；

步骤4 在图像中放入一个9×9的窗，将这个窗中跳变参数为1的像素的个数记为窗中心像素的异构参数(Heterogeneous Parameter, HP)；

步骤5 移动窗的位置，再次执行步骤4，直到得到所有像素点的异构参数为止；

步骤6 异构参数大于10的像素点被归入需要进行POA重新计算的区域。

以图4(b)中标记的A, B两个像素点作为示范，进一步解释算法流程。参考像素A呈蓝色，根据图4(a)被标注为4。参考像素A相邻的4个像素点分别呈橙、红、蓝、绿色，表示它们分别被标注为1, 5,4, 3。其中左、右两个像素点的编号与参考像素A的编号相邻，下面的像素点与参考点的编号相同，位于参考点A上方的像素的编号与A的编号既不相邻也不相同，所以像素点A相邻像素中存在不连续的点，故A的跳变参数为1。同理，参考点B也按照同样的步骤进行分析。参考像素B的编号为1，它相邻4个像素点的编号分别为2, 2, 1, 5。对比图4(a)可以看出参考像素B的编号与其4个相邻像素的编号相邻或者相同，所以像素点B的跳变参数为0。

跳变参数的数值象征了像素点是否为跳变点，异构参数的大小表征的是像素点的POA随机性，异构参数的范围为[0, 81]。理论上，异构参数在城市建筑物密集地区较大，在平坦的草原等地区较小。图5展示了都江堰地区数据的异构参数值。图中像素点越亮，异构参数越大。可以看出在城市区域像素点的异构参数较大，在平原地区，高亮度的像素点明显减少。这与都江堰地区光学图像(图6)展现的地物相符。

在算法的第1部分中存在两个关键点。首先，POA种类数量的确定十分关键，种类的多少，如何划分会影响到跳变点个数的多少与跳变的可信程度。如果划分的种类数量过少，城市地区就无法体现POA频繁变化的特征。如果划分种类的数量较多，平坦区域将显示出与城市区域相同的随机性。因此，组的数量和它们之间的边界是至关重要的。图2(a)中3个区域分别代表不同的地物，研究这3个区域的POA分布对确定POA种类与分类原则十分重要。

图7展示了A, B, C 3个区域的POA分布情况。其中图7(b)代表了平原区域的取向角分布。可以看出，在平原区域，大多数像素点的POA都处于–15°到15°之间。因此，POA不属于该范围内的像素点基本属于城市建筑区。代表区域B的直方图的顶点位于0°，而在图7(a)与图7(c)中，直方图的顶点分别位于–3°和2°处。因此，POA处于–3°到3°的像素点，多属于B区域所代表的类型，而A区域所代表地物类型的像素点的POA大多处于–15°到–3°之间，C区域所代表地物类型的像素点的POA大多处于3°到15°之间。综合以上所述的情况，POA以24°, 15°, 3°, –3°, –15°和–24°为分界线分为5大类是合理的。

算法第1部分中另一个关键点是异构参数的阈值。根据上述步骤，在理想平原区域内，异构参数应为0；在另一极端情况下，异构参数应达到81的最大值。如果阈值设置得太低，会将平原区域一些较高大且密集的植被像素归类为城市像素；相反，若阈值设置太高，建筑边缘的像素将不被算入城市区域。所以确定阈值非常重要。因此需要对典型城区和平原区进行研究，图8显示了A, B和C区域的异构参数分布情况。

通过比较图8中的3个直方图可以看出，B区域像素的异构参数多处于0到10之间，而在图8(a)与图8(c)中异构参数相对平均地分布在0到30之间。为了确定阈值，在整个图像上尝试了从7到12的异构参数阈值，并将结果与光学照片(图6)进行了比较。图9中给出了7到12的异构参数阈值时，都江堰地区图像的状态。考虑到相对集中的大面积城市建筑区域以及平原地区的零散分布的建筑物，最终选择了10作为异构参数的阈值。

第2部分：在上文中提到当|DPOA|＞22.5°时，使用传统取向角计算方法获得的取向角无法使得交叉极化项T33达到最小值。

在使用传统去取向方法是，无论是相干分解还是非相干分解都无法在网格状的城市建筑物地区得到正确的分解结果。如图2所示，图像中的大部分建筑物都显示出典型体散射机制的特点。所以，为了找到一种使得交叉极化项最小化的方法，需要采用迭代逼近法来获得最适合的取向角。针对每个像素点的具体操作流程如下：

步骤1 在–24°到24°之间以1°的步长逐次对相干矩阵进行旋转，获得不同角度旋转下的相干矩阵；

步骤2 比较根据不同角度旋转后的所有相干矩阵的交叉极化T33项，记录下使得T33最小的两个角度值标记为α1,α2；

步骤3 将α1,α2之间的角度等分为3份，将两个等分点记做β1和β2；

步骤4 将初始相干矩阵根据α1,α2,β1和β2进行旋转，比较4个矩阵的T33值。将最小的两个T33值所对应的角度定义为新的α1,α2；

步骤5 返回步骤3进行计算，直到|α1-α2|＜0.1°跳出循环；

步骤6 定义POAnew。

4 实验结果分析

通过本文提出的取向角校正方法计算得到的取向角与传统方法得到的结果在城区的差别较大，如图10所示。

使用POAnew进行去取向操作后对都江堰地区全极化SAR图像进行分解，结果如图11所示。与使用传统算法计算POA后使用Pauli分解(图2(a))与Yamaguchi(图2(b))分解算法后得到的分解结果相比，新算法针对体散射过估的城市建筑物区域进行了POA校正，优化了分解结果，有效减轻了分解后城区显示体散射特性的现象。

本文提出的校正算法使得分解结果在保持平原、草坪以及植被地区的体散射机制处于主导地位的同时，将城市建筑物区域受到的POA影响降到最小，使得城区的二面角散射分量处于主导地位。图12展示了A, B, C 3个具有代表性的区域使用传统方法和本文方法计算取向角，旋转后的交叉极化T33项。

3个区域的T33平均值如表1所示。

表1 A, B, C区域中交叉极化T33项平均值Tab.1 AverageT33of districts A, B and C after rotation using POAs calculated by the traditional method and the new method

通过图12和表1对T33的分析可以看出，在经过本算法进行取向角校正之后，T33在建筑物区域(A,B区域)有明显的下降，在平原地区(C区域)基本保持不变。验证了本文提出的算法能够有针对性地对城市区域的POA进行校正，并获取最小交叉极化分量的结果。

使用新算法进行取向角校正后的分解结果与传统算法的分解结果对比如图13所示。

首先针对城市区域进行结果分析，区域A与区域C是城市区域。在使用Pauli分解算法时，区域A完全体现了体散射特性，呈绿色(图13(a))；在Pauli分解的区域C中，右半部分的的建筑物呈紫红色，体现了二面角散射机制的散射特性，而左下角的建筑物呈现绿色(图13(c))。通过本文提出的POA校正算法后进行分解，区域A与区域C都呈现了红色，即建筑物区域都正确地显示出了二面角散射特性(图13(d)，图13(f))。

区域B是典型的平原地区，在区域B中存在少数独立的建筑，在图13(b)中可以看出，在Pauli分解结果中，区域B呈现了绿色；在使用本文提出的POA校正算法后进行分解，区域B的大部分像素点仍然呈绿色，体现了体散射特性，但是少量的建筑物呈现红色，准确地表现了建筑物的二面角散射特性(图13(e))。

图11(b)展示了使用本文算法进行取向角校正后的Yamaguchi分解结果，通过与图2(b)的对比可以看出，本文方法使得城市区域显示出了较明显的二面角散射特性，呈红色。下面对3个特定区域进行定量分析。表2给出了A, B, C 3个区域表面散射(Ps)，体散射(Pv)和二面角散射(Pd)占总功率的百分比。从表中可以看出，通过本文提出的POA校正算法后进行分解，3个区域的二面角散射功率占总功率百分比得到了提升，在建筑物较多的A, C两个区域提升幅度较大，变化尤为明显。

和传统的POA计算方法进行比较，对所获得的分解结果进行分析可以看出，本文提出的对高分辨率区域进行POA校正的算法，能够使得分解结果更准确。本算法可以有效缓解|DPOA|＞22.5°区域分解错误的问题。

表2 A, B, C区域中各散射分量占总功率的百分比Tab.2 Percentage of scattering powers with traditional method and the new method

5 结论

基于高分辨率图像取向角跳变的特性，本文提出了一种适用于城区体散射过估问题的取向角校正方法。该方法以取向角随机性为基础，考虑高分辨率图像地物分界明确，散射特性清晰对取向角的影响，分两步对城区取向角进行校正。本文使用X波段全极化机载SAR实验数据对算法进行了验证，获得了与真实地貌相符的分解结果，有效缓解了城区体散射分量占主导的问题，验证了算法的有效性。