CN106127236A - 基于狄利克雷mrf混合模型的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于狄利克雷MRF混合模型的极化SAR图像分类方法，主要解决现有技术无法自动地确定极化SAR图像中的类别数目的问题。其方案为：1.提取并归一化极化SAR图像的N个极化散射特征，建立归一化的极化散射特征空间；2.对归一化的极化散射特征空间逐点进行降噪，建立极化散射特征空间；3.利用极化散射特征空间初始化MRF模型；4.根据初始化的MRF模型估计极化SAR图像的先验参数和似然参数；5.根据估计参数估计新的标记场，直到达到最大迭代次数，并将新的标记场确定为极化SAR图像的分类结果。本发明提高了分类精度和同质区域分类的平滑度，更好地保持边缘信息，可用于极化SAR图像的目标检测与识别。

Description

基于狄利克雷MRF混合模型的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，特别涉及一种极化SAR图像的分类方法，可用于极化SAR图像的目标检测与识别。

背景技术

极化合成孔径雷达SAR是一种高分辨率成像雷达。其在民用和军事领域的广泛使用需要极化SAR图像解译技术作为支撑，而极化SAR图像分类是机器学习和数据挖掘领域的一项重要技术，也是图像解译的重要内容之一，它可以提供极化SAR图像的整体结构信息，揭示极化SAR图像的本质。近年来，极化SAR图像分类方法一直是该领域前沿研究的热点，其中随机场模型被认为是处理极化SAR图像分类问题的强有力工具，其优势在于能够在极化SAR图像分类过程中充分考虑各像素点间的空域相关性，从而提高算法的抗噪性，目前，基于随机场模型的极化SAR图像分类方法应用较多的是基于马尔科夫随机场MRF的极化SAR图像分类方法。

Wanying Song等人发表的论文“The WGΓDistribution for Multi-LookPolarimetric SAR Data and Its Application”(IEEE Geosci.Remote Sens.Lett.,2014,12(10):2056-2060)提出了将极化SAR数据的乘性模型引入MRF框架中实现极化SAR图像分类的方法。该论文利用极化SAR数据的乘性模型构建MRF框架中的似然概率，能够更精确地描述极化SAR图像的统计特征，从而克服了Wishart分布不具备描述极化SAR图像中复杂场景统计特征的能力。但是该方法的不足是，此类MRF模型需要训练数据学习模型参数，无法在缺乏极化SAR图像先验信息的前提下精确地确定极化SAR图像中的类别数目，只限于处理监督的极化SAR图像分类。

Peter Orbanz等人发表的论文“Nonparametric Bayesian image segmentation”(Int.J.Comput.Vision,2008,77:25-45)提出了利用狄利克雷MRF混合模型实现SAR图像的直方图聚类算法。该论文有效地结合了MRF模型和狄利克雷过程混合模型DPMM，并应用到SAR图像的直方图聚类算法中，视SAR图像中的类别数目为模型中的一个随机变量，根据模型和数据不断地对该随机变量进行更新，从而使得SAR图像的直方图聚类算法能够自动地确定SAR图像中的类别数目；且该论文中SAR图像空域结构信息的引入大大提高了模型的抗噪性能。但是该方法的不足是，此类模型是通过SAR图像特征向量的欧式距离来引入空域结构信息的，对极化SAR数据并不适用。

发明内容

本发明的目的在于克服上述已有问题的缺点，提出一种基于狄利克雷MRF混合模型的极化SAR图像分类方法，以在缺乏图像先验信息的前提下自动地确定极化SAR图像中的类别数目，并实现利用狄利克雷MRF混合模型对极化SAR图像的分类。

为达到上述目的，本发明的技术方案包括如下：

(1)输入极化SAR图像；

(2)从极化SAR图像中提取N个极化散射特征，并对其进行归一化，得到N个归一化的极化散射特征，建立归一化的极化散射特征空间F₁；

(3)以极化SAR图像中的每个像素点为中心像素点，选取大小为3×3的窗口，使用广义均值方法，对归一化的极化散射特征空间F₁逐点进行降噪处理，得到极化散射特征空间F₂；

(4)初始化模型参数

将极化SAR图像的每个像素点的类别标记的集合设定为原始的标记场X＝{X_i}_{i＝1,2,…,n}，其中X_i表示极化SAR图像中像素点i的类别标记，n表示极化SAR图像的像素点个数；初始化极化SAR图像的每个像素点的类别标记为1，即X_i＝1；

用极化散射特征空间F₂作为观测场Y＝{Y_i}_{i＝1,2,…,n}，其中Y_i表示极化SAR图像中像素点i的N维降噪的极化散射特征向量；

初始化形状参数a＝1，尺度参数b＝5，核宽参数σ＝1，平衡参数λ＝1，维度参数q＝3，第一正态参数μ₀初始化为N维零向量，第二正态参数Σ₀初始化为N维单位阵；

初始化中心参数α₀，即对Flevoland地区的极化SAR图像设定中心参数α₀＝2；对Oberpfaffenhofen地区的极化SAR图像设定中心参数α₀＝10；

将当前的迭代次数设为1，设定最大迭代次数为20；

(5)根据原始的标记场X、观测场Y和初始化的模型参数，估计极化SAR图像的先验参数s_ij和似然参数{μ_k,Σ₁}；

(6)估计极化SAR图像新的标记场

6a)根据先验参数s_ij、原始的标记场X和初始化的模型参数，得到类别标记的先验概率p₃：

$p_3=\left\{\begin{array}{cc}\frac{\underset{l=1,l\≠i}{\overset{n}{\Σ}}\mathrm{\δ}(X_i,X_l)}{{\mathrm{\α}}_0+n-1}\exp (-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\Σ}{s}_{ij}\left(1-2\mathrm{\δ}\left(X_i,X_j\right)\right)),& X_i\∈X_{-i}\\ \frac{{\mathrm{\α}}_0}{{\mathrm{\α}}_0+n-1},& X_i\∉X_{-i}\end{array}\right.,$

其中X_j表示极化SAR图像中像素点j的类别标记；X_l表示极化SAR图像中像素点l的类别标记；表示极化SAR图像中像素点i的邻域系统；X_-i表示极化SAR图像中除了像素点i的类别标记以外的其他像素点的类别标记的集合；δ(X_i,X_l)表示Delta函数，当X_i与X_l相等时，δ(X_i,X_l)的值为1，否则为0；δ(X_i,X_j)表示Delta函数，当X_i与X_j相等时，δ(X_i,X_j)的值为1，否则为0；s_ij表示极化SAR图像的先验参数；exp(·)表示求指数操作；

6b)根据似然参数{μ_k,Σ₁}、原始的标记场X和观测场Y，得到类别标记的似然概率p₄：

$p_4=\frac{1}{{\left(2\mathrm{\π}\right)}^{N/2}|{\mathrm{\Σ}}_1{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_k\right)}^T{\mathrm{\Σ}}_1^{-1}(Y_i-{\mathrm{\μ}}_k\left)\right),X_i=k,$

其中μ_k表示极化SAR图像中类别标记为k的像素点的均值向量参数；Σ₁表示极化SAR图像的协方差矩阵参数；(·)^-1表示求矩阵逆的操作；(·)^T表示求矩阵转置的操作；|·|表示求矩阵行列式的操作；X_i＝k表示极化SAR图像中像素点i的类别标记为k；

6c)根据类别标记的先验概率p₃和似然概率p₄，得到类别标记的后验概率p₅；

$p_5=\left\{\begin{array}{cc}{z}^{-1}{p}_4\·n_k^{-i}\·\exp (-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\mathrm{\Σ}}{s}_{ij}(1-2\mathrm{\δ}(X_i,X_j)\left)\right),& k=1,2,...,C\\ z^{-1}{\mathrm{\α}}_0\frac{1}{{\left(2\mathrm{\π}\right)}^N|{\mathrm{\Σ}}_1{\mathrm{\Σ}}_0{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_0\right)}^T{\left({\mathrm{\Σ}}_1{\mathrm{\Σ}}_0\right)}^{-1}(Y_i-{\mathrm{\μ}}_0\left)\right),& k=C+1\end{array}\right.,$

其中指的是极化SAR图像中除了像素点i以外的类别标记为k的像素点的个数；z为归一化参数；C表示原始的标记场X中的类别数目；

6d)通过极大化类别标记的后验概率p₅，估计得到极化SAR图像中像素点i的新的类别标记

${\stackrel{\‾}{X}}_i=\underset{k=1,2,...,C,C+1}{max}\left\{p_5\right\};$

6e)将新的类别标记的集合确定为极化SAR图像的新的标记场

(7)判断是否达到最大迭代次数，输出分类结果：

将当前的迭代次数加1作为新的迭代次数，判定新的迭代次数是否达到步骤(4)所设定的最大迭代次数20，如果没有达到，则设定新的标记场为原始的标记场X，返回步骤(5)；如果达到，则将新的标记场确定为极化SAR图像的最终的分类结果并输出。

本发明与现有技术相比，具有以下优点：

1)相比于基于MRF的极化SAR图像的监督分类方法，即Wanying Song等人提出的方法，本发明方法将狄利克雷过程混合模型DPMM引入到MRF框架中，可以克服在缺乏极化SAR图像先验信息的情况下难以确定极化SAR图像中的类别数目这一弊端，因此本发明更适合处理复杂场景的极化SAR图像的分类。

2)相比于Peter Orbanz等人提出的方法，本发明利用极化协方差矩阵之间的相似度估计极化SAR图像的先验参数，能增强模型与极化SAR图像的数据相关性，能使同质区域分类结果更平滑且能更好地保持边缘信息，进而获得更好的分类结果，因此本发明可以显著地提高极化SAR图像分类结果的精度。

附图说明

图1是本发明的实现流程图；

图2是用本发明和现有方法对Flevoland地区的极化SAR图像的分类对比图；

图3是用本发明和现有方法对Oberpfaffenhofen地区的极化SAR图像的分类对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1，输入极化SAR图像。

本发明选取如下两幅极化SAR图像：

Flevoland地区的极化SAR图像：图像大小为320×326像素；等效视数为4；分辨率为2.5652m×2m；雷达系统为AIRSAR；

Oberpfaffenhofen地区的极化SAR图像：图像大小为500×450像素；等效视数为2；分辨率为3m×0.89m；雷达系统为ESAR。

步骤2，提取并归一化极化散射特征，建立归一化的极化散射特征空间F₁。

2a)从极化SAR图像中提取N个极化散射特征F_r，r＝1,2,...,N，N＝17表示极化散射特征的个数，本实例所选取的N个极化散射特征如表1所示：

表1

通过极化SAR数据实验验证了层次后向策略应用于极化目标分解后，所提取的极化散射特征有助于提高分类精度，因此本发明提取用结合层次后向策略的Yamaguchi-Sato分解方法对极化SAR图像进行分解所得到的螺旋散射功率、体散射功率、二面角散射功率和表面散射功率作为极化散射特征F₁、F₂、F₃、F₄；

极化SAR图像的相干矩阵T有利于解译目标的散射机制，例如相干矩阵T的对角元素分别表示Pauli分解的系数，于是极化SAR图像的相干矩阵T的各个元素T₁₁、T₁₂、T₁₃、T₂₂、T₂₃、T₃₃有助于提高分类精度，因此本实例提取T₁₁、|T₁₂|、|T₁₃|、T₂₂、|T₂₃|、T₃₃作为极化散射特征F₅、F₆、F₇、F₈、F₉、F₁₀；

极化SAR图像的相干矩阵T的特征值与极化总功率同样是极化SAR图像分类的重要特征，因此需要提取极化SAR图像的相干矩阵T的3个特征值与极化总功率作为极化散射特征F₁₁、F₁₂、F₁₃、F₁₄；

用Cloude分解方法对极化SAR图像进行分解得到的3个极化散射特征：极化熵、平均散射角和反熵是最近十年内极化SAR图像分类中最常用的3个极化散射特征，其中极化熵对植被、树林、非常粗糙地面等目标是非常敏感的，平均散射角用于确定目标的主导散射，反熵用于区分植被与裸地等散射机制差异较小的类别，因此本发明提取极化熵、平均散射角和反熵作为极化散射特征F₁₅、F₁₆、F₁₇；

2b)通过下式分别对极化SAR图像的N个极化散射特征F_r，r＝1,2,...,N，进行归一化，得到N个归一化的极化散射特征

$\stackrel{\‾}{F_r}=\frac{F_r-min\left(F_r\right)}{\max \left(F_r\right)-min\left(F_r\right)};$

2c)将归一化的极化散射特征中像素点i的特征值的集合确定为极化SAR图像中像素点i的归一化的极化散射特征向量F_i，i＝1,2,...,n：

$F_i={\left\{\stackrel{\‾}{F_r}\left(i\right)\right\}}_{r=1,2,...,N};$

2d)将归一化的极化散射特征向量F_i的集合确定为归一化的极化散射特征空间F₁：

F₁＝{F_i}_{i＝1,2,...,n}。

步骤3，使用广义均值方法对归一化的极化散射特征空间F₁逐点进行降噪处理，得到极化散射特征空间F₂。

3a)以极化SAR图像中的像素点i为中心像素点，选取大小为3×3的窗口N_3×3(i)，使用广义均值方法，根据归一化的极化散射特征在3×3的窗口N_3×3(i)内提取极化SAR图像中像素点i的第r个降噪的极化散射特征G_r(i)：

$G_r\left(i\right)={\left(\underset{t\∈N_{3\×3}\left(i\right)}{\Π}\stackrel{\‾}{F_r}\right(t\left)\right)}^{\frac{1}{9}},$

其中表示第r个归一化的极化散射特征中像素点t的特征值；N_3×3(i)表示以像素点i为中心像素点的3×3的窗口；

3b)将N个降噪的极化散射特征G_r(i)的集合确定为极化SAR图像中像素点i的N维降噪的极化散射特征向量G_i：

G_i＝{G_r(i)}_{r＝1,2,...,N}；

3c)将N维降噪的极化散射特征向量G_i的集合确定为极化散射特征空间F₂：

F₂＝{G_i}_{i＝1,2,...,n}。

步骤4，初始化模型参数。

4a)初始化MRF模型(X,Y)：

将极化SAR图像的每个像素点的类别标记的集合设定为原始的标记场X＝{X_i}_{i＝1,2,…,n}，其中X_i表示极化SAR图像中像素点i的类别标记；初始化极化SAR图像的每个像素点的类别标记为1，即X_i＝1；

用极化散射特征空间F₂作为观测场Y＝{Y_i}_{i＝1,2,…,n}，其中Y_i表示极化SAR图像中像素点i的N维降噪的极化散射特征向量；

4b)初始化形状参数a＝1，尺度参数b＝5，核宽参数σ＝1，平衡参数λ＝1，维度参数q＝3，第一正态参数μ₀初始化为N维零向量，第二正态参数Σ₀初始化为N维单位阵；

4c)初始化中心参数α₀，即对Flevoland地区的极化SAR图像设定中心参数α₀＝2；对Oberpfaffenhofen地区的极化SAR图像设定中心参数α₀＝10；

4d)将当前的迭代次数设为1，设定最大迭代次数为20。

步骤5，估计极化SAR图像的先验参数s_ij和似然参数{μ_k,Σ₁}。

先验参数s_ij的估计方法一般有Derin方法、极化散射特征向量的欧氏距离方法等，但是Derin方法运算量较大，欧氏距离方法不能精确地描述极化SAR数据，因此本发明利用极化协方差矩阵的相似度估计极化SAR图像的先验参数s_ij，既能降低算法运算量，又能增强与极化SAR图像的数据相关性；

似然参数{μ_k,Σ₁}常采用最大似然方法进行估计，但是最大似然方法不能体现狄利克雷模型的随机性优势，因此本发明利用Gibbs采样算法对似然参数{μ_k,Σ₁}的后验分布进行采样，即可估计得到似然参数{μ_k,Σ₁}。

5a)从极化SAR图像中选取像素点i的极化协方差矩阵C_i与邻域像素点的极化协方差矩阵C_j，根据初始化的模型参数，计算C_i与C_j之间的对称的改进型Wishart距离d(C_i,C_j)：

$d(C_i,C_j)=\frac{1}{2}\left(Tr\right(C_j^{-1}{C}_i-q)+Tr(C_i^{-1}{C}_j-q\left)\right),$

其中Tr(·)表示求矩阵迹的操作；

5b)根据对称的改进型Wishart距离d(C_i,C_j)和初始化的模型参数，计算极化SAR图像的先验参数s_ij：

$s_{ij}=\exp \left(\frac{-d^2(C_i,C_j)}{2{\mathrm{\σ}}^2}\right),j\∈\∂\left(i\right);$

5c)根据原始的标记场X、观测场Y和初始化的模型参数，得到均值向量参数μ_k的后验概率p₁：

$p_1\∝\underset{X_i=k}{\Π}N\left(Y_i\right|{\mathrm{\μ}}_k,{\mathrm{\Σ}}_1)N\left({\mathrm{\μ}}_k\right|{\mathrm{\μ}}_0,{\mathrm{\Σ}}_0),$

其中N(·|·,·)表示正态分布；Π表示求各项乘积的操作；

5d)根据均值向量参数μ_k的后验概率p₁，利用Gibbs采样算法，估计得到极化SAR图像的均值向量参数μ_k；

5e)根据原始的标记场X、观测场Y和初始化的模型参数，得到方差参数的后验概率p₂：

$p_2\∝\mathrm{\Γ}(a+\frac{n}{2},b+\frac{1}{2}\underset{i=1,X_i=k}{\overset{n}{\Σ}}{\left(Y_{ir}-{\mathrm{\μ}}_{kr}\right)}^2),$

其中Γ(·,·)表示伽马分布；Y_ir表示极化SAR图像中像素点i的N维降噪的极化散射特征向量Y_i的第r个分量；μ_kr表示极化SAR图像中类别标记为k的像素点的均值向量参数μ_k的第r个分量；r＝1,2,…,N；

5f)根据方差参数的后验概率p₂，利用Gibbs采样算法，估计极化SAR图像的方差参数即可得到协方差矩阵参数

5g)由估计得到的均值向量参数μ_k和协方差矩阵参数Σ₁，得到似然参数{μ_k,Σ₁}。

步骤6，估计极化SAR图像新的标记场

马尔科夫随机场MRF模型在分类问题中根据极大化类别标记的后验概率不断地估计新的标记场直到达到最大迭代次数。因此本发明首先根据先验参数s_ij和似然参数{μ_k,Σ₁}得到类别标记的先验概率p₃和似然概率p₄，然后根据先验概率p₃和似然概率p₄得到后验概率p₅，通过极大化类别标记的后验概率p₅即可估计得到极化SAR图像新的标记场

6a)根据先验参数s_ij、原始的标记场X和初始化的模型参数，得到类别标记的先验概率p₃：

$p_3=\left\{\begin{array}{cc}\frac{\underset{l=1,l\≠i}{\overset{n}{\Σ}}\mathrm{\δ}(X_i,X_l)}{{\mathrm{\α}}_0+n-1}\exp (-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\Σ}{s}_{ij}\left(1-2\mathrm{\δ}\left(X_i,X_j\right)\right)),& X_i\∈X_{-i}\\ \frac{{\mathrm{\α}}_0}{{\mathrm{\α}}_0+n-1},& X_i\∉X_{-i}\end{array}\right.,$

其中X_j表示极化SAR图像中像素点j的类别标记；X_l表示极化SAR图像中像素点l的类别标记；表示极化SAR图像中像素点i的邻域系统；X_-i表示极化SAR图像中除了像素点i的类别标记以外的其他像素点的类别标记的集合；δ(X_i,X_l)表示Delta函数，当X_i与X_l相等时，δ(X_i,X_l)的值为1，否则为0；δ(X_i,X_j)表示Delta函数，当X_i与X_j相等时，δ(X_i,X_j)的值为1，否则为0；s_ij表示极化SAR图像的先验参数；exp(·)表示求指数操作；

6b)根据似然参数{μ_k,Σ₁}、原始的标记场X和观测场Y，得到类别标记的似然概率p₄：

$p_4=\frac{1}{{\left(2\mathrm{\π}\right)}^{N/2}|{\mathrm{\Σ}}_1{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_k\right)}^T{\mathrm{\Σ}}_1^{-1}(Y_i-{\mathrm{\μ}}_k\left)\right),X_i=k,$

其中μ_k表示极化SAR图像中类别标记为k的像素点的均值向量参数；Σ₁表示极化SAR图像的方差矩阵参数；(·)^-1表示求矩阵逆的操作；(·)^T表示求矩阵转置的操作；|·|表示求矩阵行列式的操作；X_i＝k表示极化SAR图像中像素点i的类别标记为k；

6c)基于Baysian准则，根据类别标记的先验概率p₃和似然概率p₄，得到类别标记的后验概率p₅；

$p_5=\left\{\begin{array}{cc}{z}^{-1}{p}_4\·n_k^{-i}\·\exp (-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\mathrm{\Σ}}{s}_{ij}(1-2\mathrm{\δ}(X_i,X_j)\left)\right),& k=1,2,...,C\\ z^{-1}{\mathrm{\α}}_0\frac{1}{{\left(2\mathrm{\π}\right)}^N|{\mathrm{\Σ}}_1{\mathrm{\Σ}}_0{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_0\right)}^T{\left({\mathrm{\Σ}}_1{\mathrm{\Σ}}_0\right)}^{-1}(Y_i-{\mathrm{\μ}}_0\left)\right),& k=C+1\end{array}\right.,$

其中指的是极化SAR图像中除了像素点i以外的类别标记为k的像素点的个数；z为归一化参数；C表示原始的标记场X中的类别数目；

6d)通过极大化类别标记的后验概率p₅，估计得到极化SAR图像中像素点i的新的类别标记

${\stackrel{\‾}{X}}_i=\underset{k=1,2,...,C,C+1}{max}\left\{p_5\right\};$

6e)将新的类别标记的集合确定为极化SAR图像的新的标记场

(7)判断是否达到最大迭代次数，输出分类结果：

将当前的迭代次数加1作为新的迭代次数，判定新的迭代次数是否达到步骤(4)所设定的最大迭代次数20，如果没有达到，则设定新的标记场为原始的标记场X，返回步骤(5)；如果达到，则将新的标记场确定为极化SAR图像的最终的分类结果并输出。

下面结合仿真实验对本发明的效果做进一步说明：

1)实验条件

实验仿真环境为：MATLAB R2010b,Intel(R)Core(TM)i7-4770CPU 3.4GHz，Window7旗舰版。

2)实验内容

实验1、分别应用基于MRF的极化SAR图像的监督分类方法，即Wanying Song等人提出的方法、Peter Orbanz等人提出的方法及本发明方法对Flevoland地区的极化SAR图像进行分类，结果如图2所示，其中：

图2(a)是Flevoland地区的极化SAR图像的Pauli分解图，

图2(b)是Wanying Song等人提出的方法得到的分类结果，

图2(c)是Peter Orbanz等人提出的方法得到的分类结果，

图2(d)是本发明方法得到的分类结果，

图2(e)是Flevoland地区的极化SAR图像的分类结果参考图，

图2(f)是Flevoland地区的极化SAR图像各个类别对应的地物类型。

从图2可以看到，本发明在抗噪性能和边界信息保持上均优于Wanying Song等人提出的方法和Peter Orbanz等人提出的方法，而且相比于Peter Orbanz等人提出的方法，本发明得到的分类结果中误分类数目明显减少。

为了验证分类结果的质量，以分类精度为性能指标参数，评价实验1中本发明和现有方法得到的分类结果，以图2(e)为分类结果参考图。本发明和现有方法得到的分类结果的分类精度(％)如表2所示。

表2

从表2可以看出，本发明方法得到的分类结果中各个类别的分类精度和总精度都比较高，虽然Wanying Song等人提出的方法也能获得较高的分类精度，但是该方法最大的弊端是不能自动地确定极化SAR图像中的类别数目，而Peter Orbanz等人提出的方法由于缺乏与极化SAR数据的相关性，误分类数目较多，因此分类精度较低。

实验2、分别应用基于MRF的极化SAR图像的监督分类方法，即Wanying Song等人提出的方法、Peter Orbanz等人提出的方法及本发明方法对Oberpfaffenhofen地区的极化SAR图像进行分类，结果如图3所示，其中：

图3(a)是Oberpfaffenhofen地区的极化SAR图像的Pauli分解图，

图3(b)是Wanying Song等人提出的方法得到的分类结果，

图3(c)是Peter Orbanz等人提出的方法得到的分类结果，

图3(d)是本发明方法得到的分类结果。

从图3可以看到，相比于Wanying Song等人提出的方法和Peter Orbanz等人提出的方法，本发明方法得到的分类结果中受噪声影响较小，而且细节信息的保持更加完整，最重要的是能够对相似的地物目标进行区分，误分类区域较少。

综上，本发明能使同质区域分类结果更平滑且能更好地保持边缘信息，分类精度明显提高。

Claims (5)

1.一种基于狄利克雷MRF混合模型的极化SAR图像分类方法，包括如下步骤：

(1)输入极化SAR图像；

(2)从极化SAR图像中提取N个极化散射特征，并对其进行归一化，得到N个归一化的极化散射特征，建立归一化的极化散射特征空间F₁；

(3)以极化SAR图像中的每个像素点为中心像素点，选取大小为3×3的窗口，使用广义均值方法，对归一化的极化散射特征空间F₁逐点进行降噪处理，得到极化散射特征空间F₂；

(4)初始化模型参数

将极化SAR图像的每个像素点的类别标记的集合设定为原始的标记场X＝{X_i}_{i＝1,2,…,n}，其中X_i表示极化SAR图像中像素点i的类别标记，n表示极化SAR图像的像素点个数；初始化极化SAR图像的每个像素点的类别标记为1，即X_i＝1；

用极化散射特征空间F₂作为观测场Y＝{Y_i}_{i＝1,2,…,n}，其中Y_i表示极化SAR图像中像素点i的N维降噪的极化散射特征向量；

初始化形状参数a＝1，尺度参数b＝5，核宽参数σ＝1，平衡参数λ＝1，维度参数q＝3，第一正态参数μ₀初始化为N维零向量，第二正态参数Σ₀初始化为N维单位阵；

初始化中心参数α₀，即对Flevoland地区的极化SAR图像设定中心参数α₀＝2；对Oberpfaffenhofen地区的极化SAR图像设定中心参数α₀＝10；

将当前的迭代次数设为1，设定最大迭代次数为20；

(5)根据原始的标记场X、观测场Y和初始化的模型参数，估计极化SAR图像的先验参数s_ij和似然参数{μ_k,Σ₁}；

(6)估计极化SAR图像新的标记场

6a)根据先验参数s_ij、原始的标记场X和初始化的模型参数，得到类别标记的先验概率p₃：

$p_3=\left\{\begin{array}{cc}\frac{\underset{l=1,l\≠i}{\overset{n}{\mathrm{\Σ}}}\mathrm{\δ}\left(X_i,X_l\right)}{{\mathrm{\α}}_0+n-1}\exp \left(-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\mathrm{\Σ}}{s}_{ij}\left(1-2\mathrm{\δ}\left(X_i,X_j\right)\right)\right),& X_i\∈X_{-i}\\ \frac{{\mathrm{\α}}_0}{{\mathrm{\α}}_0+n-1},& X_i\∉X_{-i}\end{array}\right.,$

其中X_j表示极化SAR图像中像素点j的类别标记；X_l表示极化SAR图像中像素点l的类别标记；表示极化SAR图像中像素点i的邻域系统；X_-i表示极化SAR图像中除了像素点i的类别标记以外的其他像素点的类别标记的集合；δ(X_i,X_l)表示Delta函数，当X_i与X_l相等时，δ(X_i,X_l)的值为1，否则为0；δ(X_i,X_j)表示Delta函数，当X_i与X_j相等时，δ(X_i,X_j)的值为1，否则为0；s_ij表示极化SAR图像的先验参数；exp(·)表示求指数操作；

6b)根据似然参数{μ_k,Σ₁}、原始的标记场X和观测场Y，得到类别标记的似然概率p₄：

$p_4=\frac{1}{{\left(2\mathrm{\π}\right)}^{N/2}|{\mathrm{\Σ}}_1{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_k\right)}^T{\mathrm{\Σ}}_1^{-1}(Y_i-{\mathrm{\μ}}_k\left)\right),X_i=k,$

其中μ_k表示极化SAR图像中类别标记为k的像素点的均值向量参数；Σ₁表示极化SAR图像的协方差矩阵参数；(·)^-1表示求矩阵逆的操作；(·)^T表示求矩阵转置的操作；|·|表示求矩阵行列式的操作；X_i＝k表示极化SAR图像中像素点i的类别标记为k；

6c)根据类别标记的先验概率p₃和似然概率p₄，得到类别标记的后验概率p₅；

$p_5=\left\{\begin{array}{cc}{z}^{-1}{p}_4\·n_k^{-i}\·\exp (-\mathrm{\λ}\underset{j\∈\∂\left(i\right)}{\mathrm{\Σ}}{s}_{ij}(1-2\mathrm{\δ}\left(X_i,X_j\right)\left)\right),& k=1,2,...,C\\ z^{-1}{\mathrm{\α}}_0\frac{1}{{\left(2\mathrm{\π}\right)}^N|{\mathrm{\Σ}}_1{\mathrm{\Σ}}_0{|}^{1/2}}\exp (-\frac{1}{2}{\left(Y_i-{\mathrm{\μ}}_0\right)}^T{\left({\mathrm{\Σ}}_1{\mathrm{\Σ}}_0\right)}^{-1}(Y_i-{\mathrm{\μ}}_0\left)\right),& k=C+1\end{array}\right.,$

其中指的是极化SAR图像中除了像素点i以外的类别标记为k的像素点的个数；z为归一化参数；C表示原始的标记场X中的类别数目；

6d)通过极大化类别标记的后验概率p₅，估计得到极化SAR图像中像素点i的新的类别标记

${\stackrel{\‾}{X}}_i=\underset{k=1,2,...,C,C+1}{max}\left\{p_5\right\};$

6e)将新的类别标记的集合确定为极化SAR图像的新的标记场

(7)判断是否达到最大迭代次数，输出分类结果：

将当前的迭代次数加1作为新的迭代次数，判定新的迭代次数是否达到步骤(4)所设定的最大迭代次数20，如果没有达到，则设定新的标记场为原始的标记场X，返回步骤(5)；如果达到，则将新的标记场确定为极化SAR图像的最终的分类结果并输出。

2.根据权利要求1所述的基于狄利克雷MRF混合模型的极化SAR图像分类方法，其特征在于，所述步骤(2)的实现如下：

2a)从极化SAR图像中提取N个极化散射特征F_r，r＝1,2,...,N，N＝17表示极化散射特征的个数，其中：

F₁，F₂，F₃，F₄分别表示用结合层次后向策略的Yamaguchi-Sato分解方法对极化SAR图像进行分解得到的螺旋散射功率、体散射功率、二面角散射功率、表面散射功率；

F₅＝T₁₁，T₁₁表示极化SAR图像的相干矩阵T的第1行第1列的元素；

F₆＝|T₁₂|，T₁₂表示极化SAR图像的相干矩阵T的第1行第2列的元素；

F₇＝|T₁₃|，T₁₃表示极化SAR图像的相干矩阵T的第1行第3列的元素；

F₈＝T₂₂，T₂₂表示极化SAR图像的相干矩阵T的第2行第2列的元素；

F₉＝|T₂₃|，T₂₃表示极化SAR图像的相干矩阵T的第2行第3列的元素；

F₁₀＝T₃₃，T₃₃表示极化SAR图像的相干矩阵T的第3行第3列的元素；

F₁₁，F₁₂，F₁₃分别表示极化SAR图像的相干矩阵T的3个特征值；

F₁₄表示极化SAR图像的相干矩阵T的对角线元素和，即极化总功率；

F₁₅，F₁₆，F₁₇分别表示用Cloude分解方法对极化SAR图像进行分解得到的极化散射特征，即F₁₅表示极化熵，F₁₆表示平均散射角，F₁₇表示反熵；

2b)通过下式分别对极化SAR图像的N个极化散射特征F_r，r＝1,2,...,N，进行归一化，得到N个归一化的极化散射特征

$\stackrel{\‾}{F_r}=\frac{F_r-min\left(F_r\right)}{\max \left(F_r\right)-min\left(F_r\right)};$

2c)将归一化的极化散射特征中像素点i的特征值的集合确定为极化SAR图像中像素点i的归一化的极化散射特征向量F_i，i＝1,2,...,n：

$F_i={\left\{\stackrel{\‾}{F_r}\left(i\right)\right\}}_{r=1,2,...,N};$

2d)将归一化的极化散射特征向量F_i的集合确定为归一化的极化散射特征空间F₁：

F₁＝{F_i}_{i＝1,2,...,n}。

3.根据权利要求1所述的基于狄利克雷MRF混合模型的极化SAR图像分类方法，其特征在于，步骤(3)中对归一化的极化散射特征空间F₁逐点进行降噪处理，按如下步骤进行：

3a)根据归一化的极化散射特征提取极化SAR图像中像素点i的第r个降噪的极化散射特征G_r(i)：

$G_r\left(i\right)={\left(\underset{t\∈N_{3\×3}\left(i\right)}{\Π}\stackrel{\‾}{F_r}\right(t\left)\right)}^{\frac{1}{9}},$

其中表示第r个归一化的极化散射特征中像素点t的特征值；N_3×3(i)表示以像素点i为中心像素点的3×3的窗口；

3b)将N个降噪的极化散射特征G_r(i)的集合确定为极化SAR图像中像素点i的N维降噪的极化散射特征向量G_i：

G_i＝{G_r(i)}_{r＝1,2,...,N}；

3c)将N维降噪的极化散射特征向量G_i的集合确定为极化散射特征空间F₂：

F₂＝{G_i}_{i＝1,2,...,n}。

4.根据权利要求1所述的基于狄利克雷MRF混合模型的极化SAR图像分类方法，其特征在于，步骤(5)中估计极化SAR图像的先验参数s_ij，按如下步骤进行：

5a)从极化SAR图像中选取像素点i的极化协方差矩阵C_i与邻域像素点的极化协方差矩阵C_j，根据初始化的模型参数，计算C_i与C_j之间的对称的改进型Wishart距离d(C_i,C_j)：

$d(C_i,C_j)=\frac{1}{2}\left(Tr\right(C_j^{-1}{C}_i-q)+Tr(C_i^{-1}{C}_j-q\left)\right),$

其中Tr(·)表示求矩阵迹的操作；

5b)根据对称的改进型Wishart距离d(C_i,C_j)和初始化的模型参数，计算极化SAR图像的先验参数s_ij：

$s_{ij}=\exp \left(\frac{-d^2(C_i,C_j)}{2{\mathrm{\σ}}^2}\right),j\∈\∂\left(i\right).$

5.根据权利要求1所述的基于狄利克雷MRF混合模型的极化SAR图像分类方法，其特征在于，步骤(5)中估计极化SAR图像的似然参数{μ_k,Σ₁}，按如下步骤进行：

5c)根据原始的标记场X、观测场Y和初始化的模型参数，得到均值向量参数μ_k的后验概率p₁：

$p_1\∝\underset{X_i=k}{\Π}N\left(Y_i\right|{\mathrm{\μ}}_k,{\mathrm{\Σ}}_1)N\left({\mathrm{\μ}}_k\right|{\mathrm{\μ}}_0,{\mathrm{\Σ}}_0),$

其中N(·|·,·)表示正态分布；Π表示求各项乘积的操作；

5d)根据均值向量参数μ_k的后验概率p₁，利用Gibbs采样算法，估计得到极化SAR图像的均值向量参数μ_k；

5e)根据原始的标记场X、观测场Y和初始化的模型参数，得到方差参数的后验概率p₂：

$p_2\∝\mathrm{\Γ}(a+\frac{n}{2},b+\frac{1}{2}\underset{i=1,X_i=k}{\overset{n}{\Σ}}{\left(Y_{ir}-{\mathrm{\μ}}_{kr}\right)}^2),$

其中Γ(·,·)表示伽马分布；Y_ir表示极化SAR图像中像素点i的N维降噪的极化散射特征向量Y_i的第r个分量；μ_kr表示极化SAR图像中类别标记为k的像素点的均值向量参数μ_k的第r个分量；r＝1,2,…,N；

5f)根据方差参数的后验概率p₂，利用Gibbs采样算法，估计极化SAR图像的方差参数即可得到协方差矩阵参数

5g)由估计得到的均值向量参数μ_k和协方差矩阵参数Σ₁，得到似然参数{μ_k,Σ₁}。