CN101833093A - 星载多时相合成孔径雷达图像的自动变化检测方法 - Google Patents

Abstract

本发明公开了一种星载多时相合成孔径雷达图像的自动变化检测方法，涉及合成孔径雷达技术，其包括：1：成像处理后，两个时相获取的合成孔径雷达图像数据的各自标准化算术平均趋近于标准正态分布；2：分别对趋近于正态分布的两个时相获取的合成孔径雷达图像数据的算术平均采用标准正态分布和Hermite多项式逼近；3：采用归一化的交叉熵计算两个时相合成孔径雷达图像的差异值，得到一幅差异指数图像；4：对差异指数图像采用恒虚警率检测方法进行变化区域的自动提取。本发明方法能够有效提取SAR图像上的变化区域，可在多个领域内广泛应用。

Description

星载多时相合成孔径雷达图像的自动变化检测方法

技术领域

本发明涉及合成孔径雷达(SAR)技术领域，是一种星载多时相合成孔径雷达图像的自动变化检测方法，对不同时间获取的卫星SAR图像进行变化检测的方法。

背景技术

随着合成孔径雷达(SAR)的发展，应用多时相的SAR数据进行地表和地物变化的动态监视成为十分重要的遥感应用领域。由于星载SAR遥感数据具有全天候、全天时的独特优势，同时具有固定的重访周期，所以它比光学遥感更适合于变化检测。但是SAR作为主动遥感，其成像机理与光学遥感有本质的差异，数据处理要复杂得多，因此SAR变化检测要困难得多，尤其是SAR固有的相干斑噪声影响，更加重了SAR变化检测的难度。尽管人们对SAR变化检测进行了深入的研究，但主要沿用了光学遥感变化检测方法，检测精度和效果一直不理想。由于SAR变化检测技术在国内起步较晚，加之SAR数据源不足，因此在变化检测的方法和技术上没有根本的突破，误检和漏检现象比较严重，大大影响了检测效果和精度。

发明内容

本发明的目的是提供一种星载多时相合成孔径雷达(SAR)图像的自动变化检测方法，以解决星载SAR图像的变化检测问题，提高检测效果和精度。

为达到上述目的，本发明的的技术方案是：

一种星载多时相合成孔径雷达图像的自动变化检测方法，其包括：

步骤1：成像处理后，根据中心极限定理，两个时相获取的合成孔径雷达图像数据的各自标准化算术平均趋近于标准正态分布；

步骤2：基于Edgeworth展开原理，分别对趋近于正态分布的两个时相获取的合成孔径雷达图像数据的算术平均采用标准正态分布和Hermite多项式逼近；

步骤3：分别对两个时相获取的合成孔径雷达图像在Edgeworth逼近的基础上，采用归一化的交叉熵计算两个时相合成孔径雷达图像的差异值，得到一幅差异指数图像；

步骤4：对差异指数图像采用基于瑞利分布的恒虚警率检测方法进行变化区域的自动提取。

所述的检测方法，其所述步骤3中得到的一幅差异指数图像，显示了两个时相合成孔径雷达图像之间的差异程度，若根据预先设定的阈值对差异指数图像进行分割，即分割成变化的区域和未变化的区域。

所述的检测方法，其所述步骤3中得到的一幅差异指数图像，其直方图分布接近瑞利分布，瑞丽分布的分布密度函数为：

$f\left(x\right)=\frac{x}{b^2}{e}^{\frac{{-x}^2}{{2b}^2}}$

其中b为形状参数，X为自变量，e为常数2.718281828459。

所述的检测方法，其所述步骤4的基于瑞利分布的恒虚警率检测方法，其阈值计算公式为：

$T=\frac{(\sqrt{-2\log \left(P_{\text{fa}}\right)}-\sqrt{\frac{\mathrm{\π}}{2}})}{\sqrt{2-\frac{\mathrm{\π}}{2}}}\mathrm{\δ}+\mathrm{\μ}$

P_fa为给定的虚警率，δ为差异指数图像计算的方差，μ为差异指数图像的均值；π为圆周率常数，3.141592653589793…；

对于差异指数图像，分别计算出均值、方差，然后给定虚警率，代入上式计算出阈值，然后对差异图像进行二值化分割，得变化区域。

本发明方法充分考虑了SAR图像的统计特性，综合应用了Edgeworth展开原理、归一化交叉熵和恒虚警率检测方法，能够有效地进行不同时相SAR图像的变化检测。

本发明方法有助于推动我国遥感动态监测的技术水平提高，促进SAR遥感数据在土地覆盖/变化、农作物长势、海冰/积雪、植被、各种灾害监测、军事目标的动态监视等领域的广泛应用。

附图说明

图1是本发明星载多时相SAR图像变化检测方法的示意框图。

具体实施方式

下面将结合附图1对本发明加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

一种星载多时相合成孔径雷达(SAR)图像的自动变化检测方法，包括：

步骤1：成像处理后，根据中心极限定理，两个时相获取的SAR图像数据的各自标准化算术平均趋近于标准正态分布；

步骤2：基于Edgeworth展开原理，分别对趋近于正态分布的两个时相获取的SAR图像数据的算术平均采用标准正态分布和Hermite多项式逼近；

步骤3：分别在两个时相取的SAR图像在Edgeworth逼近的基础上，采用归一化的交叉熵计算两个时相SAR图像的差异值，得到一幅差异指数图像；

步骤4：对差异指数图像采用基于瑞利分布的恒虚警率检测方法进行变化区域的自动提取。

根据本发明的实施例，所述变化检测方法，是对不同时相获取的卫星SAR图像进行自动变化检测。

成像处理后，SAR图像数据的每个像元都是地面多个小面元后向散射的矢量合成，是多个随机分布变量的综合。而包含一定区域的SAR图像数据是大量像元组成的分布区域，即分布目标。分布目标内也包含大量的散射体，一般没有十分突出的强散射体，其图像回波数据的分布还是比较均匀的，通常认为数据的分布具有一定的统计特性，可以用某种概率分布模型进行描述，即所谓的杂波模型。而实际应用中，不同的场景，如山地、森林、海洋、河流等场景会对应不同的分布模型，即使同一场景采用不同频段获取不同分辨率的SAR图像，其统计分布特性也会有明显的差异，因此为了准确描述SAR图像的统计分布模型，需要根据SAR图像数据自适应地来逼近分布模型。

根据中心极限定理可知，若X为随机变量，它的均值μ，方差σ，则标准化的算术平均Y为：

$Y=\frac{\frac{1}{n}\mathrm{\ΣX}-\mathrm{\μ}}{\frac{\mathrm{\σ}}{\sqrt{n}}}---\left(1\right)$

其中，n为样本的个数，∑表示对n个X求和。

趋近于标准正态分布，而根据Edgeworth级数展开原理可知，如果一个分布偏离标准正态分布不远，则该分布函数可以用标准正态分布和多项式进行逼近，因此变量Y的分布可以用下式逼近：

f(y)＝G(y)(1+v(y))   (2)

而： $v\left(y\right)=\frac{1}{6}{\mathrm{\ρ}}_3{H}_3\left(y\right)+\frac{1}{24}{\mathrm{\ρ}}_4{H}_4\left(y\right)+\frac{1}{72}{{\mathrm{\ρ}}_3}^2{H}_6\left(y\right)---\left(3\right)$

其中，G(y)表示标准正态分布，ρ₃，ρ₄为变量Y的3阶和4阶累计量按公式(5)和(6)进行计算，H₃，H₄，H₆为变量y的Hermite多项式，即：

$\left.\begin{array}{c}{H}_3=y^3-3y\\ H_4=y^4-{6y}^2+3\\ H_6=y^6-{15y}^4+{45y}^2-15\end{array}\right\}---\left(4\right)$

利用变量Y的生矩函数，可以得到，变量Y的3阶和4阶累计量为：

${\mathrm{\ρ}}_3=\frac{K_3\left(x\right)}{\sqrt{n}{\mathrm{\σ}}^3}---\left(5\right)$

${\mathrm{\ρ}}_4=\frac{K_4\left(x\right)}{n{\mathrm{\σ}}^4}-3---\left(6\right)$

其中K₃(x)，K₄(x)为变量X的的3阶和4阶累计量，可以根据X的各阶原点矩计算：

K₃(x)＝a₃-3a₂a₁+2a₁ ³                (7)

K₄(x)＝a₄-4a₃a₁-3a₂ ²+12a₂a₁ ²-6a₁ ⁴   (8)

a_i表示X的i阶原点矩。

SAR变化检测就是要给出两幅图像的差异情况，根据Edgeworth展开可以得到SAR数据的密度分布函数。交叉熵(Kullback-Leibler Divergence)准确地描述了两个密度分布函数的之间的差异，假设两个分布函数分别为：p(x)和q(x)，则它们的差异K(Q|P)可以用下式表示：

$K\left(Q\right|P)=\∫\log \frac{p\left(x\right)}{q\left(x\right)}p\left(x\right)\mathrm{dx}---\left(9\right)$

log表示对数，∫表示积分，dx表示对变量x进行积分。

为了解决交叉熵的不对称问题，可以采用如下的归一化的聚合交叉熵KLD：

KLD＝[K(Q|P)-K(P|Q)]²    (10)

这样保证了它的非负、对称和规范性。对于两幅SAR图像，可以采用上述的归一化交叉熵进行计算，从而实现对两幅SAR图像差异的计算。

假设两幅不同时相的SAR数据，样本个数为n，分别用变量X和

表示，它们的均值和方差分别为：μ，σ和

则变量

和

将分别趋近于正态分布，可以用Edgeworth展开进行逼近表示：

f(y)＝G(y)(1+v(y))    (11)

$f\left(\tilde{y}\right)=G\left(\tilde{y}\right)(1+u\left(\tilde{y}\right))---\left(12\right)$

其中G(y)和

表示正态分布，而：

$v\left(y\right)=\frac{1}{6}{\mathrm{\ρ}}_3{H}_3\left(y\right)+\frac{1}{24}{\mathrm{\ρ}}_4{H}_4\left(y\right)+\frac{1}{72}{{\mathrm{\ρ}}_3}^2{H}_6\left(y\right)---\left(13\right)$

$u\left(\tilde{y}\right)=\frac{1}{6}{\widetilde{\mathrm{\ρ}}}_3{H}_3\left(\tilde{y}\right)+\frac{1}{24}{\widetilde{\mathrm{\ρ}}}_4{H}_4\left(\tilde{y}\right)+\frac{1}{72}{{\widetilde{\mathrm{\ρ}}}_3}^2{H}_6\left(\tilde{y}\right)---\left(14\right)$

ρ₃，ρ₄和分别对两幅图像数据按照前面的公式(5)和(6)进行计算。

将它们带入交叉熵的计算公式，

$\∫f\left(y\right)\log \frac{f\left(y\right)}{f\left(\tilde{y}\right)}\mathrm{dy}=\∫f\left(y)\right)\log \frac{f\left(y)\right)}{G\left(y\right)}\mathrm{dy}+\∫f\left(y)\right)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}+\∫f\left(y)\right)\log \frac{G\left(\tilde{y})\right)}{f\left(\tilde{y}\right)}\mathrm{dy}---\left(15\right)$

公式(15)右边中第1项为负熵

$J_E(f,G)=\∫f\left(y)\right)\log \frac{f\left(y)\right)}{G\left(y\right)}\mathrm{dy}\≈\frac{1}{12}{{\mathrm{\ρ}}_3}^2---\left(16\right)$

公式(15)中右边第2项：

$\∫f\left(y)\right)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}=\∫G\left(y)\right)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}+\∫G(y)v(y)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}---\left(17\right)$

经过推导公式(17)中右边的第1项为：

$\∫G\left(y)\right)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}=\frac{1}{2}[{\mathrm{\β}}^2-2\log \mathrm{\β}-1+{\mathrm{\α}}^2{\mathrm{\β}}^2]---\left(18\right)$

其中：

$\mathrm{\α}=\frac{\mathrm{\μ}-\widetilde{\mathrm{\μ}}}{\mathrm{\σ}},$ $\mathrm{\β}=\frac{\mathrm{\σ}}{\widetilde{\mathrm{\σ}}},$ $\tilde{y}=\mathrm{y\β}+\mathrm{\α\β}$

而公式(18)中右边的第2项

$\∫G\left(y\right)v\left(y\right)\log \frac{G\left(y)\right)}{G\left(\tilde{y}\right)}\mathrm{dy}\≈0$

公式(16)中右边第3项：

$\∫f\left(y)\right)\log \frac{G\left(\tilde{y})\right)}{f\left(\tilde{y}\right)}\mathrm{dy}=\∫G\left(y\right)(1+v\left(y\right))\log \frac{G\left(\tilde{y}\right)}{G\left(\tilde{y}\right)(1+u\left(\tilde{y}\right))}$

$\≈\∫G\left(y\right)u\left(\tilde{y}\right)\mathrm{dy}+\∫G\left(y\right)v\left(y\right)u\left(\tilde{y}\right)\mathrm{dy}$

$\∫G\left(y\right)u\left(\tilde{y}\right)\mathrm{dy}=\∫G\left(y\right)[\frac{1}{6}{\widetilde{\mathrm{\ρ}}}_3{H}_3\left(\tilde{y}\right)+\frac{1}{24}{\widetilde{\mathrm{\ρ}}}_4{H}_4\left(\tilde{y}\right)+\frac{1}{72}{{\widetilde{\mathrm{\ρ}}}_3}^2{H}_6\left(\tilde{y}\right)]\mathrm{dy}=b_1+b_2+b_3$

其中：

$b_1=\frac{1}{6}{\widetilde{\mathrm{\ρ}}}_3\∫G\left(y\right)H_3\left(\tilde{y}\right)\mathrm{dy}=\frac{1}{6}{\widetilde{\mathrm{\ρ}}}_3\∫G\left(y\right)[{(\mathrm{y\β}+\mathrm{\α\β})}^3-3(\mathrm{y\β}+\mathrm{\α\β})]\mathrm{dy}=\frac{K_3\left(\tilde{x}\right)}{6{\widetilde{\mathrm{\σ}}}^3}(3{\mathrm{\β}}^3\mathrm{\α}+{\mathrm{\α}}^3{\mathrm{\β}}^3-3\mathrm{\α})$

$b_2=\frac{1}{24}{\widetilde{\mathrm{\ρ}}}_4\∫G\left(y\right)H_4\left(\tilde{y}\right)\mathrm{dy}=\frac{1}{24}{\widetilde{\mathrm{\ρ}}}_4\∫G\left(y\right)[{(\mathrm{y\β}+\mathrm{\α\β})}^4-6{(\mathrm{y\β}+\mathrm{\α\β})}^2+3]\mathrm{dy}$

$=\frac{1}{24}(\frac{K_4\left(\tilde{x}\right)}{{\widetilde{\mathrm{\σ}}}^4}-3)(3{\mathrm{\β}}^4\mathrm{\δ}-6{\mathrm{\β}}^2\mathrm{\γ}+3)$

$b_3=\frac{1}{72}{{\widetilde{\mathrm{\ρ}}}_3}^2\∫G\left(y\right)H_6\left(\tilde{y}\right)\mathrm{dy}=\frac{1}{72}{{\widetilde{\mathrm{\ρ}}}_3}^2\∫G\left(y\right)[{(\mathrm{y\β}+\mathrm{\α\β})}^6-15{(\mathrm{y\β}+\mathrm{\α\β})}^4+45{(\mathrm{y\β}+\mathrm{\α\β})}^2-15]\mathrm{dy}$

$=\frac{{K_3}^2\left(\tilde{x}\right)}{72{\widetilde{\mathrm{\σ}}}^{6}n}(15{\mathrm{\β}}^6\mathrm{\η}-15{\mathrm{\β}}^4\mathrm{\δ}+45{\mathrm{\β}}^2\mathrm{\γ}-15)$

$\∫G\left(y\right)v\left(y\right)u\left(\tilde{y}\right)\mathrm{dy}=\frac{K_3\left(x\right)K_3\left(\tilde{x}\right)}{{6\mathrm{\σ}}^3{\widetilde{\mathrm{\σ}}}^{3}n}{\mathrm{\β}}^3$

其中：

γ＝α²+1

δ＝α⁴+6α²+3

η＝α⁶+15α⁴+45α²+15

而式中的 $\mathrm{\α}=\frac{\mathrm{\μ}-\widetilde{\mathrm{\μ}}}{\mathrm{\σ}}$

类似的，也可以同样得到

的具体计算公式。将两个计算结果带入公式(9)，就可以得到归一化的交互熵，它表示了两幅图像的差异程度。

由两幅SAR图像得到一幅差异图象，显示了两个时相之间SAR图像的差异程度，根据差异图像采用阈值分割，就可以将差异图像分割成变化的区域和未变化的区域。

基于直方图进行阈值分割是比较通用而有效的方法，但是根据直方图如何确定阈值一直是难点，手动选择，依赖于人为因素，而且要不断进行调整。自动选择是人们希望实现的方向，为此我们引入恒虚警率(CFAR)的方法来进行变化区域的自动分割和提取。经过处理得到的差异图像，其直方图分布十分接近瑞利分布，因此我们采用了瑞利分布模型进行CFAR来检测变化的区域。

瑞丽分布的分布密度函数为：

$f\left(x\right)=\frac{x}{b^2}{e}^{\frac{-x^2}{{2b}^2}}$

其中，b为形状参数，e为常数2.718281828459。基于瑞利分布的CFAR检测的阈值计算公式：

$T=\frac{(\sqrt{-2\log \left(P_{\text{fa}}\right)}-\sqrt{\frac{\mathrm{\π}}{2}})}{\sqrt{2-\frac{\mathrm{\π}}{2}}}\mathrm{\δ}+\mathrm{\μ}$

P_fa为给定的虚警率，δ为差异图像计算的方差，μ为差异图像的均值，π为圆周率常数，为3.141592653589793…。

对于差异图像，分别计算出均值、方差，然后给定虚警率，就可以根据上式计算出阈值，然后对差异图像进行二值化分割：

如果x＞T，则x为变化的像素；否则，x为没有变化的像素。

Claims (4)

1.一种星载多时相合成孔径雷达图像的自动变化检测方法，其特征在于，包括：

步骤1：成像处理后，根据中心极限定理，两个时相获取的合成孔径雷达图像数据的各自标准化算术平均趋近于标准正态分布；

步骤2：基于Edgeworth展开原理，分别对趋近于正态分布的两个时相获取的合成孔径雷达图像数据的算术平均采用标准正态分布和Hermite多项式逼近；

步骤3：分别对两个时相获取的合成孔径雷达图像在Edgeworth逼近的基础上，采用归一化的交叉熵计算两个时相合成孔径雷达图像的差异值，得到一幅差异指数图像；

步骤4：对差异指数图像采用基于瑞利分布的恒虚警率检测方法进行变化区域的自动提取。

2.如权利要求1所述的检测方法，其特征在于，所述步骤3中得到的一幅差异指数图像，显示了两个时相合成孔径雷达图像之间的差异程度，若根据预先设定的阈值对差异指数图像进行分割，即分割成变化的区域和未变化的区域。

3.如权利要求1所述的检测方法，其特征在于，所述步骤3中得到的一幅差异指数图像，其直方图分布接近瑞利分布，瑞丽分布的分布密度函数为：

$f\left(x\right)=\frac{x}{b^2}{e}^{\frac{-x^2}{{2b}^2}}$

其中，b为形状参数，而x为函数自变量，e为常数2.718281828459。

4.如权利要求1或2所述的检测方法，其特征在于，所述步骤4的基于瑞利分布的恒虚警率检测方法，其阈值T计算公式为：

$T=\frac{(\sqrt{-2\log \left(P_{\text{fa}}\right)}-\sqrt{\frac{\mathrm{\π}}{2}})}{\sqrt{2-\frac{\mathrm{\π}}{2}}}\mathrm{\δ}+\mathrm{\μ}$

P_fa为给定的虚警率，δ为差异指数图像计算的方差，μ为差异指数图像的均值；π为圆周率常数，为3.141592653589793…；

对于差异指数图像，分别计算出均值、方差，然后给定虚警率，代入上式计算出阈值，然后对差异图像进行二值化分割，得到变化区域。

Images