基于混合矩的极化SAR图像K分布模型参数估计新方法

崔浩贵 刘 涛 蒋宇中 高 俊

1 引言

极化合成孔径雷达(Polarimetric Synthetic Aperture Radar, PolSAR)具有全天时全天候的观测能力和多维信息遥感能力，是重要的遥感信息源。杂波统计建模及其参数估计方法是PolSAR图像目标检测和分类识别等图像解译手段的基础课题[1,2]。在高分辨情形下，地面或者海面散射杂波一般服从某种非高斯分布。针对PolSAR图像杂波非高斯建模问题，目前研究最多的是多变量乘积模型。该模型将目标的散射表示为一个表征目标雷达散射截面的纹理分量和一个表征相干斑的复高斯矢量的乘积，并且这两个变量之间相互统计独立。在多变量模型中，Gamma分布、逆 Gamma分布、Fisher分布和α稳定分布等非高斯分布常常用来对纹理分量进行建模[37]-。

K分布模型是目前应用最广泛的PolSAR图像统计模型，它是在纹理分量为Gamma分布时根据多变量乘积模型得到的。寻找快速、准确的参数估计方法是K分布研究的核心问题。传统的PolSAR图像K分布参数估计方法是利用单极化SAR的方法分别估计出每个极化通道强度图像的参数，然后将各通道的参数求平均。常用的单通道K分布参数估计方法为矩估计法[8]和最大似然估计法[9]。文献[10]的研究表明在矩估计法中采用分数阶矩能取得更好的估计效果。文献[11]提出了基于 zrlg(z)的K分布参数估计方法，该方法在 z l g(z)情形下有解析的形式。文献[12]对文献[11]的方法在不同r值下的性能进行了仿真比较，指出0

上述基于单通道的估计方法实际上只用了PolSAR图像协方差矩阵的对角元素，没有用到各极化通道间的相关信息。在单通道方法的基础上，针对 K分布下 PolSAR图像协方差矩阵的分布特性，以矩阵为变量对参数估计方法进行了研究。文献[13]将文献[10]的方法扩展到多维情形，给出了基于协方差矩阵行列式值二阶矩的参数估计方法。文献[14,15]提出了基于协方差矩阵的对数累积量(Matrix Log-Cumulants, MLC)的参数估计方法，该方法本质上是协方差矩阵行列式值的Mellin变换。研究表明基于对数累积量的估计方法在已有方法中是估计精度最优的，并已获得广泛应用[13,16]。该方法的缺陷在于没有解析的表达式，虽然根据表达式的单调性能较快地找到方程的解，但是在处理海量的PolSAR图像数据时，其运算速度仍有一定的局限性。本文将文献[11]的估计方法扩展到全极化多视数据的情形，提出了基于|z |rlg|z|混合矩(ZrLZ)的参数估计方法，该方法在r=1/d时有解析表达式。并且对ZrLZ方法在r取值不同值时的估计性能进行了仿真分析。最后用仿真数据和实测数据对本文方法与已有方法进行了比较，结果表明 ZrLZ方法对不同α值区域的适应能力更好，并且在r=1/d时的计算速度明显优于对数累积量方法。

2 PolSAR图像K分布统计模型

PolSAR图像中每个像素点可用目标矢量表示，在互易条件下可写为

多变量乘积模型是由一个表征目标雷达散射截面积的纹理分量 t和一个表征相干斑分量的复高斯矢量y组成，它们是相互独立的，即目标矢量为

其中复高斯矢量y的概率密度函数为

这里d是复高斯矢量的维数， Γ = E (y yH)，E(⋅)表示随机变量的数学期望。

实际中往往通过多视处理来抑制相干斑，该过程可由式(4)表示。

这里L为视图数，即进行平均的像素点的个数，k代表第k个像素点。如果视图数不是很多，往往假设多视处理的像素点具有相同的纹理分量，即 ()X k独立于k，此时式(4)可简化为

其中

式中Y为复高斯矢量y对应的协方差矩阵，Y服从Wishart分布，即 Y ～ Wd(L,Γ)，其概率密度函数为

其中d是Y的维数， Γ = E (y yH)，函数 Γd(L)为复数形式的多变量Gamma函数

当式(5)中纹理分量服从均值为1的Gamma分布时，其乘积模型服从K分布，即 Z ～ K (L, Σ, α )，其概率密度函数为

其中L为视图数， Σ = E (Z), α为Gamma分布的形状参数，Kα-Ld(.)为第α-Ld阶的第2类修正Bessel函数。

3 基于 |z |r lg|z |混合矩的 K 分布参数估计方法

假定高斯假设下协方差矩阵测量值Z和期望值Σ的行列式值分别为γ=Z,Υ=Σ , d是协方差矩阵的维数，L是视图数。已知随机变量η= ( 2L)dγ/Υ 为d个相互独立的 χ2分布变量的乘积构成，其自由度分别为2L, 2(L - 1),…,2(L-d+ 1 )[14,15]，考虑如式(5)所示的乘积模型，第 2类(Mellin变换)特征函数满足：

根据 χ2分布变量的Mellin变换[14]，可得行列式值随机变量 η = ( 2L)dγ/Υ 的Mellin变换为

将式(11)进行Mellin逆变换，可知高斯假设下协方差矩阵的行列式值的高阶矩为

Gamma分布 t ～ γ(1,α)的矩特征为

对于式(5)所示的乘积模型，由于纹理分量与相干斑分量是相互独立的，因此其矩特征存在如下关系：

因此可知K分布下协方差矩阵行列式值的矩特征为E(| |)r

Z 关于r求导，可得

将式(15)代入式(16)，可得

其中ψ(0)(x)为Digamma函数，并且ψ(0)(x)=Γ'( x ) /Γ( x )。对式(17)取不同的r值，并将两个式子相减，可以得到α的参数估计方法如式(18)所示

其中r1≠r2。当r=0时，式(17)实际上就是一阶对数累积量的表达式：

已有的研究表明将式(19)用于参数估计，其估计性能较好[7]，同时为了便于对估计方法进行分析，这里将 r2的值取0，此时式(18)可写为

我们将式(20)的估计方法称为基于|z |rlg|z|混合矩的参数估计法(ZrLZ)。对式(20)求导可得求解方程具有单调递减性，可以利用基于中值的搜索方法进行快速求解。另外我们发现当r=1/d时，由式(20)可得到ZrLZ(r=1/d)估计α的解析表达式

式(21)的推导应用了等式

4 实验结果与分析

这里对以下几种常用的K分布模型的参数估计方法性能进行比较：

(1)文献[15]提出的基于协方差矩阵行列式值二阶对数累积量的参数估计方法(Second-order MLC,SMLC)；

(2)文献[13]提出的基于二阶矩特征(Second MOMent, SMOM)的参数估计方法；

(3)文献[10]提出的基于单极化通道强度数据的分数阶矩的参数估计方法(Fractional MOMent,FMOM)；

(4)本文给出的ZrLZ方法。文献[12]指出在单极化通道情形下 0

上述几种方法中，FMOM为基于单极化通道的方法，其它的是基于协方差矩阵的估计方法。ZrLZ(r=1/d)和 SMOM 方法具有解析的表达式。SMLC和ZrLZ方法(r≠1/d)可以用基于中值的搜索方法进行求解，这里搜索精度设为 0.001。FMOM方法的求解方程不具备单调性，因此用遍历搜索的方法求解，选取的遍历范围为 α ∈ ( 0,100)，搜索精度为0.001。

4.1 仿真数据分析

K 分布仿真数据由如式(5)所示的乘积模型得到，视图数取 L = 1 0。形状参数α的取值范围一般在0.1(极不均匀区域)到10(海洋、平原等均匀区域)之间。图1给出了在不同的α取值下( α ∈ [ 0.1,10])，不同参数估计方法在样本数为512时的相对估计偏差与相对估计方差，仿真次数为 N = 1 0000次。其中相对估计偏差定义为对 | α︿ -α|/α求平均，相对估计方差为 v ar(α︿ ) /α2。首先比较ZrLZ方法在不同参数r下的估计性能，可以看到r=0.60时的估计性能较差，而r=0.05, r=0.20以及r=1/d这几种方法估计性能非常接近。因此下面在与常用的几种K分布模型的参数估计方法相比较时 ZrLZ方法只讨论r=1/d的情形。

同时从图1中可以看出当α值较小时，SMOM方法相对估计偏差和方差都非常大，此时该方法几乎失效。FMOM方法在α很小时效果最好，但是在α较大时，估计性能最差。另外可以看出本文的ZrLZ(r=1/d)方法在α值较小时，估计性能优于SMLC方法。当α值增大时，ZrLZ方法与 SMLC方法性能相当。因此可以得出ZrLZ方法对于不同α值的估计性能要更稳健。

图2为不同估计方法运算时间的比较，这里数据源为α=6的K分布仿真数据。计算机CPU为Inter E5700，双核3 GHz，内存大小为2 G。从图2可以看出在样本数小于512时，ZrLZ(r=1/d)方法的运算速度要明显快于其它方法。实际中往往采用滑窗的方法来对实测图像进行估计，为了适应不同的地物类型，每个滑窗包含的像素点一般都不会很多，因此该速度优势具有实用价值。

4.2 实测数据分析

传统的PolSAR图像统计模型的假设检验方法是对各极化通道的幅度或强度数据分别进行假设检验。该方法的缺点在于只考虑了协方差矩阵的对角元素，没有用到协方差矩阵中非对角元素的信息。最近，文献[16]提出了基于矩阵对数累积量的假设检验方法，实验结果表明该假设检验方法有效并且直观。对于模型参数未知需要从样本值估计得到的复合假设检验，统计量可由式(23)得到

假设检验中概率值(probability value, p值)是指在由 H0所规定的总体中随机抽样，获得等于及大于现有统计量的概率。根据式(23)所示的统计量可以通过蒙特卡罗仿真的方法得到p值[16]。进行假设检验时，向量k的维数要求大于m。对于K分布，其阶数大于1的MLC与参数Σ无关。另外L采用名义视图数代替，此时只剩下一个未知参数α，式(23)中取k=[,]T。

这里取的实测数据为PolSAR Pro软件中常用的AIRSAR系统San Francisco地区的L波段PolSAR图像(10 m×10 m)和 ESAR 系统 Oberpfaffenhofen地区的L波段PolSAR图像(3 m×3 m)，其中多视视图数假定为 L = 4 (2×2)。在两幅图像中各选取了4个大小为20×20的区域进行参数估计与假设检验。

两幅图像中 4个区域的选择及其（k2,k3）散点图如图3和图4所示。该散点图是在400个样本中随机选取200个计算（k2,k3），并且重复该过程50次得到的。图中的两条线分别为K分布和G0分布的理论（k2, k3）曲线[16]。San Francisco图像中选取的4个区域分别为城区、海洋、植被A和植被B，由图3可知城区的MLC靠近G0分布，海洋的MLC靠近Wishart分布，而植被A和植被B靠近K分布。Oberpfaffenhofen图像中选取的4个区域分别为城区、森林、植被A和植被B，由图4可知城区MLC在K分布和G0分布之间，实际上这是Fisher分布的区域[16]，而其它的3个区域都靠近K分布。

图1 在不同的α取值下，不同参数估计方法在样本数为512时的相对估计偏差与相对方差

图2 本文方法与原有方法在不同样本数下的运算时间比较

图3 San Francisco区域选取及MLC散点图

图4 Oberpfaffenhofen区域选取及MLC散点图

由于 FMOM 估计性能较差，这里不做比较。由仿真比较可知 ZrLZ无解析解的估计方法中r=0.60时估计性能较差，r=0.05与r=0.20估计性能相差不大。因此这里参数估计方法选取ZrLZ(r=0.20和r=1/d), SMLC和SMOM方法。下面我们在显著性水平α下，对各个区域在不同的参数估计方法下检验假设 H0：区域内样本值服从K(L = 4 ,Σ, α︿ )。

表1给出了San Francisco 4个区域不同算法下的K分布模型的参数估计结果及假设检验p值，从表中可以看出，对于城区，所有估计方法得到的 p值都无法通过显著性水平α=0.05的假设检验。从图3中散点图中可看到这是因为该区域的样本MLC更靠近G0分布，离K分布较远。对于海洋区域，所有的估计方法都能通过α=0.05的假设检验，但是相比较来说SMLC方法的p值最大，SMOM方法p值最小，而其它两种方法p值相近。在植被A区域，4种方法的p值都较大。而在植被B区域，SMOM方法的p值特别小，无法通过α=0.05的假设检验，而其它3种方法的p值都很大。

表2给出了Oberpfaffenhofen 4个区域不同算法下的 K分布模型的参数估计结果及假设检验 p值，从表中可以看出，对于城区，同样地所有估计方法得到的p值都无法通过显著性水平α=0.05的假设检验。对于森林区域，所有估计方法的p值都较大，但是相比较来说SMOM方法最小，而其它3种方法结果相差不大。对于植被A区域，SMOM方法无法通过α=0.05的假设检验，而其它 3种方法的p值都很大。对于植被B区域，SMOM方法无法通过α=0.05的假设检验，而其它方法都能通过α=0.05的假设检验。

表1 San Francisco 4个区域不同算法下的K分布模型的参数估计结果及假设检验p值

总结来说对于样本分布本身就离K分布较远的区域，例如两幅图中的城区，在α=0.05下四种方法都拒绝假设0H。对于接近 K分布的实测数据区域，SMOM方法的稳健性较差，其它3种方法的都能通过α=0.05的假设检验，并且p值相差不大。但是另一方面从算法复杂度来说，ZrLZ(r=1/d)具有解析的表达式，其计算速度要明显优于其它方法(见图 2)。

表2 Oberpfaffenhofen 4个区域不同算法下的K分布模型的参数估计结果及假设检验p值

5 结论

本文以|z |rlg|z|混合矩为基础推导了K分布的参数估计方法，并且该方法在r=1/d时有解析表达式。研究了最优r值的选择问题，仿真结果表明在选定的几个r值中，r取0.20时最优。最后用仿真数据和实测数据对本文方法与原有方法进行了比较，实验结果表明本文方法对不同α值的适应性最好，并且ZrLZ方法在r=1/d时运算速度的优势明显。但是最优r值的推导问题还有待研究。另外可以将基于|z |rlg|z|混合矩的参数估计方法的推导应用于G0分布和U分布等PolSAR图像模型，其推导结果及估计效果也值得研究。

[1] Collins M J, Denbina M, and Atteia G. On the reconstruction of quad-pol SAR data from compact polarimetry data for ocean target detection[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 591-600.

[2] Akbari V, Doulgeris A P, Moser G, et al.. A textural–contextual model for unsupervised segmentation of multipolarization synthetic aperture radar images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(4): 2442-2453.

[3] Lee J S, Schuler D, Lang R, et al.. K-distribution for multi-look processed polarimetric SAR imagery[C].Preceedings of the IEEE Geoscience and Remote Sensing Symposium, Pasadena, 1994: 2179-2181.

[4] Khan S and Guida R. On fractional moments of multilook polarimetric whitening filter for polarimetric SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing, 2014,52(6): 3502-3512.

[5] Cheng J, Gao G, Ding W, et al.. An improved scheme for parameter estimation of G0 distribution model in highresolution SAR images[J]. Progress in Electromagnetics Research, 2013, 134: 23-46.

[6] Sei T, Shibata H, Takemura A, et al.. Properties and applications of Fisher distribution on the rotation group[J].Journal of Multivariate Analysis, 2013, 116: 440-445.

[7] Liu T, Cui H G, Mao T, et al.. Modeling multilook polarimetric SAR images with heavy-tailed rayleigh distribution and novel estimation based on matrix log-cumulants[J]. Science China Information Sciences, 2013, 56(6): 1-14.

[8] Iskander D R and Zoubir A M. Estimating the parameters of K-distribution using higher-order and fractional moments[J].IEEE Transactions on Aerospace and Electronic Systems,1999, 35(4): 1453-1457.

[9] Blacknell D. Comparison of parameter estimators for K-distribution[J]. IEE Proceedings-Radar, Sonar and Navigation, 1994, 141(1): 45-52.

[10] Frery A C, Correia A H, and Da Freitas C. Classifying multifrequency fully polarimetric imagery with multiple sources of statistical evidence and contextual information[J].IEEE Transactions on Geoscience and Remote Sensing, 2007,45(10): 3098-3109.

[11] Blacknell D and Tough R. Parameter estimation for the K-distribution based on [zlog(z)][J]. IEE Proceedings-Radar,Sonar and Navigation, 2001, 148(6): 309-312.

[12] 胡文琳, 王永良, 王首勇. 基于 zrlog(z)期望的K分布参数估计[J]. 电子与信息学报, 2008, 30(1): 203-205.Hu Wen-lin, Wang Yong-liang, and Wang Shou-yong.Estimation of the parameters of K-distribution based on zrlog(z)[J]. Jounal of Electronics & Information Technology,2008, 30(1): 203-205.

[13] Doulgeris A P, Anfinsen S N, and Eltoft T. Classification with a non-Gaussian model for PolSAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008,46(10): 2999-3009.

[14] Anfinsen S N, Doulgeris A P, and Eltoft T. Estimation of the equivalent number of looks in polarimetric synthetic aperture radar imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(11): 3795-3809.

[15] Anfinsen S N and Eltoft T. Application of the matrix-variate Mellin transform to analysis of polarimetric radar images[J].IEEE Transactions on Geoscience and Remote Sensing, 2011,49(6): 2281-2295.

[16] Anfinsen S N, Doulgeris A P, and Eltoft T. Goodness-of-fit tests for multilook polarimetric radar data based on the Mellin transform[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(7): 2764-2781