高分辨率极化SAR影像KummerU分布非监督分类方法

朱 腾， 何汉武， 黄铁兰*， 张 坡

(1. 广东工贸职业技术学院测绘遥感信息学院， 广州 510510； 2. 广东工业大学机电工程学院， 广州 510039；3. 广东工贸职业技术学院计算机与信息工程学院， 广州 510510)

影像分类是遥感数据应用的重要研究内容之一[1]. 由于全极化合成孔径雷达(Fully Polarimetric Synthetic Aperture Radar)能够测定地物目标对组合极化波的相干成像，从而获得更多的地物散射信息，因此，无论在民用或军用领域，极化SAR影像分类解译技术均具有重要意义[2].

在早期的全极化SAR数据分布拟合研究中，LEE等[3]假设一个完全发育的相干斑是由大量服从高斯分布的散射体组合而成，进而根据中心极限定律推导出极化相干矩阵服从Wishart分布，并结合最大似然准则提出了经典的H/α-Wishart非监督分类算法. 然而，随着极化SAR系统成像分辨率的不断提高，影像分辨单元内的散射体数目急剧减小，基于中心极限定理的散射模型难以满足SAR影像的精确统计建模需求，造成分割、分类、目标识别与跟踪等后续应用的精度不断下降[4]. 针对这一问题，研究人员进行了大量的研究，如：使用Fisher分布对城市区域的高分辨率极化SAR数据进行了分布拟合，混合了Gamma分布与逆Gamma分布的复合函数对都市区域的回波拟合表现出了较好的适应性[5]；在Fisher分布的基础上，层次分割模型对高分辨率极化SAR数据的异质区域分割问题取得了不错的结果[6]；进一步将Fisher分布推广到协方差矩阵C的KummerU分布，并使用KummerU分布对模拟数据和L波段真实极化SAR数据进行了分割分类实验，实验结果表明使用KummerU分布可以对重头或长拖尾的分布进行建模，得到精度更高的分割结果[7]；利用KummerU分布与水平集方法实现了高分辨率极化SAR影像分割方法，在各向同性与各向异性区域均取得良好划分结果[8]；结合KummerU分布与马尔科夫随机场方法对CONVAIR数据进行了小范围分类，得到了区域一致性更好的结果[9]；使用KummerU分布对X波段TerraSAR和C波段AIRSAR数据进行了拟合实验，实验结果也表明KummerU分布模型的拟合精度高于Gamma分布、K分布、G0分布等传统拟合方法[10]. 已有的研究主要从分布拟合的角度论证了KummerU分布的有效性，着重关注分布函数、直方图和据估计值等统计信息的相似性，对于涉及空间信息、参数优化和全局寻优效率等方面的分类应用考虑较少.

针对以上研究现状，本文结合最大似然准则与KummerU分布提出适用于高分辨率极化SAR影像的KummerU距离非监督分类方法，同时在聚类寻优部分引入集群寻优算法——PSO算法[11]，以提高高分辨率化SAR数据非监督分类过程的全局寻优可靠性；并选择我国自研的X波段机载极化SAR影像，进行H/α-Wishart、PSO-Wishart、PSO-KummerU方法的分类对比试验，以验证PSO-KummerU方法在高分辨率极化SAR非监督分类方面的有效性.

1 极化数据分布与分类距离度量准则

全极化SAR数据的原始格式为Sinclair散射矩阵，根据发射与接收的极化波方向的4种组合[12]，其矩阵形式为：

(1)

其中，V表示垂直极化，H表示水平极化，如SVH表示以水平极化方式发射、垂直极化方式接收的回波强度. 为方便分析，一般将Sinclair矩阵以Pauli基表示为向量形式，然后转换为协方差矩阵C，对于多视极化SAR数据[13]，有：

(2)

其中，Sp表示Pauli基下的Sinclair向量，*T表示共轭转置，n表示多视处理的视数.

根据LEE等[3]的研究，当影像分辨率不高时，对于同质区域可以假设Sinclair向量由大量近似复向量叠加组成，根据大数定律有协方差矩阵C服从Wishart分布：

(3)

其中，Σ表示样本协方差矩阵的期望，q表示通道数(单基互易系统下q=3，双基系统时q=4)，n表示视数，K(n,q)为Gamma函数的连乘积[14]. 根据式(3)，利用最大似然估计即可得到协方差矩阵分类时的Wishart距离度量：

(4)

其中，Cω与C0分别表示聚类时第ω类的聚类中心、某一像素对应的协方差矩阵.

然而，随着极化SAR影像的分辨率不断提高，单像元对应的散射矢量相对减少，在诸如城区、茂密植被和森林等复杂区域，相干斑一致性假设已不再适用. 针对这一问题，HARANT等[7]使用非高斯的积模型进行拟合：

C=μΣ，

(5)

其中，μ表示纹理变量，Σ表示同质区域的协方差矩阵期望. 在积模型假设下，当μ取常数时，协方差矩阵C分布退化为Wishart分布；在最新的高分辨率极化SAR数据分布研究中，μ被表示为Gamma分布与逆Gamma分布的梅林卷积[15]：

(6)

其中，n和q分别表示视数和通道数；Γ(·)为Gamma函数；U(a,b,z)为KummerU函数积分项；a=np+M，b=1+np-n，z=nTr(Σ-1C)/Mm. 在n、M、m3个参数的调节下，KummerU分布可以更全面地描述纹理与同质区域. 如当M→∞时，Fisher分布退化为Gamma分布，KummerU分布退化为K分布；当L→∞且M→∞时，大数定律重新适用，KummerU分布退化为Wishart分布[7]. 依据式(6)与最大似然准则,可得到KummerU距离：

(7)

2 极化SAR影像PSO-KummerU非监督分类算法

针对高分辨率极化SAR影像不再适合Wishart分布的问题，本文引入KummerU分布作为SAR影像聚类过程的距离度量准则. 同时，聚类寻优过程采用粒子群优化算法框架[16]，减少H/α-Wishart分类采用的K均值迭代模式带来的噪声影响大、初始值敏感和全局寻优能力差等问题.

粒子群算法的核心是使用D维向量xi作为粒子i的位置参数，用于记录寻优问题在D维特征空间中的一组可行解，同时利用D维向量vi记录粒子的速度参数、D维向量Pi表示粒子群记录的个体最优、Gi表示全局最优记录，然后反复迭代更新粒子群的速度与位置，最终找到寻优问题的全局最优解[17]. 粒子群的速度与位置更新公式如下：

vid=w·vid+c1·r1·(Pid-xid)+c2·r2·(Gd-xid)，

x(i+1)d=xid+vid，

(8)

其中，r1与r2为步长控制参数，服从(0,1)上的均匀分布；c1、c2分别为个体、全局学习因子，默认值为2；w表示惯性权重，用于控制粒子速度在下次迭代中的比例[18].

针对全极化SAR影像聚类问题，本文采用的PSO-KummerU聚类方法包括初始分类与PSO迭代2个部分(图1). 其中分类迭代步骤的适应度函数决定了各个粒子的优劣程度，进而影响整个粒子群的寻优方向与效率，本文结合KummerU分布设计的适应度函数为：

(9)

其中，Ci为像素i的协方差矩阵，Σj表示像素i所属第j类的聚类中心对应的协方差矩阵，d(·)为距离度量准则.

图1 PSO-KummerU方法分类流程Figure 1 The workflow of the PSO-KummerU classification method

综上所述，本文采用PSO-KummerU非监督分类方法的具体流程如下：

(1)利用3×3Boxcar滤波对高分辨率SAR数据进行预处理，适当减小相干斑噪声影响；

(2)对散射矩阵进行Pauli分解，并计算得到协方差矩阵C与熵H、散射角α等极化特征参数；

(3)利用H/α分类对影像进行初始划分，同时初始化粒子群；

(4)对于每一个粒子，根据粒子的位置参数与KummerU距离对每个像元进行归类；

(5)计算每个粒子适应度，然后更新PSO的个体最优与全局最优记录；

(6)判断是否满足终止条件，是则终止运算并依据全局最优记录得到分类结果；否则更新粒子的位置与速度参数，迭代进入步骤(5).

3 实验结果及分析

为验证KummerU分布与PSO寻优算法对高分辨率极化SAR影像的有效性，本文采用了一组真实X波段高分辨率极化SAR影像[19]进行分类对比实验. 该数据采自中国海南陵水黎族自治县周边农田区域，由中国电子科技集团38研究所研制的X波段原型样机拍摄于2010年，影像分辨率可达0.5 m×0.5 m，实验采用的影像大小为804 px×573 px. 该区域主要包括植被、水塘、裸地和道路等地物，其中农作物主要包括条带状禾本植物(稻田耕地)、混杂状矮小植株(辣椒田)和木本经济作物(芒果林)，实验区的PauliRGB合成图与真实地物参照如图2.

图2 海南陵水实验区的PauliRGB合成图与真实地物参照Figure 2 The PauliRGB and real terrain reference of Hainan Lingshui experimental area

为验证粒子群智能算法与KummerU距离在高分辨率极化SAR数据聚类分析中的寻优能力，分别采用H/α-Wishart、PSO-Wishart、PSO-KummerU 3种非监督分类方法进行实验(图3). 为方便对比，3种方法在初分类步骤均采用H/α划分，将实验数据初始化为8个类别；粒子群参数也采用相同设置，粒子规模为50个，最大迭代次数为20次，学习因子c1=c2=2.0，惯性因子w从0.9逐渐变动到0.4.

图3 海南陵水实验区分类结果Figure 3 The classification results of Hainan Lingshui experimental area

由PauliRGB合成图(图2A)可知：(1) 相干成像的极化SAR数据呈现明显的斑点特征，而X波段高分辨率传感器极化SAR影像中的斑点则愈发细碎. (2)整体看来，田间的水塘因回波较弱而呈块状黑色；田间道路和裸土同样以低熵表面散射为主，颜色较暗；耕地中的裸土表面更为粗糙，因而交杂着更多的明暗斑点. (3)从植被形态来区分，作物区域类型主要表现为3种典型类别：植株较为矮小的辣椒田、呈条带状的耕地以及能够形成阴影区的芒果树林.

H/α-Wishart方法分类结果(图3A)与PSO-Wishart方法的分类结果(图3B)区别不大，整体来看：地物划分散乱，树林阴影与水体被合并为品红类，耕地的条带、农田区的裸土都未得到较好的划分，使用PSO算法进行解空间寻优也仅仅减少了树林阴影与水体的误分类;对于X波段高分辨率极化SAR影像，使用Wishart距离难以得到较好的聚类效果，引入PSO方法对聚类中心的优化效果也非常有限.

由PSO-KummerU方法的分类结果(图3C)可知：采用KummerU距离以后，虽然部分区域仍存在细小斑点，但从总体分类效果来看，PSO-KummerU方法较好地划分了树林阴影与水塘，阴影区仅有少量像点被划分为水体；作物按植株高度、密集程度可大致分为3类，裸土则分为了道路与耕地裸土2类. 存在少量错分、混分的细小斑点的原因为：一方面是由于非监督分类算法精度有限；另一方面是由于高分辨率数据对地面的采样能力更为精细，虽然部分区域调绘标注为同一类别，但该区域中的部分像素确实代表其他类型的地物，如耕地区域中实际包含植株阴影、强反射点和混杂体散射区等.

综上所述，采用KummerU分布进行高分辨率SAR影像非监督分类可以取得更好的分类效果.

对照图2B的真实地物参照，选取了5类典型地物进行分类精度分析，由对比结果(表1)可知：(1)使用传统的Wishart距离进行分类度量时，对于辣椒田、幼苗耕地这类复杂散射区域，Wishart距离难以进行有效的区分，分类精度为48%～62%，而用KummerU距离进行类间划分后，精度普遍提高到70%以上. (2)对于不同散射机制的大类区分，如划分粗糙表面散射的裸土、镜面反射的水体与复杂体散射的植被时，Wishart距离与KummerU距离都能做到基本的区分，如3种方法对于水体的区分精度都达到98.5%以上，茂密作物类型分类精度也能达到78%以上. 可见PSO-KummerU方法的总体精度提升主要由于KummerU分布能够更好地拟合与区分异质区域，从而能对作物进行精准分类.

表1 海南实验数据分类精度对比Table 1 The comparison of classification accuracy of Hainan data

4 小结

针对高分辨率数据不再满足传统的匀质区域假设的问题，本文通过采用异构杂波模型对数据分布进行拟合，利用KummerU距离对高分辨率极化SAR数据的异质区域进行准确分类；采用海南陵水农作物试验区国产X波段高分辨率极化SAR影像，对比了H/α-Wishart、PSO-Wishart与PSO-KummerU方法的分类结果. 实验结果表明：KummerU距离在高分辨率全极化SAR数据分类中的有效性明显强于Wishart距离，主要原因为KummerU分布能够更好地拟合和区分高分辨率极化SAR数据的异质区域.

目前，对于作物区域的分类主要参考其植株大小、茂密程度，后续研究将尝试根据作物的具体类别，分析不同种类作物、不同地物类型在KummerU分布下的参数估计，以监督分类的形式提高KummerU距离分类方法的分类精度.