CN104166128A - 基于广义似然比的多航过sar相干变化检测方法 - Google Patents

Abstract

本发明公开了一种基于广义似然比的多航过SAR相干变化检测方法，包括以下步骤：S1：选取多航过SAR图像对记为{f₁,f₂,…,f_K}；S2：选取多航过SAR图像像素对；S3：对协方差矩阵进行最大似然估计；S4：进行似然比假设检验；S5：依次选取多航过SAR图像像素对，重复步骤S4，得到变化检测结果。本发明假设SAR成像区域在多幅图像采集期间发生变化和未发生变化时其对应的复像素对分别服从不同的圆对称复高斯分布，对假设中圆对称复高斯分布的协方差矩阵进行估计，然后确定检测统计量，并将该检测统计量和门限比较，检验上述两个假设的成立与否，也就是检测成像区域有无变化，可以实现微弱变化检测和变化过程的观测。

Description

基于广义似然比的多航过SAR相干变化检测方法

技术领域

本发明属于合成孔径雷达(SAR)技术领域，特别涉及一种基于广义似然比的多航过SAR相干变化检测方法。 

背景技术

变化检测技术可以广泛应用于监测森林植被、土壤水分、积雪量的变化；监测农作物的长势变化、土地覆盖的变化；监测各种灾害发生前后的变化，如地震区域的定位和灾害评估；监测海冰的运动、山岳冰川的移动以及滑坡运动；军事目标区域的动态监测、战场打击评估等。但是在云、雾等恶劣天气条件下难以通过光学图像获取变化信息。由于SAR是一种全天时、全天候的现代高分辨率微波遥感成像雷达，通过对获取的SAR图像进行变化检测，能够在第一时间为决策提供有力的支持。然而当目标变化微弱时，传统SAR图像变化检测方法检测不到变化。另外，受森林、海浪等微动杂波的影响，传统SAR图像变化检测方法使目标的变化淹没在周围背景的变化中，区分不出目标变化与背景变化的差异，使得变化检测性能急剧下降。 

多航过SAR是目前微波遥感技术研究与应用的一个重要领域，它深入挖掘信息有效地提高了雷达对地物信息的获取能力，对多航过雷达遥感图像进行信息处理的研究具有重要的理论价值和广阔的应用前景。 

近年来，国外很多学者开始探索利用多航过极化微波遥感获取同一地区的多幅图像，经过信号处理检测目标微弱变化。在文献“Leslie M Novak，Change detection for multi-polarization,multi-pass SAR，in Defense and Security.International Society for Optics and Photonics,2005,pp.234–246”中提出了一种基于广义似然比的变化检测方法，该方法利用多航过图像的幅度信息提高了变化检测性能，但是此非相干变化检测方法通过场景的后向散射功率的变化来进行变化检测，没有利用相位信息，对SAR图像变化的灵敏度低。然而雷达系统是相干系统，回波的幅度和相位对目标的变化都非常敏感。单纯利用幅度信息进行变化检测的方法丢失了大量的有用的相位信息已经不能满足精确变化检测需求。在文献“Mark Preiss and Nicholas JS Stacy,Coherent Change Detection：Theoretical Description and Experimental Results,Tech.Rep,DTIC Document,2006”中提出了一种相干变化检测方法，该方法同时利用图像的幅度和相位信息，能够检测图像微弱变化。但是该方法只利用两幅图像的信息，没有利用多幅图像的时相信息。 

发明内容

本发明的目的在于克服两航过非相干变化检测方法中只利用幅度信息以及仅利用两航 过SAR图像信息导致变化检测概率低或检测不到变化的问题，提供一种假设SAR成像区域在多幅图像采集期间发生变化和未发生变化时其对应的复像素对分别服从不同的圆对称复高斯分布，再确定检测统计量并与门限比较以检验上述两个假设的成立与否的基于广义似然比的多航过SAR相干变化检测方法。 

本发明的目的是通过以下技术方案来实现的：基于广义似然比的多航过SAR相干变化检测方法，包括以下步骤： 

S1：选取多航过SAR图像对：选取不同时间对同一地区多次观测获得的K幅SAR图像并进行配准后记为{f₁,f₂,…,f_K}，K≥3； 

S2：选取多航过SAR图像像素对，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的值，并记为向量假设向量X_mn～CN(0,Γ)，其中Γ＝E[X_mnX_mn ^H]为协方差矩阵，设H₀表示目标区域未发生变化，得到协方差矩阵Γ₀，设H₁表示目标区域发生变化，得到协方差矩阵Γ₁； 

S3：对协方差矩阵进行最大似然估计：分别对协方差矩阵Γ₀和Γ₁进行最大似然估计； 

S4：进行似然比假设检验； 

S5：依次选取多航过SAR图像像素对，重复步骤S4，得到变化检测结果。 

进一步地，所述的步骤S2的具体实现方法为：由于图像存储为矩阵形式，因此f₁,f₂,…,f_K均为复图像，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的值分别记为 令复向量m＝1,2,…,M，n＝1,2,…,N，其中[]^T表示转置运算，M、N为图像对应的复矩阵的大小； 

设向量X_mn服从K维圆对称复高斯分布，即X_mn～CN(0,Γ)，其中协方差矩阵Γ＝E[X_mnX_mn ^H]，因此X_mn的概率密度函数表示为： 

$p\left(X_{\mathrm{mn}}\right)=\frac{1}{{\mathrm{\π}}^K\left|\mathrm{\Γ}\right|}\exp (-{X_{\mathrm{mn}}}^H\mathrm{\Γ}{X_{\mathrm{mn}}}^H)$

其中，E[X_mnX_mn ^H]表示求X_mnX_mn ^H的均值，X_mn ^H表示求X_mn的复共轭转置，|Γ|表示Γ的行列式，exp表示指数运算； 

设H₀表示目标区域未发生变化，此时对应的向量X_mn～CN(0,Γ₀)，其中， 

其中，j表示复数单元，ρ_ab≈1，Φ_ab≈0°，a＝1，2，…，K，b＝1，2，…，K，a＜b； 

设H₁表示目标区域发生变化，此时对应的向量X_mn～CN(0,Γ₁)，其中， 

ρ′_ab≈0，Φ′_ab≠0°，由于Γ₀、Γ₁未知，因此为复合假设。 

进一步地，所述的步骤S3对协方差矩阵Γ₀和Γ₁进行最大似然估计的具体方法为： 

S31：对Γ₀进行最大似然估计，根据最大似然估计理论，对于假设H₀，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显未变化区域得到参数为Γ₀的复高斯随机向量Y的L个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₀)，Γ₀的似然函数lik(Γ₀)表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_0\right)=\underset{i=1}{\overset{L}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_0)=\frac{1}{{\mathrm{\π}}^{\mathrm{KL}}{\left|{\mathrm{\Γ}}_0\right|}^L}{\exp }^{-\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_0}^{-1}{Y}_i}$

令求得协方差矩阵Γ₀的最大似然估计： 

$\widehat{{\mathrm{\Γ}}_0}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y}_i{Y_i}^H;$

S32：对Γ₁进行最大似然估计，对于假设H₁，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显变化区域得到的参数为Γ₁的复高斯随机向量Y的S个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₁)，Γ₁的似然函数lik(Γ₁)可以表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_1\right)=\underset{i=1}{\overset{S}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_1)=\frac{1}{{\mathrm{\π}}^{\mathrm{KL}}{\left|{\mathrm{\Γ}}_1\right|}^S}{\exp }^{-\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_1}^{-1}{Y}_i}$

令求出协方差矩阵Γ₁的最大似然估计： 

$\widehat{{\mathrm{\Γ}}_1}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{Y}_i{Y_i}^H.$

进一步地，所述的步骤S4包括以下子步骤： 

S41：确定检验统计量，根据步骤S2的假设和步骤S3的参数估计得到： 

$p\left(X_{\mathrm{mn}}\right|H_0)=\frac{1}{{\mathrm{\π}}^K\left|\widehat{{\mathrm{\Γ}}_0}\right|}\exp (-{X_{\mathrm{mn}}}^H{\widehat{{\mathrm{\Γ}}_0}}^{-1}{X}_{\mathrm{mn}})$

$p\left(X_{\mathrm{mn}}\right|H_1)=\frac{1}{{\mathrm{\π}}^K\left|\widehat{{\mathrm{\Γ}}_1}\right|}\exp (-{X_{\mathrm{mn}}}^H{\widehat{{\mathrm{\Γ}}_1}}^{-1}{X}_{\mathrm{mn}})$

取K幅多航过SAR图像中第m行第n列像素点的Q个相互独立的像素对根据似然比假设检验理论，令似然比： 

$z=\frac{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^Q;H_0)}{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^Q;H_1)}=\underset{i=1}{\overset{Q}{\mathrm{\Π}}}\frac{p(X_{\mathrm{mn}}^i;H_0)}{p(X_{\mathrm{mn}}^i;H_1)}$

将代入上式化简得， 

$z={\left(\frac{\left|\widehat{{\mathrm{\Γ}}_1}\right|}{\left|\widehat{{\mathrm{\Γ}}_0}\right|}\right)}^Q\exp (-\mathrm{Tr}\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\left\}\right)$

取对数并忽略常数项，得检测统计量： 

$Z=\mathrm{Tr}\left\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\right\}=\mathrm{Tr}\left\{\widehat{{\mathrm{\Γ}}_{d}G}\right\}$

其中表示矩阵G的迹，也就是矩阵对角元素的和， 

${\widehat{\mathrm{\Γ}}}_d={{\widehat{\mathrm{\Γ}}}_0}^{-1}-{{\widehat{\mathrm{\Γ}}}_1}^{-1},G=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H;$

S42：进行假设检验，定义一个和用来检测的图像一样大小的矩阵，记为R，选择门限T，判断门限T与检测统计量Z的大小：当Z>T时，判定假设H₁成立，变化检测结果就是该像素对应的区域发生变化，令R对应的像素值为255；否则判定假设H₀成立，变化检测结果就是该像素对应的区域未发生变化，令R对应的像素值为0。 

进一步地，所述的步骤S5中得到变化检测结果的方法为：依次选取多航过SAR图像像素对，重复步骤S4，直到确定矩阵R中所有的像素值，R即为变化检测结果。 

本发明的有益效果是： 

本发明利用广义似然比假设检验理论，假设SAR成像区域在多幅图像采集期间发生变 化和未发生变化时其对应的复像素对分别服从不同的圆对称复高斯分布，对假设中圆对称复高斯分布的协方差矩阵进行估计，然后确定检测统计量，并将该检测统计量和门限比较，检验上述两个假设的成立与否，也就是检测成像区域有无变化。 

本发明的检测方法充分利用多航过SAR复图像丰富的幅度和相位信息，与现有的两航过相干变化检测方法和多航过非相干变化检测方法相比，可以实现微弱变化检测和变化过程的观测，在地球遥感和地质灾害监测等领域十分适用。 

附图说明

图1为本发明的检测方法的流程图； 

图2为本发明实施例采用的某一地区的SAR图像f₁； 

图3为本发明实施例采用的某一地区的SAR图像f₂； 

图4为本发明实施例采用的某一地区的SAR图像f₃； 

图5为本发明实施例变化检测结果图像； 

图6为本发明实施例不同航过相干变化检测的ROC曲线。 

具体实施方式

下面结合附图进一步说明本发明的技术方案，但本发明所保护的内容不局限于以下所述。 

为了方便描述本发明的内容，首先作以下解释： 

1、复高斯分布 

a)复高斯分布的定义 

假设X和Y是k维实空间的随机向量，向量vect[X Y]是2k维正态随机向量。那么以X为实部、Y为虚部的复随机向量Z＝X+Yj具有复高斯分布，记为Ζ～CN(μ,Γ,C)。 

μ＝E[Z],Γ＝E[(Z-μ)(Z-μ)^H],C＝E[(Z-μ)(Z-μ)']， 

其中：j是复数单元，即—1开根，E[Z]表示求Z的均值，均值μ可以是任意k维的复矢量；Z^H表示求Z的复共轭转置，协方差矩阵Γ必须是厄尔米特的和非负定的；Z'表示求Z矩阵转置，相关矩阵C是对称的。 

b)圆对称复高斯分布 

圆对称复高斯分布对应的参数为μ＝0，C＝0。若k维随机向量Z＝X+iY服从圆对称复高斯分布，通常写做Ζ～CN(0,Γ)，其概率密度函数为: 

$f\left(Z\right)=\frac{1}{{\mathrm{\π}}^k\left|\mathrm{\Γ}\right|}{\exp }^{-Z^H{\mathrm{\Γ}}^{-1}Z}$

其中|Γ|表示Γ的行列式，exp表示指数运算，Γ^-1表示Γ的逆。 

2、复合假设 

虚无与对立假设中不只包含一个母数值。在一般的假设检验中还存在一个未知的参数。 

3、最大似然估计 

假设得到参数α的随机变量y的m个观测。每个观测y_i有pdfp(y_i|α)，并且由于m个观测是独立的，联合称之为似然函数，可以表示为 

$p\left(\stackrel{\→}{y}\right|\mathrm{\α})=\underset{i=1}{\overset{m}{\mathrm{\Π}}}p\left(y_i\right|\mathrm{\α})$

其中Π表示连乘运算。那么通过求似然函数的最大值可以确定参数α的最大似然估计计算时，这个最大值通过令似然函数的导数为零求得，也就是 

$\frac{\∂}{\∂\mathrm{\α}}p\left(\stackrel{\→}{y}\right|\mathrm{\α})=0$

由于对数函数是单调函数，所以常用一个更简单的过程来求似然函数对数的最大值(称为对数-似然)，也就是 

$\frac{\∂}{\∂\mathrm{\α}}\ln p\left(\stackrel{\→}{y}\right|\mathrm{\α})=0$

最大似然估计的前提是似然函数导数存在。 

4、虚警概率、检测概率 

本发明中虚警概率是指场景未发生变化但被检测为变化的概率。检测概率指场景发生变化同时该变化也被检测出的概率。 

5、ROC曲线 

通过改变门限得到多组检测概率和虚警概率，以虚警概率为横坐标，检测概率为纵坐标，根据这些点连成的曲线称为ROC曲线。 

本发明提供一种基于广义似然比的多航过SAR相干变化检测方法，如图1所示，包括以下步骤： 

S1：选取多航过SAR图像对：选取不同时间对同一地区多次观测获得的K幅SAR图像并进行配准后记为{f₁,f₂,…,f_K}，K≥3； 

S2：选取多航过SAR图像像素对，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的 值，并记为向量假设向量X_mn～CN(0,Γ)，其中Γ＝E[X_mnX_mn ^H]为协方差矩阵，设H₀表示目标区域未发生变化，得到协方差矩阵Γ₀，设H₁表示目标区域发生变化，得到协方差矩阵Γ₁； 

S3：对协方差矩阵进行最大似然估计：分别对协方差矩阵Γ₀和Γ₁进行最大似然估计； 

S4：进行似然比假设检验； 

S5：依次选取多航过SAR图像像素对，重复步骤S4，得到变化检测结果。 

进一步地，所述的步骤S2的具体实现方法为：由于图像存储为矩阵形式，因此f₁,f₂,…,f_K均为复图像，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的值分别记为 令复向量m＝1,2,…,M，n＝1,2,…,N，其中[]^T表示转置运算，M、N为图像对应的复矩阵的大小； 

设向量X_mn服从K维圆对称复高斯分布，即X_mn～CN(0,Γ)，其中协方差矩阵Γ＝E[X_mnX_mn ^H]，因此X_mn的概率密度函数表示为： 

$p\left(X_{\mathrm{mn}}\right)=\frac{1}{{\mathrm{\π}}^K\left|\mathrm{\Γ}\right|}\exp (-{X_{\mathrm{mn}}}^H\mathrm{\Γ}{X_{\mathrm{mn}}}^H)$

其中，E[X_mnX_mn ^H]表示求X_mnX_mn ^H的均值，X_mn ^H表示求X_mn的复共轭转置，|Γ|表示Γ的行列式，exp表示指数运算； 

设H₀表示目标区域未发生变化，此时对应的向量X_mn～CN(0,Γ₀)，其中， 

其中，j表示复数单元，ρ_ab≈1，Φ_ab≈0°，a＝1，2，…，K，b＝1，2，…，K，a＜b； 

设H₁表示目标区域发生变化，此时对应的向量X_mn～CN(0,Γ₁)，其中， 

ρ′_ab≈0，Φ′_ab≠0°，由于Γ₀、Γ₁未知，因此为复合假设。 

进一步地，所述的步骤S3对协方差矩阵Γ₀和Γ₁进行最大似然估计的具体方法为： 

S31：对Γ₀进行最大似然估计，根据最大似然估计理论，对于假设H₀，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显未变化区域得到参数为Γ₀的复高斯随机向量Y的L个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₀)，Γ₀的似然函数lik(Γ₀)表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_0\right)=\underset{i=1}{\overset{L}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_0)=\frac{1}{{\mathrm{\π}}^{\mathrm{KL}}{\left|{\mathrm{\Γ}}_0\right|}^L}{\exp }^{-\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_0}^{-1}{Y}_i}$

令求得协方差矩阵Γ₀的最大似然估计： 

$\widehat{{\mathrm{\Γ}}_0}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y}_i{Y_i}^H;$

S32：对Γ₁进行最大似然估计，对于假设H₁，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显变化区域得到的参数为Γ₁的复高斯随机向量Y的S个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₁)，Γ₁的似然函数lik(Γ₁)可以表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_1\right)=\underset{i=1}{\overset{S}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_1)=\frac{1}{{\mathrm{\π}}^{\mathrm{KS}}{\left|{\mathrm{\Γ}}_1\right|}^S}{\exp }^{-\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_1}^{-1}{Y}_i}$

令求出协方差矩阵Γ₁的最大似然估计： 

$\widehat{{\mathrm{\Γ}}_1}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{Y}_i{Y_i}^H.$

进一步地，所述的步骤S4包括以下子步骤： 

S41：确定检验统计量，根据步骤S2的假设和步骤S3的参数估计得到： 

$p\left(X_{\mathrm{mn}}\right|H_0)=\frac{1}{{\mathrm{\π}}^K\left|\widehat{{\mathrm{\Γ}}_0}\right|}\exp (-{X_{\mathrm{mn}}}^H{\widehat{{\mathrm{\Γ}}_0}}^{-1}{X}_{\mathrm{mn}})$

$p\left(X_{\mathrm{mn}}\right|H_1)=\frac{1}{{\mathrm{\π}}^K\left|\widehat{{\mathrm{\Γ}}_1}\right|}\exp (-{X_{\mathrm{mn}}}^H{\widehat{{\mathrm{\Γ}}_1}}^{-1}{X}_{\mathrm{mn}})$

取K幅多航过SAR图像中第m行第n列像素点的Q个相互独立的像素对根据似然比假设检验理论，令似然比： 

$z=\frac{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^Q;H_0)}{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^Q;H_1)}=\underset{i=1}{\overset{Q}{\mathrm{\Π}}}\frac{p(X_{\mathrm{mn}}^i;H_0)}{p(X_{\mathrm{mn}}^i;H_1)}$

将代入上式化简得， 

$z={\left(\frac{\left|\widehat{{\mathrm{\Γ}}_1}\right|}{\left|\widehat{{\mathrm{\Γ}}_0}\right|}\right)}^Q\exp (-\mathrm{Tr}\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\left\}\right)$

取对数并忽略常数项，得检测统计量： 

$Z=\mathrm{Tr}\left\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\right\}=\mathrm{Tr}\left\{\widehat{{\mathrm{\Γ}}_{d}G}\right\}$

其中表示矩阵G的迹，也就是矩阵对角元素的和， 

${\widehat{\mathrm{\Γ}}}_d={{\widehat{\mathrm{\Γ}}}_0}^{-1}-{{\widehat{\mathrm{\Γ}}}_1}^{-1},G=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H;$

S42：进行假设检验，定义一个和用来检测的图像一样大小的矩阵，记为R，选择门限T，判断门限T与检测统计量Z的大小：当Z>T时，判定假设H₁成立，变化检测结果就是该像素对应的区域发生变化，令R对应的像素值为255；否则判定假设H₀成立，变化检测结果就是该像素对应的区域未发生变化，令R对应的像素值为0。 

进一步地，所述的步骤S5中得到变化检测结果的方法为：依次选取多航过SAR图像像素对，重复步骤S4，直到确定矩阵R中所有的像素值，R即为变化检测结果。 

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2013a上验证正确，下面结合具体实施例进一步说明本发明的技术方案：本实施方式中以三航过变化检测为例，首先选取不同时间对同一地区多次观测获得的3幅SAR图像，如图2、图3和图4所示，经过配准后，记为{f₁,f₂,f₃}，对比f₁可知f₂在图中黑色圆圈所表示区域发生较小的变化，对比f₁可知f₃在图中黑色椭圆所表示区域发生较大范围的变化。 

依次选取f₁,f₂,f₃复图像对应的复矩阵第m行第n列的值分别记为令复向量m＝1,2,…,451，n＝1,2,…,451，假设向量X_mn服从3维圆对称复高斯分布，X_mn的概率密度函数可表示为： 

$p\left(X_{\mathrm{mn}}\right)=\frac{1}{{\mathrm{\π}}^3\left|\mathrm{\Γ}\right|}\exp (-{X_{\mathrm{mn}}}^H\mathrm{\Γ}{X_{\mathrm{mn}}}^H)$

假设H₀为目标区域未发生变化，此时对应的向量X_mn～CN(0,Γ₀)，假设H₁为目标区域发生变化，此时对应的向量X_mn～CN(0,Γ₁)。 

对Γ₀进行最大似然估计：根据最大似然估计理论，对于假设H₀，根据{f₁,f₂,f₃}中明显未变化区域我们得到参数为Γ₀的随机向量Y的1681个相互独立的观测，每个观测Y_i的概率密度函数为p(Y_i|Γ₀)，Γ₀的似然函数lik(Γ₀)可以表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_0\right)=\underset{i=1}{\overset{1681}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_0)=\frac{1}{{\mathrm{\π}}^{5043}{\left|{\mathrm{\Γ}}_0\right|}^{1681}}{\exp }^{-\underset{i=1}{\overset{1681}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_0}^{-1}{Y}_i}$

令求得协方差矩阵Γ₀的最大似然估计：

对Γ₁进行最大似然估计：对于假设H₁，根据{f₁,f₂,f₃}中明显变化区域我们得到的参数为Γ₁的随机向量X的121个相互独立的观测，每个观测Y_i的概率密度函数为p(Y_i|Γ₁)，Γ₁的似然函数lik(Γ₁)可以表示为： 

$\mathrm{lik}\left({\mathrm{\Γ}}_1\right)=\underset{i=1}{\overset{121}{\mathrm{\Π}}}p\left(Y_i\right|{\mathrm{\Γ}}_0)=\frac{1}{{\mathrm{\π}}^{363}{\left|{\mathrm{\Γ}}_1\right|}^{121}}{\exp }^{-\underset{i=1}{\overset{121}{\mathrm{\Σ}}}{Y_i}^H{{\mathrm{\Γ}}_1}^{-1}{Y}_i}$

令求出协方差矩阵Γ₁的最大似然估计：

取3幅多航过SAR图像中第m行第n列像素点的9个相互独立的像素对 根据似然比假设检验理论，令似然比： 

$z=\frac{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^9;H_0)}{p(X_{\mathrm{mn}}^1,X_{\mathrm{mn}}^2,...X_{\mathrm{mn}}^9;H_1)}=\underset{i=1}{\overset{9}{\mathrm{\Π}}}\frac{p(X_{\mathrm{mn}}^i;H_0)}{p(X_{\mathrm{mn}}^i;H_1)}$

将代入上式得， $z={\left(\frac{\left|\widehat{{\mathrm{\Γ}}_1}\right|}{\left|\widehat{{\mathrm{\Γ}}_0}\right|}\right)}^9\exp (-\mathrm{Tr}\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\left\}\right)$ 取对数并忽略常数项，得检测统计量： $Z=\mathrm{Tr}\left\{({\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1})\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H\right\}=\mathrm{Tr}\left\{\widehat{{\mathrm{\Γ}}_{d}G}\right\}$ 其中 $\widehat{{\mathrm{\Γ}}_d}={\widehat{{\mathrm{\Γ}}_0}}^{-1}-{\widehat{{\mathrm{\Γ}}_1}}^{-1},G=\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{X}_{\mathrm{mn}}^i{X_{\mathrm{mn}}^i}^H.$

定义一个451×451的矩阵，记为R。选择门限T＝0，当Z>0时，判定假设H₁成立，令R对应的像素值为255，否则判定假设H₀成立，令R对应的像素值为0。 

依次选取像素对，重复上述操作，直到确定矩阵R中所有的像素值，得到如图5所示的变化检测结果图像。 

下面对本发明得到的变化检测方法进行性能分析： 

(1)计算一组检测概率、虚警概率：分别生成25个二维圆对称复高斯随机变量和三维圆对称复高斯随机变量其中i＝1,2,…,25， 

$X_i^{3\′}~\mathrm{CN}(0,{\mathrm{\Γ}}_{30}^{\′}),Y_i^{3\′}~\mathrm{CN}(0,{\mathrm{\Γ}}_{31}^{\′}),{\mathrm{\Γ}}_{20}^{\′}=10\left[\begin{array}{cc}2.2686& 0.45\·2.2686\\ 0.45\·2.2686& 2.2686\end{array}\right],$

${\mathrm{\Γ}}_{21}^{\′}=10\left[\begin{array}{cc}2.2686& 0\\ 0& 0.95070\end{array}\right],$

${\mathrm{\Γ}}_{30}^{\′}=10\left[\begin{array}{ccc}2.2686& 0.45\·2.2686& 0.45\·2.2686\\ 0.45\·2.2686& 2.2686& 0.45\·2.2686\\ 0.45\·2.2686& 0.45\·2.2686& 2.2686\end{array}\right],{\mathrm{\Γ}}_{30}^{\′}=10\left[\begin{array}{ccc}2.2686& 0& 0\\ 0& 1.0847& 0\\ 0& 0& 0.95070\end{array}\right].$

可以看出Γ'₂₀,Γ'₃₀满足步骤2中的H₀假设，Γ'₂₁,Γ'₃₁满足步骤2中的H₁假设。令门限T＝0，由计算检测统计量如果大于门限，则记为虚警。由计算检测统计量如果大于门限，则记为正确检测。同样地，由计算检测统计量如果大于门限，则记为虚警。由计算检测统计量如果大于门限，则记为正确检测。重复10000次上述步骤，统计虚警次数和正确检测的次数并分别除以10000，即得门限为0时两航过和三航过变化检测的虚警概率和检测概率。 

(2)用蒙特卡洛仿真方法求ROC曲线：在[-200,200]以1为间隔均匀改变步骤(1)中门限的值，计算401组检测概率和虚警概率。在以虚警概率为横轴，检测概率为纵轴的坐标系中，每个门限求得的虚警概率和检测概率对应于坐标轴中的点，用光滑曲线连接这些点得到ROC，如图6所示，可以看出和两航过相比，三航过相干变化检测在相同的虚警概率的条件下，检测概率明显提高。 

Claims (5)

1.基于广义似然比的多航过SAR相干变化检测方法，其特征在于：包括以下步骤： 

S1：选取多航过SAR图像对：选取不同时间对同一地区多次观测获得的K幅SAR图像并进行配准后记为{f₁,f₂,…,f_K}，K≥3； 

S2：选取多航过SAR图像像素对，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的值，并记为向量假设向量X_mn～CN(0,Γ)，其中Γ＝E[X_mnX_mn ^H]为协方差矩阵，设H₀表示目标区域未发生变化，得到协方差矩阵Γ₀，设H₁表示目标区域发生变化，得到协方差矩阵Γ₁； 

S3：对协方差矩阵进行最大似然估计：分别对协方差矩阵Γ₀和Γ₁进行最大似然估计； 

S4：进行似然比假设检验； 

S5：依次选取多航过SAR图像像素对，重复步骤S4，得到变化检测结果。 

2.根据权利要求1所述的基于广义似然比的多航过SAR相干变化检测方法，其特征在于：所述的步骤S2的具体实现方法为：由于图像存储为矩阵形式，因此f₁,f₂,…,f_K均为复图像，依次选取f₁,f₂,…,f_K对应的复矩阵第m行第n列的值分别记为令复向量m＝1,2,…,M，n＝1,2,…,N，其中[]T表示转置运算，M、N为图像对应的复矩阵的大小； 

设向量X_mn服从K维圆对称复高斯分布，即X_mn～CN(0,Γ)，其中协方差矩阵Γ＝E[X_mnX_mn ^H]，因此X_mn的概率密度函数表示为： 

其中，E[X_mnX_mn ^H]表示求X_mnX_mn ^H的均值，X_mn ^H表示求X_mn的复共轭转置，|Γ|表示Γ的行列式，exp表示指数运算； 

设H₀表示目标区域未发生变化，此时对应的向量X_mn～CN(0,Γ₀)，其中， 

其中，j表示复数单元，ρ_ab≈1，Φ_ab≈0°，a＝1，2，…，K，b＝1，2，…，K，a＜b； 

设H₁表示目标区域发生变化，此时对应的向量X_mn～CN(0,Γ₁)，其中， 

ρ′_ab≈0，Φ′_ab≠0°，由于Γ₀、Γ₁未知，因此为复合假设。 

3.根据权利要求2所述的基于广义似然比的多航过SAR相干变化检测方法，其特征在于：所述的步骤S3对协方差矩阵Γ₀和Γ₁进行最大似然估计的具体方法为： 

S31：对Γ₀进行最大似然估计，根据最大似然估计理论，对于假设H₀，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显未变化区域得到参数为Γ₀的复高斯随机向量Y的L个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₀)，Γ₀的似然函数lik(Γ₀)表示为： 

令求得协方差矩阵Γ₀的最大似然估计： 

S32：对Γ₁进行最大似然估计，对于假设H₁，根据多航过SAR图像对{f₁,f₂,…,f_K}中明显变化区域得到的参数为Γ₁的复高斯随机向量Y的S个相互独立的观测Y_i，每个观测Y_i的概率密度函数为p(Y_i|Γ₁)，Γ₁的似然函数lik(Γ₁)可以表示为： 

令求出协方差矩阵Γ₁的最大似然估计： 

4.根据权利要求3所述的基于广义似然比的多航过SAR相干变化检测方法，其特征 在于：所述的步骤S4包括以下子步骤： 

S41：确定检验统计量，根据步骤S2的假设和步骤S3的参数估计得到： 

取K幅多航过SAR图像中第m行第n列像素点的Q个相互独立的像素对根据似然比假设检验理论，令似然比： 

将代入上式化简得， 

取对数并忽略常数项，得检测统计量： 

其中表示矩阵G的迹，也就是矩阵对角元素的和， 

S42：进行假设检验，定义一个和用来检测的图像一样大小的矩阵，记为R，选择门限T，判断门限T与检测统计量Z的大小：当Z>T时，判定假设H₁成立，变化检测结果就是该像素对应的区域发生变化，令R对应的像素值为255；否则判定假设H₀成立，变化检测结果就是该像素对应的区域未发生变化，令R对应的像素值为0。 

5.根据权利要求4所述的基于广义似然比的多航过SAR相干变化检测方法，其特征在于：所述的步骤S5中得到变化检测结果的方法为：依次选取多航过SAR图像像素对，重复步骤S4，直到确定矩阵R中所有的像素值，R即为变化检测结果。