CN105321163A - 检测全极化sar图像的变化区域的方法和装置 - Google Patents

Abstract

本发明公开一种检测全极化SAR图像的变化区域的方法和装置，该方法包括：使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到所述主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据所述主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，该两幅图像为包括同一对象区域的高分辨率全极化SAR图像；根据该两幅图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比该两幅图像得到差异信息，基于差异信息确定相对于该两幅图像的差异图像；将差异图像与该两幅图像进行比对，以检测出变化区域。

Description

检测全极化SAR图像的变化区域的方法和装置

技术领域

本发明涉及一种检测全极化SAR图像的变化区域的方法和装置。

背景技术

进入21世纪以来，世界人口急剧增加，土地沙漠化逐渐严重，资源日益匮乏，自然灾害频发。2008年汶川大地震、2010年青海玉树大地震及灾难性印尼海啸、2011年日本福岛大地震和2013年的台风“海燕”多次唤起了全人类对遥感在灾害监测与评估方面的重视。在实际应用中，上述大多数监测问题都是以变化检测的形式提出，即通过不同时段的观测，来识别城市、森林、水域、农田等地物状态的变化过程。雷达对地观测技术是解决上述问题的重要手段之一，合成孔径雷达(SyntheticApertureRadar,SAR)作为现代遥感领域的一项突破性成就，其使雷达的基本功能发生了重大改变。近些年，SAR凭借其相较于光学遥感影像，在全天时、全天候，尤其是多云多雨天气状况下也能有效施行对地观测的独特优势，成为了目前遥感领域中的一种有效态监测手段，其被越来越多地应用到如震后灾害评估、海啸受灾状况评估、农作物生长状况监测等的变化检测工作中。在城镇化建设上，SAR图像变化检测技术作为研究城区变化的重要手段，已在城区扩建、非法建筑及土地滥用、湿地生态面积缩减状况、城中耕地面积监测和湖波变化中发挥重要作用。

但是，由于受到数据源和SAR技术的限制，目前已提出的SAR影像变化检测方法基本上都是针对单通道(单一微波波段、单一极化方式)的SAR影像设计的，而很少引入极化信息。2007年，加拿大Radarsat-2、德国TerrraSAR-X和意大利Cosmo-SkyMed新型雷达系统的出现，引起了全极化雷达遥感研究和应用的热潮。与传统的单极化SAR系统相比，极化SAR(PolarimetricSyntheticApertureRadar,PolSAR)能在不同收发极化组合下，测量地物目标的极化散射特性(即，极化SAR能够通过发射和接收不同极化方式的电磁波，来探测地面目标对电磁波的调制特性)，能够获取同一地物目标更全面信息，从而使其在目标检测、识别、分类以及目标参数反演等方面也具有显著的优势。受雷达系统侧视成像和斑噪特性的影响，通过单极化的后向散射系数很难对照射地物进行全面解译和识别。利用多极化信息可以得到相互补充的不同信息，增加信息含量，获取同一地物目标更全面信息，有利于地物变化的检测与识别。

最近几年，SAR成像技术日益成熟，图像质量逐渐提高，分辨率不断增强，获取数据的能力和精度越来越高。相对于高分辨率光学影像、中低分辨率SAR影像，高分辨率SAR影像特点主要表现在：1)高空间分辨率使得SAR图像纹理结构更加清晰，能够在较小的空间尺度上探测细节变化，实现对建筑物、道路等地面目标细微特征的探测；2)地物空间拓扑信息更加丰富，更充分地反应在了影像地物的尺寸、形状、邻域地物关系。然而，对于极化SAR来说，米级分辨率影像在带来地物丰富细节的同时，也带来了数据存储量倍增、独立地物成像特征破碎、统计特征复杂化等诸多高分辨率SAR影像的解译难题。这时，传统的基于纹理均匀假设的极化SAR杂波统计模型和相应的相似性度量，已不再适用于高分辨率全极化SAR图像分割、分类和目标检测等应用需求。目前，全极化SAR变化检测还存在下面2个不足：

①在PolSAR杂波模型方面，单视情况下描述全极化散射矢量的球不变随机矢量模型(SphericallyInvariantRandomVectors,SIRV)模型已基本上与单极化SAR模型的发展一一对应。但是，为了方便SAR数据压缩、初步去噪和局部散射体分析，全极化SAR数据常常以多视方式存储。多视后的SAR数据一般用极化协方差矩阵或相干矩阵表征。目前用于拟合极化相干矩形的多视极化SAR乘性模型中，往往假设模型中纹理的均值为1，这种假设在高分辨率图像上是不成立的。对于高分辨率数据而言，其纹理结构十分清晰，包含了极化数据的功率信息，纹理的特征已成为地物识别和检测的重要信息。一般情况下，除非对相干矩阵进行归一化，否则其纹理均值是不等于1的。

②在全极化SAR变化测度研究上，现在主要的测度Wishart距离和SIRV距离中均没有包含纹理的变化信息。对于SAR变化检测而言，图像之间最明显的变化一般就是强度的变化，这个变化反映在全极化SAR图像上就是图像功率的变化，而图像功率又是由纹理表征的。因此，纹理的变化是全极化SAR地物类型变化的主要标志之一，应该是变化检测算法重点检测的信息。在定义高分辨率全极化SAR变化检测时，既要考虑地物极化散射特性(斑噪)的变化，也要同时考虑地物本身散射功率(纹理)的变化。

发明内容

本发明所要解决的问题是现有的全极化SAR模型不能提供精确的区域变化区域，提供一种基于极不均匀纹理的多视极化SAR乘性模型能检测出准确的变化区域。

为了解决上述问题，本发明提供了一种检测全极化SAR图像的变化区域的方法，方法包括：步骤S2：使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，主、辅SAR图像为包括同一对象区域的高分辨率全极化SAR图像；步骤S4：根据主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比主、辅SAR图像得到差异信息，基于差异信息确定相对于主、辅SAR图像的差异图像；步骤S6：将差异图像与主SAR图像和/或辅SAR图像进行比对，以检测出变化区域。

作为优选，在步骤S2之前，方法进一步包括：步骤S1：对接收的主、辅SAR图像进行同一对象区域的对应匹配。

作为优选，步骤S2中的使用基于极不均匀纹理的多视极化SAR乘性模型对主、辅SAR图像进行拟合包括：使用基于极不均匀纹理的乘性模型的概率密度函数推导主、辅SAR图像的纹理图像、斑噪协方差图像。

基于极不均匀纹理的乘性模型的概率密度函数为：

$p_W(W;L,\mathrm{\τ},\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{\∫}_0^{\∞}\frac{1}{{\mathrm{\τ}}^{\mathrm{pL}}}\exp \{-\mathrm{tr}\left(\frac{{\mathrm{\Σ}}^{-1}W}{\mathrm{\τ}}\right)\}\·p_{\mathrm{\τ}}\left(\mathrm{\τ}\right)\mathrm{d\τ}$

其中，L代表多视为L视，W为L视下未标准化的全极化相干矩阵，Pw为W的概率密度函数，τ表示图像的纹理，Σ为斑噪协方差矩阵，Γ_p(L)为L视下概率密度函数的复数形式的多元Gamma函数。

作为优选，步骤S2中的分别根据主、辅SAR图像的纹理图像确定其对应的纹理图像的概率密度函数的分布参数包括：设置主、辅SAR图像的纹理图像以符合Gamma函数分布。

作为优选，基于极不均匀纹理的多视极化SAR乘性模型的相似性度量为：

S_P＝MLL_P(X)+MLL_P(Y)-MLL_P(X∪Y)

其中，Sp为Pw的变化测度，X为主SAR图像内的子区域，MLLp(X)表示区域X内Pw的最大对数似然函数，Y为辅SAR图像内的子区域，MLLp(Y)表示区域Y内Pw的最大对数似然函数，区域X与区域Y为对应子区域。

作为优选，在步骤S6之前，包括：步骤S5：通过对差异图像进行阀值分割来确定其区域和边界。

作为优选，步骤S5包括：S51、对差异图像进行自适应KI阀值分割得到区别变化区域与不变区域的作为二值化图像的二值化掩膜；S52、通过连通区搜索算法去除二值化图像中的变化面积小于预定值的变化区域；S53、通过形态学的闭运算中的膨胀操作对变化区域进行膨胀得到膨胀区域，并对该膨胀区域的周界进行记录；S54、通过形态学的闭运算中的腐蚀操作，采用与膨胀操作相对应的比例对膨胀区域进行腐蚀，得到腐蚀后区域；S55、将腐蚀后区域的周界与膨胀区域的周界进行作差得到差异图像中变化区域的边界。

作为优选，在步骤S6之后进一步包括：步骤S7：将变化区域与主、辅SAR图像进行叠合，确定主、辅SAR图像内的变化部分。

本发明还提供了一种检测全极化SAR图像的变化区域的装置，其特征在于，全极化SAR图像基于极不均匀纹理，检测装置包括：参数确定模块，其配置为使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，主、辅SAR图像为包括同一对象区域的高分辨率全极化SAR图像；差异信息确定模块，其配置为根据主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比主、辅SAR图像得到差异信息，基于差异信息确定相对于主、辅SAR图像的差异图像；比对模块，其配置为将差异图像与主SAR图像和/或辅SAR图像进行比对，以检测出变化区域。

作为优选，上述装置进一步包括：区域匹配模块，其配置为对接收的主、辅SAR图像进行同一对象区域的对应匹配。

本发明相对于现有技术的有益效果在于：

1、基于极不均匀纹理的多视极化SAR乘性模型检测出全极化SAR图像的变化区域，相对于现有技术具有高的检测率和更低的虚警率，实现PolSAR影像的高可信检测，能够进一步推动SAR成像新理论、新方法的发展和应用；

2、确定一种多视PolSAR乘性模型，能够充分考虑地物纹理不均匀分布，能够实现对高分辨率全极化SAR图像的精确拟合；

3、确定一种对应于上述多视PolSAR乘性模型的全极化SAR变化测度，该变化测度中同时考虑了纹理和斑噪的变化；

4、相对于基于Wishart距离的变化检测方法，根据本发明技术方案的方法具有更高的检测率和更低的虚警率。

针对于传统极化SAR乘性模型的不足(即，现有技术中多视极化SAR乘性模型都是基于纹理均值为1的假设推导的)，本文提出一种极不均匀纹理的多视的极化SAR乘性模型，并推导了模型中相应的斑噪协方差矩阵的估计方法，该模型中充分考虑地物纹理不均匀分布，能够实现对高分辨率全极化SAR图像的精确拟合。根据本发明实施例的多视极化SAR乘性模型，基于相似性度量的方法，提出了一个相应的全极化SAR变化测度，在这一变化测度中，同时包含了纹理信息和斑噪协方差矩阵的变化信息，这更有利于提取更为全面的全极化SAR图像的变化信息，在此测度上，提出了一种基于极不均匀纹理模型的全极化SAR变化检测方法，实现PolSAR影像的高可信检测，能够进一步推动SAR成像新理论、新方法的发展和应用。

附图说明

图1为根据本发明实施例的方法的流程图；

图2为如图2所示为根据本发明的一个实施例的方法的流程图；

图3a为根据本发明实施例的某城区的影像在2009年4月9日的Pauli合成图；

图3b为图3a中的城区的影像在2010年6月15日的Pauli合成图；

图4a为根据本发明实施例的从主图像的像素斑噪协方差矩阵图像内任意选取的像素1进行迭代算法的收敛情况；

图4b为根据本发明实施例的从辅图像的像素斑噪协方差矩阵图像内任意选取的像素2进行迭代算法的收敛情况；

图5a为2009年得到的主图像通过迭代算法得到的纹理图像的示意图；

图5b为2010年得到的辅图像通过迭代算法得到的纹理图像的示意图；

图6a为基于本发明实施例的相似性度量得到的主、辅图像的差异图；

图6b为基于现有技术中的Wishart距离得到的主、辅图像的差异图；

图7为图6a和图6b所示的两个差异图的ROC曲线的示意图；

图8a为对图6a所示差异图的二值化掩膜进行KI阈值分割后的示意图；

图8b为对图8a所示图进行连通区搜索后的示意图；

图9为根据本发明的实施例对图8b所示图进行形态学运算后的示意图；

图10a为将图9变化区域与2009年该地区分类图进行融合后的示意图；

图10b为图10a所示的地区在2010年的分类示意图；

图10c为图10a和图10b所示图中的各部分所分别代表的地理形态的示意图。

具体实施方式

以下结合附图对本发明的进行详细描述。

根据本发明的实施例，提供了一种全极化SAR图像的变化区域的方法。其中，全极化SAR图像的每个像素点均具有独立的纹理参数和斑噪参数，在本发明中每个像素点的纹理参数以τ表示，τ为正随机变量，一幅全极化SAR图像的所有像素点的纹理参数τ可组成对应于该全极化SAR图像的纹理图像，该纹理图像的概率密度函数的分布参数以P_τ(τ)表示，在已知纹理图像的情况下，可以通过假设纹理图像的概率密度函数从而得出其分布参数P_τ(τ)，而在本发明中每个像素点的斑噪参数以斑噪协方差函数Σ表示，即一幅全极化SAR图像的所有像素点的斑噪参数Σ也可组成对应于该全极化SAR图像的斑噪协方差图像。

此外，通常将全极化SAR图像存储为L视下未标准化的全极化W或者T，由于在多视情况下T＝W/L，因而T与W的地位相同，其中，L代表视数。

如图1所示为根据本发明实施例的方法，该方法可以包括：

步骤S2：使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到主、辅SAR图像的纹理图像τ、斑噪协方差图像Σ，并分别根据主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数P_τ(τ)，主、辅SAR图像为包括同一对象区域的高分辨率极化SAR图像，其中，本发明的目标在于确定主、辅SAR图像的变化，因而可以根据用户自行需求设置文中“先后”的时间间隔；

步骤S4：根据主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比主、辅SAR图像的差异图像得到差异信息，基于差异信息确定相对于主、辅SAR图像的差异图像；

步骤S6：将差异图像与主SAR图像和/或辅SAR图像进行比对，以检测出变化区域。

通过该方法可以基于极不均匀纹理的多视极化SAR乘性模型检测出全极化SAR图像的变化区域，因而相对于基于均匀纹理的现有技术的检测方案来说具有更高的检测率和更低的虚警率。

并且，在步骤S2之前，根据本发明实施例的方法可以进一步包括：

步骤S1：对接收的主、辅SAR图像进行同一对象区域的对应匹配，即对于主、辅SAR图像中的相应区域进行匹配使得主辅SAR图像内的各点能够一一对应，从而使得后续比对操作时基于对应区域进行。步骤S1属于对主、辅SAR图像进行预处理，由于SAR图像容易受到斑噪的影响，对图像进行滤波处理会有利于图像的解译。但是，根据本发明的实施例所采用的多视极化SAR乘性模型(在实际应用中也可称为全极化SAR图像杂波模型)正是针对于SAR图像内斑噪的形成机理建立的，对SAR图像进行去噪反而会影响到杂波模型的拟合精度。因此，本发明的预处理在一般情况下可以省略图像滤波的步骤，除非图像的等效视数不够大，其中，等效视数为图像的一个参数，参考本领域技术人员的公知即可定义“等效视数不够大”的情况。

优选地，步骤S2中的使用基于极不均匀纹理的多视极化SAR乘性模型对主、辅SAR图像进行拟合包括：使用基于极不均匀纹理的乘性模型的概率密度函数推导主、辅SAR图像的纹理图像、斑噪协方差图像。

基于极不均匀纹理的乘性模型的概率密度函数为：

$p_W(W;L,\mathrm{\τ},\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{\∫}_0^{\∞}\frac{1}{{\mathrm{\τ}}^{\mathrm{pL}}}\exp \{-\mathrm{tr}\left(\frac{{\mathrm{\Σ}}^{-1}W}{\mathrm{\τ}}\right)\}\·p_{\mathrm{\τ}}\left(\mathrm{\τ}\right)\mathrm{d\τ}$

其中，L代表多视为L视，W为L视下未标准化的极化相干矩阵，Pw为W的概率密度函数，其推导方法在下文中具体描述，τ表示图像的纹理(也可称为主、辅SAR图像内一含像素的纹理参数)，Σ为斑噪协方差矩阵，Γ_p(L)为L视下概率密度函数的复数形式的多元Gamma函数。

进一步优选地，步骤S2中的分别根据主、辅SAR图像的纹理图像确定其对应的纹理图像的概率密度函数的分布参数可以包括：

设置主、辅SAR图像的纹理图像以符合Gamma函数分布，由于纹理是不均匀的，用户可以根据需求设置纹理图像符合本领域公知技术中的其他函数分布，Gamma函数分布仅作为一种优选。

此外，基于极不均匀纹理的多视极化SAR乘性模型的相似性度量可以为：

S_P＝MLL_P(X)+MLL_P(Y)-MLL_P(X∪Y)

其中，Sp为Pw的变化测度，X为主SAR图像内的子区域，MLLp(X)表示区域X内Pw的最大对数似然函数，Y为辅SAR图像内的子区域，MLLp(Y)表示区域Y内Pw的最大对数似然函数，区域X与区域Y为主、辅SAR图像内的对应区域。其中，在X与Y均未变化的特殊情况下，Sp＝0。该Sp的推导方法在下文中具体描述，该Sp中同时考虑了极不均匀纹理和斑噪的变化，区别于现有技术中将纹理视为均匀的情况，当然，现有技术中纹理视为均匀的情况为本发明所考虑的极不均匀纹理的一种特殊情况，因而本发明也可适用于纹理均匀的图像处理情况。

此外，在步骤S6之前，根据本发明实施例的方法还可以包括：步骤S5，即通过对差异图像进行阀值分割来确定其区域和边界，具体地，步骤S5包括：

S51、对差异图像进行自适应KI阀值分割得到区别变化区域与不变区域的作为二值化图像的二值化图像(也可以称为二值化掩膜图像)，在实际应用中，变化区域与不变区域通常在二值化图像中以两种不同的颜色或其他有视觉差异的方式进行表达，；

S52、通过连通区搜索算法去除所述二值化图像中的变化面积(实际单块面积)小于预定值(优选的，预定值为100平方米)的变化区域，在实际应用根据本发明的方法时，SAR通常为以米级为精度单位，因而根据SAR图像与实际面积的比例设置，可得到大小不同的图像，本文中的“实际单块面积”指SAR采集的真实地表的面积。在SAR图像中，二值化图像中100平方米的变化区域通常仅在肉眼中示出为一点，消除这些点使得变化区域变得连续平滑是为了后续操作的需求服务；

S53、通过形态学的闭运算中的膨胀操作对变化区域进行膨胀得到膨胀区域，并对该膨胀区域的周界进行记录；

S54、通过形态学的闭运算中的腐蚀操作，采用与所述膨胀操作相对应的比例对所述膨胀区域进行腐蚀，得到腐蚀后区域；

S55、将腐蚀后区域的周界与膨胀区域的周界进行作差得到差异图像中所述变化区域的边界。

此外，在步骤S6之后可以进一步包括步骤S7：将变化区域与主、辅SAR图像进行叠合，确定主、辅SAR图像内的变化部分，即通过两个图像的叠合可以在像素点的精度里判定哪些范围有变化，而哪些范围没有变化。并且，可以将得到的最终的变化区域与原来的主辅图像进行叠合，输出土地利用图。

根据本发明的另一个实施例，还提供一种检测全极化SAR图像的变化区域的装置，全极化SAR图像基于极不均匀纹理，检测装置包括：

参数确定模块，其配置为使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，主、辅SAR图像为包括同一对象区域的高分辨率全极化SAR图像；

差异信息确定模块，其配置为根据主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比主、辅SAR图像得到差异信息，基于差异信息确定相对于主、辅SAR图像的差异图像；

比对模块，其配置为将差异图像与主SAR图像和/或辅SAR图像进行比对，以检测出变化区域。

作为优选，上述装置进一步包括：

区域匹配模块，其配置为对主、辅SAR图像进行同一区域的对应匹配。

如图2所示为根据本发明的一个实施例的方法，具体包括：

步骤S21：预处理，即对采集的主、辅SAR图像的数据进行图像配准；

步骤S22：进行差异图像提取，即利用根据本发明实施例的改进的多视极化SAR乘性模型拟合极化主、辅SAR图像，并基于根据本发明的相似性试题进行差异图提取；

步骤S23：阈值分割，即对提取的差异图进行自适应的KI阀值分割，并得到二值化掩膜；

步骤S24：图像整合，即将得到的二值化掩膜与主和/或辅SAR图像进行图像整合，然后处理输出。

以下对极化SAR乘性模型进行定义，以说明本发明中使用的基于极不均匀纹理的乘性模型的概率密度函数的原理：

单视(即，L＝1)情况下，SIRV模型(球不变随机矢量模型)认为全极化数据的散射矢量k是由复圆高斯矢量z与一个正随机参数τ的平方根相乘得到的，因而可以称其为乘性模型。此处τ表示图像的纹理参数，z表征图像的斑噪。根据这一定义，可以把L视下未标准化的极化相干矩阵W表示为

$W=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{k}_l{k}_l^H={\mathrm{KK}}^H---\left(1\right)$

其中

式中，矩阵R是对角线元素为图像纹理的是实对角矩阵，H代表共轭转秩，矩阵Z^H是由斑噪矢量组成的均值为0多元复高斯矩阵，记为其PDF(概率密度函数)为

$p_{Z^H}(Z^H;L,\mathrm{\Σ})={\mathrm{\π}}^{-\mathrm{Lp}}{\left|\mathrm{\Σ}\right|}^{-L}\mathrm{etr}(-{\mathrm{\Σ}}^{-1}{\mathrm{ZZ}}^H)---\left(3\right)$

为了推导矩阵W的概率密度函数，首先假设矩阵R是给定的常数矩阵，即假设图像纹理是已知的。定义那么矩阵A到矩阵Z^H的雅可比行列式为

$J(A\→Z^H)=\underset{i=1}{\overset{L}{\mathrm{\Π}}}{\sqrt{{\mathrm{\τ}}_i}}^{2p}={\left|R\right|}^p---\left(4\right)$

这样，可以得到矩阵R给定情况下，矩阵A的PDF为

$p_A\left(A\right|R;L,\mathrm{\Σ})=\frac{1}{{\mathrm{\π}}^{\mathrm{Lp}}{\left|\mathrm{\Σ}\right|}^L{\left|R\right|}^p}\exp \{-\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{A}^H{R}^{-1}A\right)\}---\left(5\right)$

从式(1)可以得到W＝A^HA，因此矩阵A到矩阵W的雅可比行列式为

$J(A\→W)=\frac{{\mathrm{\π}}^{\mathrm{Lp}}{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right)}---\left(6\right)$

此处Γ_p(L)为复数形式的多元Gamma函数。

这样，联合式(6)和式(5)，就可以得到在矩阵R给定情况下，未标准化的极化相干矩阵W的PDF为

$p_W\left(W\right|R;L,\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|R\right|}^p{\left|\mathrm{\Σ}\right|}^L}\exp \{-\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{\mathrm{KR}}^{-1}{K}^H\right)\}---\left(7\right)$

上式p_W(W|R；L,Σ)对矩阵R在从0到∞进行积分，就可以得到W的PDF为

$p_W(W;L,R,\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{\∫}_0^{\∞}\frac{1}{{\left|R\right|}^p}\exp \{-\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{\mathrm{KR}}^{-1}{K}^H\right)\}{p}_R\left(R\right)\mathrm{dR}---\left(8\right)$

式(8)就是多视极化SAR乘性模型的表达式，或者叫做SIRV模型到多视情况的扩展，记为W∈S(L,R,Σ)。当图像纹理为Diracdelta函数时，即p_R(R)＝δ(R-I_L)，其中矩阵I_L为L×L的单位矩阵。这时图像中每个像素的纹理相等，多视极化SAR乘性模型就退化为复Wishart分布模型。

$p_W(W;L,\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}\mathrm{etr}(-{\mathrm{\Σ}}^{-1}W)---\left(9\right)$

记为W∈W(L,Σ)。这说明复Wishart分布模型可以看出多视极化SAR乘性模型的一个特例。

在式(8)中，由于纹理τ的分布是不确定的，因此矩阵R的分布很难求解。为了简化计算，需要对R进行近似处理。当视数不大时，窗口内像素的纹理有很强的相关性。这时，可以认为在窗口内，纹理是固定值，即R≈τI_L。这样，矩阵W变成

W＝τ·ZZ^H(10)

上式认为未标准化的相干矩阵是由纹理τ和未标准化的噪声协方差矩阵相乘得到的。在这种情况下，式(8)变为

$p_W(W;L,\mathrm{\τ},\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{\∫}_0^{\∞}\frac{1}{{\mathrm{\τ}}^{\mathrm{pL}}}\exp \{-\mathrm{tr}\left(\frac{{\mathrm{\Σ}}^{-1}W}{\mathrm{\τ}}\right)\}\·p_{\mathrm{\τ}}\left(\mathrm{\τ}\right)\mathrm{d\τ}---\left(11\right)$

以下基于上述基于极不均匀纹理的乘性模型的概率密度函数对极化SAR乘性模型的参数进行估计，包括：

在式(11)中，需要估计的参数由三类：纹理τ、斑噪协方差矩阵Σ和p_τ(τ)中的纹理分布参数。其中斑噪协方差矩阵的估计是参数估计的关键，首先把W的PDF重写为

$p_W(W;L,\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{h}_M\left[\mathrm{tr}\left({\mathrm{\Σ}}^{-1}W\right)\right]---\left(12\right)$

其中

$M=\mathrm{pL},h_M\left(x\right)={\∫}_0^{\∞}\frac{1}{{\mathrm{\τ}}^M}\exp (-\frac{x}{\mathrm{\τ}})\·p_{\mathrm{\τ}}\left(\mathrm{\τ}\right)\mathrm{d\τ}---\left(13\right)$

对于窗口内的N个像素极化相干矩阵，假设他们之间是独立且同分布(IndependentandIdenticallyDistributed,IID)的，这样N个相干矩阵的似然函数可以表示为

$L_W(W_1,...,W_L;\mathrm{\Σ},{\mathrm{\τ}}_1,...,{\mathrm{\τ}}_L)={\left(\frac{1}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}\right)}^N\underset{i=1}{\overset{N}{\mathrm{\Π}}}{\left|W_i\right|}^{L-p}{h}_M\left[\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{W}_i\right)\right]---\left(14\right)$

求解式中Σ的似然估计值，就是要找到式中似然函数L_W(·)最大时Σ的值。这个过程等价于令lnL_W(·)对Σ的一阶导数为0，即

$-\mathrm{LN}\frac{\∂\ln \left|\mathrm{\Σ}\right|}{\∂\mathrm{\Σ}}+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{g_M\left[\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{W}_i\right)\right]}{h_M\left[\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{W}_i\right)\right]}\frac{\∂\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{W}_i\right)}{\∂\mathrm{\Σ}}=0---\left(15\right)$

此处

$g_M\left(x\right)=\∂h_M\left(x\right)/\∂x---\left(16\right)$

这样，可以得到Σ的似然估计值

$\widehat{\mathrm{\Σ}}=\frac{1}{\mathrm{LN}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{h_{M+1}\left[\mathrm{tr}\left({\widehat{\mathrm{\Σ}}}^{-1}{W}_i\right)\right]}{h_M\left[\mathrm{tr}\left({\widehat{\mathrm{\Σ}}}^{-1}{W}_i\right)\right]}\·W_i---\left(17\right)$

上式为Σ的实际似然估计值。但是，该式仍然包含有h_M(x)，而h_M(x)中依赖于p_τ(τ)。也就是说，式(17)中其实有2个未知数，Σ的估计需要事先知道纹理，这显然是不可能的。

为了解决上面这个问题，需要采用逐个最大化的估计方法，分两步来求解Σ的似然估计值。首先，假设图像纹理是已知的，即认为式(14)中仅有Σ一个参数，这时求解Σ变成了解算下面的方程

$\∂\ln L_W(W_1,...W_N|{\mathrm{\τ}}_1,...,{\mathrm{\τ}}_N;\mathrm{\Σ}\left|\right)/\∂\mathrm{\Σ}=0---\left(18\right)$

在这个假设下，Σ的估计值为

$\widehat{\mathrm{\Σ}}=\frac{1}{\mathrm{LN}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{W_i}{{\widehat{\mathrm{\τ}}}_i}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{T_i}{{\widehat{\mathrm{\τ}}}_i}---\left(19\right)$

然后，再假设斑噪协方差矩阵Σ式已知的，即认为式(14)中仅有τ一个参数，这时求解τ等价于解算下面的方程

$\∂\ln L_W(W_1,...W_N|\mathrm{\Σ};{\mathrm{\τ}}_1,...,{\mathrm{\τ}}_N)/\∂\mathrm{\τ}=0---\left(20\right)$

在这个假设下，τ的估计值为

${\widehat{\mathrm{\τ}}}_i=\mathrm{tr}\left({\widehat{\mathrm{\Σ}}}^{-1}{T}_i\right)/p---\left(21\right)$

这样，联立式(19)和式(21)，就得到了Σ的似然估计值

$\widehat{\mathrm{\Σ}}=\frac{p}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{T_i}{\mathrm{tr}\left({\widehat{\mathrm{\Σ}}}^{-1}{T}_i\right)}---\left(22\right)$

可以通过下面的迭代处理解来求解

${\widehat{\mathrm{\Σ}}}_{t+1}=f\left({\widehat{\mathrm{\Σ}}}_t\right)---\left(23\right)$

式中下标t表示迭代次数。为了让极化SAR的功率信息都表现在纹理τ中，限制以下条件

tr(Σ)＝p(24)

这时，式(23)变为

${\widehat{\mathrm{\Σ}}}_{t+1}=\frac{p}{\mathrm{tr}\left({\widehat{\mathrm{\Σ}}}_t\right)}f\left({\widehat{\mathrm{\Σ}}}_t\right)---\left(25\right)$

此处还需要确立一个迭代的初始值来启动迭代，可以选择单位矩阵或Toeplitz阵来作为迭代的初始矩阵，Toeplitz阵M的表达式为

M_ij＝ρ^|i-j|(26)

此处ρ∈(0,1)。

另外，还需要定义式(25)的迭代终止条件。当和近似相等时，迭代终止，具体条件为

$N\left(t\right)=\frac{{\left|\right|{\widehat{\mathrm{\Σ}}}_{t+1}-{\widehat{\mathrm{\Σ}}}_t\left|\right|}_F}{{\left|\right|{\widehat{\mathrm{\Σ}}}_t\left|\right|}_F}---\left(27\right)$

一般情况下，当N(t)接近10^-16时，N(t)的曲线收敛。这时，可以得到Σ的似然估计值。然后通过式(21)，可以得到τ的估计值。这样就可以得到极化SAR图像的纹理图像和斑噪协方差图像。根据不同的纹理分布假设，可进一步解算p_τ(τ)的分布参数。

以下定义基于极化SAR乘性模型的相似性度量：

在极化SAR乘性模型的基础上，还需要定义一个变化测度(其为相似性度量一种)来定义不同时相极化SAR图像的变化程度。由于纹理和斑噪协方差矩阵两个参数都可能发生变化，可以使用相似性度量的方法来定义变化测度。对于两个不同时相全极化SAR图像中对应的区域X、Y来说，它们的相似性度量为

S＝MLL(X)+MLL(Y)-MLL(X∪Y)(28)

式中，MLL(X)表示区域X的最大对数似然函数，当区域X和Y均无变化时，S＝0。对于任一区域Q，他的MLL函数的表达式为

$\mathrm{MLL}\left(Q\right)=\underset{i\∈Q}{\mathrm{\Σ}}\ln \left(p_T\left(T_i\right|{\mathrm{\θ}}_Q)\right)---\left(29\right)$

θ_Q表示极化相干矩阵PDF的参数集合。

在大多数的极不均匀纹理假设情况下，假设图像纹理服从Gamma分布(使用Gamma分布可以更为精确)，对应的PDF为

$p\left(\mathrm{\τ}\right)=\frac{1}{\mathrm{\Γ}\left(v\right)}\frac{v}{\mathrm{\μ}}{\left(\frac{\mathrm{v\τ}}{\mathrm{\μ}}\right)}^{v-1}{e}^{-\frac{\mathrm{v\τ}}{\mathrm{\μ}}}---\left(30\right)$

此处v为形状参数，μ为分布的均值。这样，未标准化的极化相干矩阵W的概率分布可以表示为

$p_W(W;v,\mathrm{\μ},L,\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L\mathrm{\Γ}\left(v\right)}{\left(\frac{v}{\mathrm{\μ}}\right)}^v{\∫}_0^{\∞}\exp \{-\frac{\mathrm{tr}\left({\mathrm{\Σ}}^{-1}W\right)}{\mathrm{\τ}}-\frac{\mathrm{v\τ}}{\mathrm{\μ}}\}\·{\mathrm{\τ}}^{v-\mathrm{pL}-1}\mathrm{d\τ}---\left(31\right)$

已知第二类修正Bessel函数为

$K_n\left(x\right)=\frac{x^n}{2^{n+1}}{\∫}_0^{+\∞}{t}^{-n-1}{e}^{-t-z^2/4t}\mathrm{dt}---\left(32\right)$

那么式(31)将变为

$p_W(W;v,\mathrm{\μ},L,\mathrm{\Σ})=\frac{2{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L\mathrm{\Γ}\left(v\right)}{\left(\frac{v}{\mathrm{\μ}}\right)}^{v+\mathrm{pL}/2}{\left(\mathrm{tr}\left({\mathrm{\Σ}}^{-1}W\right)\right)}^{v-\mathrm{pL}/2}{K}_{v-\mathrm{pL}}\left(2\sqrt{\mathrm{tr}\left({\mathrm{\Σ}}^{-1}W\right)v/\mathrm{\μ}}\right)---\left(33\right)$

这样可以得到极化相干矩阵T的概率分布

$p_T(T;v,\mathrm{\μ},L,\mathrm{\Σ})=\frac{2{\left|T\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L\mathrm{\Γ}\left(v\right)}{\left(\frac{\mathrm{Lv}}{\mathrm{\μ}}\right)}^{(v+\mathrm{pL})/2}{\left(\mathrm{tr}\left({\mathrm{\Σ}}^{-1}T\right)\right)}^{(v-\mathrm{pL})/2}{K}_{v-\mathrm{pL}}\left(2\sqrt{\mathrm{tr}\left({\mathrm{\Σ}}^{-1}T\right)\mathrm{Lv}/\mathrm{\μ}}\right)---\left(34\right)$

因为，式中包含有第二类修正Bessel函数，因此该分布称为新的矩阵类型K分布。注意到当纹理的均值为1时，即μ＝1，该分布退化为传统的矩阵类型的K分布。

根据中心极限定理，随着形状参数v变大，Gamma分布将变为均值μ，方差为μ²/v的为高斯分布。当v趋于正无穷时，p(τ)变为Diracdelta函数p(τ)＝δ(τ-μ)，也就是说这时为纹理分布的方差变为0，即纹理分布不再随机，变成均匀分布。这时，定义Σ_μ＝μΣ，式(34)中的新的矩阵类型的K分布将变为复Wishart分布。

将式(34)带入到式(29)，并消去无用项，得到的最大对数似然函数为

$\begin{array}{c}{\mathrm{MLL}}_P\left(Q\right)=N\frac{v+\mathrm{pL}}{2}(\ln \mathrm{Lv}-\ln \mathrm{\μ})+\frac{v-\mathrm{pL}}{2}\underset{i\∈Q}{\mathrm{\Σ}}\ln \left[\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{T}_i\right)\right]\\ +\underset{i\∈Q}{\mathrm{\Σ}}\ln \left\{K_{v-\mathrm{pL}}\left(2\sqrt{\mathrm{tr}\left({\mathrm{\Σ}}^{-1}{T}_i\right)\mathrm{Lv}/\mathrm{\μ}}\right)\right\}-\mathrm{NL}\ln \left|\mathrm{\Σ}\right|-N\ln \mathrm{\Γ}\left(v\right)\end{array}---\left(35\right)$

该式中，Σ由式(25)迭代求解，Gamma分布的参数v、μ可根据纹理的估计值使用最大似然或者矩估计解算。

将式(35)带入到(28)中，即可得到一个新的相似性度量

S_P＝MLL_P(X)+MLL_P(Y)-MLL_P(X∪Y)(36)

当v→∞时，该测度将退化为Wishart距离

S_W＝NL(2ln|T_X∪Y|-ln|T_X|-ln|T_Y|)(37)

对比传统的Wishart距离和SIRV距离，可以发现，在根据本发明的实施例中提出的变化测度S_P中，同时包含了纹理信息和协方差矩阵的变化信息。

以下将根据本发明的技术方案应用于实际中进行说明，且便于说明根据本发明实施例的相似性度量的效率，接下来将其与传统的Wishart距离进行对比。使用2009年4月9日与2010年6月15日3视某城区的两幅RadarSat-2全极化影像作为实验数据，图像分辨率为10×8m²，大小为895×540，雷达入射角从38.37°到39.85°。两幅影像的极化灰度合成图像见图3a和图3b，其中，图3a为2009年4月9日Pauli合成图，图3b为2010年6月15日Pauli合成图。该区域中有城区、草地、森林、山地、湖波和裸地等地物，因此非常适合作为极不均匀纹理分布的实验区。

为方便验证变化检测结果，通过实地调研，并结合与两个图像成像时间接近的GoogleEearth光学图像(2009年3月15日和2010年6月19日)，在2010年的极化彩色图像，即图3b上标注出了两幅图像的主要变化区域。主要的变化类型有3类：T1：裸地到城区的变化；T2：城区到裸地的变化；T3：草地到水域的变化。

首先，根据式(25)的迭代算法，使用9×9的滑动窗口来估计两幅极化图像的纹理和斑噪协方差矩阵。一般情况下，估计每个像素斑噪协方差矩阵需要迭代的次数为10～30次，对于2009年的图像，这里平均的迭代次数为15.3。图4a示出了根据本发明实施例的从主图像的像素斑噪协方差矩阵图像内任意选取的像素1进行迭代算法的收敛情况，图4b为根据本发明实施例的从辅图像的像素斑噪协方差矩阵图像内任意选取的像素2进行迭代算法的收敛情况，其中，RelativeError为相对误差，而NumberofIterations为迭代次数，而RelativeErrorconvergence为相对误差的收敛性，可以看到式(27)的N(t)收敛的很快，到达10^-16后，N(t)趋于平稳。

图5a和图5b给出了通过迭代算法得到的两副纹理图像，其中，图5a为2009年得到的影像(即主图像)通过迭代算法得到的纹理图像，图5b为2010年得到的影像(即辅图像)通过迭代算法得到的纹理图像。由于纹理中只包含极化SAR数据的功率信息，而除去了斑噪干扰，因此求取纹理的过程相对于极化SAR图像去噪的过程。比较纹理图与全极化相干矩阵的灰度合成图(图3a和图3b)发现，纹理图明显噪声明显减少，图像更为清晰，对比度也更高。对比两个时相的纹理图，可以很容易找到两副图像之间的变化区域。

当对两幅极化SAR图像拟合后，还需要使用根据本发明实施例的相似性度量S_P来求解这两幅图像的差异图。为了方便对比，也计算了Wishart距离所得到的差异图，两幅差异图如图6a和图6b所示，其中，图6a为基于本发明实施例的相似性度量得到的主、辅图像的差异图，图6b为基于现有技术中的Wishart距离得到的主、辅图像的差异图。两个方法都检测出了图3b中标记的主要变化区域，但对于图中4个画圈的区域来说，Wishart距离得到的检测值偏小。从图5a和图5b可知，画圈的区域纹理变化并不大，但从图3a和图3b可以发现这些区域的极化散射特性发生了变化。在极化SAR乘性模型中，斑噪协方差矩阵表示图像的极化信息。因此，Wishart距离对斑噪协方差矩阵的变化不敏感。

由于基于本发明实施例的相似性度量可以很好检测极化散射特性的变化，这使得S_P得到的差异图显得比Wishart距离的差异图有更强的“斑噪”，但其实基于S_P的差异图的等效视数是基于Wishart距离的差异的2倍，这将有利于差异图的解译和后面阈值分割。

为了进一步比较这两种变化测度，绘制了上述两个差异图的ROC(即相对误差的收敛性)曲线。已知在ROC性能评价中，ROC曲线下方所覆盖的面积越大，则表示该算子性能越好，用其对这两个变化测度分别进行性能评估，结果如图7所示，即，图7为两个变化测度(即差异图)的ROC曲线。Proposedsimilaritymeasure表示根据本发明实施例所提出的相似性度量，其对应的曲线既为根据本发明的相似性度量的ROC曲线，Wishartsimilaritymeasure表示Wishart距离对应相似性度量，即现有技术中基于Wishart距离的ROC曲线，Falseararmrate为虚警率，Detectionrate为检测率。从图中可以看出，根据本发明实施例的变化测度要优于Wishart距离。

差异图提取后，为在原图上准确定位变化区域，需将差异图提取为变化与未变化的二值化掩膜，将掩膜与原图融合得到检测结果图。图8a为对根据本发明实施例的相似性度量得出的差异图(即图6a所示的图)的二值化掩膜进行KI阈值分割后的示意图，图8b为对图8a进行连通区搜索后的示意图。

在得到二值图像后，对其进行形态学闭运算。对闭运算中膨胀运算和腐蚀运算的结果求差，就可以得到变化区域的边界，如图9所示为根据本发明的实施例对图8b所示图进行形态学运算之后的示意图，有利于变化区域面积和周长的统计。

由于阈值分割和形态学处理相当于对差异图进行二值化，即将差异图分为“变”和不变的两类。因此，可以用混淆矩阵还对其进行精度评价，评价结果见表1。可以发现，基于根据本发明实施例的相似性度量的方法的精度和Kappa系数都比基于Wishart距离的方法高。但是，在总精度上，两个方法比较接近。这主要是因为在变化检测中，未变化的区域一般都远大于变化区域，这导致混淆矩阵的对角线元素很大，从而增加了分类的总体精度。因此，用分类的总体精度并不能有效反应检测方法的优劣。

表1分类精度评价

在变化检测中，两个更为重要的精度指标是检测率和检测虚警率。检测率即是为表1中变化区域的生产者精度，而虚警率为变化区域的用户精度与1的差。两个方法的检测率和虚警率见表2。相对于传统的基于Wishart距离的检测方法，根据本发明实施例的检测方法有更高的检测率和更低的虚警率。

表2检测精度评价

当得到变化区域后，就可以把变化区域与变化前后的图像融合，制成土地利用图。由于极化SAR可以提供地物的类别的信息，将得到的变化区域与2009年的分类图叠合，如图10a、图10b和图10c，其中，图10a为将图9中的差异图与2009年该地区的分类图(即主图像的分类图)进行融合后的示意图，图10b为图10a所示的示意图在2010年(即辅图像的分类图)的示意图，以及图10c为图10a和图10b所示图中的各部分所分别代表的地理形态的示意图。从图10a和图10b中可以直接看到图9中的差异图所反应的变化区域的实际位置，并可以确定变化前后的地物类型。这种土地利用图可以为土地资源管理、城市变化以及农田耕地变迁等应用服务。

由上可知，根据国内外极化SAR影像定量应用的巨大需求，围绕高分辨率SAR可提升传统应用中的精度、准确度要求和对新型应用扩展、创新的要求，本发明挖掘高分辨率极化SAR的应用潜力，着重解决目前极化SAR图像杂波拟合和变化测度模型中存在的问题，深入开展极化SAR影像高精度时序分析研究。

为了验证提出变化检测方法的有效性，使用某地区的两景3视Radarsat-2全极化数据开展相关试验，比较了根据本发明实施例的检测方法与基于Wishart距离的检测方法，实验表明根据本发明实施例的检测方法有更高的检测率和更低的虚警率。

以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。

Claims (10)

1.一种检测全极化SAR图像的变化区域的方法，其特征在于，所述方法包括：

步骤S2：使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到所述主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据所述主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，所述主、辅SAR图像为包括同一对象区域的高分辨率全极化SAR图像；

步骤S4：根据所述主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过所述基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比所述主、辅SAR图像得到差异信息，基于所述差异信息确定相对于所述主、辅SAR图像的差异图像；

步骤S6：将所述差异图像与所述主SAR图像和/或所述辅SAR图像进行比对，以检测出所述变化区域。

2.根据权利要求1所述的方法，其特征在于，在步骤S2之前，所述方法进一步包括：

步骤S1：对接收的所述主、辅SAR图像进行所述同一对象区域的对应匹配。

3.根据权利要求2所述的方法，其特征在于，步骤S2中的使用所述基于极不均匀纹理的多视极化SAR乘性模型对所述主、辅SAR图像进行拟合包括：使用基于极不均匀纹理的乘性模型的概率密度函数推导所述主、辅SAR图像的纹理图像、斑噪协方差图像。

所述基于极不均匀纹理的乘性模型的概率密度函数为：

$p_W(W;L,\mathrm{\τ},\mathrm{\Σ})=\frac{{\left|W\right|}^{L-p}}{{\mathrm{\Γ}}_p\left(L\right){\left|\mathrm{\Σ}\right|}^L}{\∫}_0^{\∞}\frac{1}{{\mathrm{\τ}}^{\mathrm{pL}}}\exp \{-\mathrm{tr}\left(\frac{{\mathrm{\Σ}}^{-1}W}{\mathrm{\τ}}\right)\}\·p_{\mathrm{\τ}}\left(\mathrm{\τ}\right)\mathrm{d\τ}$

其中，L代表所述多视为L视，W为所述L视下未标准化的全极化相干矩阵，Pw为W的概率密度函数，τ表示图像的纹理，Σ为斑噪协方差矩阵，Γ_p(L)为L视下概率密度函数的复数形式的多元Gamma函数。

4.根据权利要求3所述的方法，其特征在于，步骤S2中的分别根据所述主、辅SAR图像的纹理图像确定其对应的纹理图像的概率密度函数的分布参数包括：

设置所述主、辅SAR图像的纹理图像以符合Gamma函数分布。

5.根据权利要求4所述的方法，其特征在于，所述基于极不均匀纹理的多视极化SAR乘性模型的相似性度量为：

S_P＝MLL_P(X)+MLL_P(Y)-MLL_P(X∪Y)

其中，Sp为所述Pw的变化测度，X为所述主SAR图像内的子区域，MLLp(X)表示区域X内Pw的最大对数似然函数，Y为所述辅SAR图像内的子区域，MLLp(Y)表示区域Y内Pw的最大对数似然函数，区域X与区域Y为对应子区域。

6.根据权利要求1所述的方法，其特征在于，在步骤S6之前，包括：

步骤S5：通过对所述差异图像进行阀值分割来确定其区域和边界。

7.根据权利要求6所述的方法，其特征在于，步骤S5包括：

S51、对所述差异图像进行自适应KI阀值分割得到区别变化区域与不变区域的作为二值化图像的二值化掩膜；

S52、通过连通区搜索算法去除所述二值化图像中的变化面积小于预定值的所述变化区域；

S53、通过形态学的闭运算中的膨胀操作对所述变化区域进行膨胀得到膨胀区域，并对该膨胀区域的周界进行记录；

S54、通过形态学的闭运算中的腐蚀操作，采用与所述膨胀操作相对应的比例对所述膨胀区域进行腐蚀，得到腐蚀后区域；

S55、将所述腐蚀后区域的周界与所述膨胀区域的周界进行作差得到所述差异图像中所述变化区域的边界。

8.根据权利要求6所述的方法，其特征在于，在步骤S6之后进一步包括：

步骤S7：将所述变化区域与所述主、辅SAR图像进行叠合，确定所述主、辅SAR图像内的变化部分。

9.一种检测全极化SAR图像的变化区域的装置，其特征在于，所述全极化SAR图像基于极不均匀纹理，所述检测装置包括：

参数确定模块，其配置为使用基于极不均匀纹理的多视极化SAR乘性模型对先后接收的全极化主、辅SAR图像进行拟合，以分别得到所述主、辅SAR图像的纹理图像、斑噪协方差图像，并分别根据所述主、辅SAR图像的纹理图像确定其对应的概率密度函数的分布参数，所述主、辅SAR图像为包括同一对象区域的高分辨率全极化SAR图像；

差异信息确定模块，其配置为根据所述主、辅SAR图像各自的纹理图像、斑噪协方差图像和概率密度函数的分布参数，并通过所述基于极不均匀纹理的多视极化SAR乘性模型的相似性度量对比所述主、辅SAR图像得到差异信息，基于所述差异信息确定相对于所述主、辅SAR图像的差异图像；

比对模块，其配置为将所述差异图像与所述主SAR图像和/或所述辅SAR图像进行比对，以检测出所述变化区域。

10.根据权利要求9所述的装置，其特征在于，进一步包括：

区域匹配模块，其配置为对接收的所述主、辅SAR图像进行所述同一对象区域的对应匹配。