CN1760888A - 利用星载sar多时相图像识别地表变化的方法 - Google Patents

Abstract

本发明属空间遥感和空对地观测信息处理技术领域，具体为一种利用星载SAR多时相图像识别地表变化的方法。本发明用不同时间SAR的图像差值，通过求解双阈值期望极大化和马尔科夫随机场EM－MRF算法，识别不同时间的地表发生散射增强、减弱与不变三类变化，用1996和2002年上海市的欧洲遥感卫星ERS－2 SAR图像作了实例研究，证明本发明方法的有效性。

Description

利用星载SAR多时相图像识别地表变化的方法

发明领域

本发明属于空间遥感与空对地观测信息处理技术领域，具体涉及一种利用星载SAR多时相图像识别地表变化的方法。

背景技术

城市与城市群体的发展是本世纪人类文明发展的重大内容。城市建筑、绿地、道路的发展，水体、植物生态、气候的变化，人口的迁移等等，所引发产生的城市地域与内涵的变化都需要科学的评估。多时相(不同时间)的空间遥感观测为监测与评估这一发展提供了便捷可行的技术手段，它的一个直接应用就是地表变化信息的识别与分类(Rignot andvan Zyl 1993，Singh 1989，Bruzzone and Prieto 2000，Merril and Jiajun 1998，Kasetkasem andVarshney 2002)。全天候高分辨率SAR成像近二十年的发展与空间遥感计划(如SIR-CSAR，ERS-1，2，RADARSAT-1，2)的执行能提供多年观测的数据图像，为地表变化信息的获取与评估提供了条件。

比较不同时相的海量的遥感图像数据，从直观的图像灰度值变化作人工定性的评估是十分粗糙的，难以准确地确定区域的变化，判断的阈值不可靠，极为容易失去重要的信息。因此，如何自动获取与分析地表变化信息，一直是遥感研究的重要课题。

地表变化信息的检测与分析有两类方法：有监督的方法和无监督的方法。前者需要多时相的大量实况数据作为分类过程的训练样本，后者则不需要附加信息，直接对多时相的遥感图像进行分析。本发明采用无监督方法。

Rignot和van Zyl(1993)曾用二幅不同时间的SAR图像之间的差值作变化分析。Singh(1989)用红外多通道数据进行变化矢量分析，由各对应像素计算得到的变化矢量幅值的统计来检测变化。这些方法相对简单易用，但不能对差值图像进行自动或非启发式的分析，对差值图像设定的阈值通常是基于经验或人为的，而阈值选取适当与否在很大程度上决定了检测结果的可靠与准确。

Bruzzone和Prieto(2000)提出了基于Bayes理论和马尔科夫随机场(MRF，MarkovRandom Fields)模型的自动检测区域变化的方法，通过期望极大化(EM，ExpectationMaximum)算法来自动获取一个阈值。他们将差值图像分成两类：变化区和无变化区。但是，他们并没有用真实SAR图像数据来说明这一算法的可靠性与可行性。同时，即使是变化区，也可以有完全不同的变化。比如楼群变成草地，后向散射可能变弱；草地上有树或变成双向反射增强的楼群，后向散射可能增强。后向散射同一db的增强或减弱完全可以有不同的地表变化的信息内涵。随着多通道、多极化、多类遥感器多时相SAR图像的获取，自动检测、识别、分类与评估地表变化更多的信息一定会成为重要的研究方向。

发明内容

本发明的目的在于提出一种利用星载SAR多时相图像正确、有效识别地表变化的方法。

本发明提出的利用星载SAR多时相图像识别地表变化的方法，是一种基于星载SAR多时相微波遥感图像的EM-MRF方法。其基本步骤为：

(1)选取不同时间的两幅星载SAR微波遥感图像，计算其差值图像，该差值图像是两幅图像的相减，或者是两张图像的相除；

(2)采用EM算法将差值图像分成三类：增强区、减弱区和不变区，并获得区分的两个阈值；

(3)根据图像上相邻像素的空间结构相关性，将EM算法的结果作为初值，利用MRF-ICM算法，进行迭代计算，使差值图像的Gibbs能量达到局部最小值，从而得到关于城市变化的散射增强、减弱、不变的分类结果。

本发明用1996和2002年上海市的ERS-2SAR图像作了实例研究，用实地情况佐证。本发明同样可应用于我国长江三角洲、珠江流域等城市链发展的自动识别变化信息评估。在多通道、多极化、多时相SAR数据获取后，可很容易地推广到更多分类区域的分类检测。

附图说明

图1.差值图像直方图。其中，(a)为差值图像直方图上阈值划分示意，(b)为实际差值图像归一化直方图上三类分布。

图2.MRF的二阶邻域系统和组系。其中，(a)二阶邻域系统，(b)单像素组系L₁及其权值α，(c)双像素组系L₂及其权值β。

图3.1996年6月4日ERS-2SAR图像。

图4.2002年4月9日ERS-2SAR图像。

图5.二者的差值图像。

图6.不同迭代次数下各参数值(γ＝0.5)。

图7.双阈值EM算法地表变化分类。

图8.MRF-ICM算法自动检测地表变化分类。

图9.上海市变化图。

图10为本发明流程框图。

具体实施方式

1、EM算法估计变化分类的先验、条件概率

考虑不同时间(t₁，t₂)的两张已经对准的大小为I×J图像X₁和X₂。其差值图像X₀＝(X₂-X₁)(也可用X₀＝(X₂/X₁))，设X是一个随机变量，代表了差值图像中I×J个像素的值，X_D＝{X(i，j)，l≤i≤I，l≤j≤J}。I、J分别为图像的高度和宽度。

检测与分析差值图像是要区分各像素属于不变化类ω_n或变化类ω_c。差值图像X_D的像素值的概率分布函数p(X)包含了这两类的条件概率分布，即有

         p(X)＝p(X|ω_n)P(ω_n)+p(X|ω_c)P(ω_c)                (1)这里P(ω_n)，P(ω_c)是“未变化”和“变化”两类的先验概率分布函数，p(X|ω_k)(k＝n，c)是条件概率分布函数。以Bruzzon和Prieto(2000)的工作作为用期望极大化(EM)算法的基础，对p(X|ω_n)、p(X|ω_c)、p(ω_n)和P(ω_c)进行无监督估计。

EM算法是估计不完整性数据问题的最大似然估计法的一种普遍的方法(Takajo andTakahashi 1998，Ghiglia and Romero 1994，Pritt and Shipman 1994，Pritt 1996)。它包括了一个估计步骤和一个极大化步骤，两个步骤都需要迭代到收敛。估计步骤用受到观测值控制的参数估计作为未知变量来进行计算。极大化步骤提供了对参数的新的估计。假设p(X|ω_n)和p(X|ω_c)服从高斯(Gauss)分布，ω_n对应的概率密度函数由均值μ_n和方差σ_n ²来描述，同样，与ω_c相关的概率密度函数由均值μ_c和方差σ_c ²描述。可以证明(Render andWalker 1984)，估计ω_n的各参数与先验概率分布有公式如下：

$P^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}\·\·\·\left(2\right)$

${\mathrm{\μ}}_n^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j){\∈X}_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}\·\·\·\left(3\right)$

${\left({\mathrm{\σ}}_n^2\right)}^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t(X(i,j)}{[X(i,j)-{\mathrm{\μ}}_n^t]}^2}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}\·\·\·\left(4\right)$

上标t和t+1代表当前t和下一次t+1的迭代。类似地，这些方程也用来估计ω_c的先验概率和条件概率密度函数的均值μ_c和方差σ_c ²。即将式(2)、(3)、(4)中的n改为C即可，依次记分(2)’、(3)’、(4)’。

$P^{t+1}\left({\mathrm{\ω}}_c\right)\frac{\underset{X(i,j){\∈X}_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}\·\·\·{\left(2\right)}^{,}$

${\mathrm{\μ}}_c^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}\·\·\·{\left(3\right)}^{,}$

${\left({\mathrm{\σ}}_c^2\right)}^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t(X(i,j)}{[X(i,j)-{\mathrm{\μ}}_c^t]}^2}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}\·\·\·{\left(4\right)}^{,}$

参数估计先从其初始值开始通过迭代上述公式到收敛来获得。可以证明(Takajo andTakahashi 1998，Ghiglia and Romero 1994)，每迭代一次，估计的参数获得一个对数函数形式L(θ)＝ln p(X_D|θ)的增加量，其中 $\mathrm{\θ}=\{P\left({\mathrm{\ω}}_n\right),P\left({\mathrm{\ω}}_c\right),{\mathrm{\μ}}_n,{\mathrm{\μ}}_c,{\mathrm{\σ}}_n^2,{\mathrm{\σ}}_c^2\},$ 收敛时将获得对数函数形式的本地最大值。

参数估计的初始值可以通过分析差值图像直方图h(X)来获得。在直方图上确定两个阈值T_n和T_c来获得一个比较接近ω_n的像素子集S_n和另一个比较接近ω_c的像素子集S_c。我们将T_n和T_c写成T＝M_D(1-γ)和T＝M_D(1+γ)，其中M_D是h(X)的中值。(即M_D＝[max{X_D}+min{X_D}]/2)，γ∈(0，1)是一个用来定义那些不容易识别的像素在M_D周围变化范围的初始化参数。那么集合S_n＝{X(i，j)|X(i，j)＜T_n}和S_c＝{X(i，j)|X(i，j)＞T_c}就可以用来计算两类ω_n和ω_c各自的初始统计参数。

在各像素值相互独立的假设下，由Bayes最小误差原理，要使差值图像X_D中的像素X(i，j)属于类别ω_k则须使得其后验条件概率最大，即写为

${\mathrm{\ω}}_k=\arg \underset{{\mathrm{\ω}}_i\∈\{{\mathrm{\ω}}_n,{\mathrm{\ω}}_c\}}{\max }\left\{P\left({\mathrm{\ω}}_i\right|X(i,j))\right\}=\underset{{\mathrm{\ω}}_i\∈\{{\mathrm{\ω}}_n,{\mathrm{\ω}}_c\}}{\arg \max \left\{P\left({\mathrm{\ω}}_i\right)p\left(X(i,j)\right|{\mathrm{\ω}}_i)\right\}}\·\·\·\left(5\right)$

这等价于在差值图像上对两类ω_n和ω_c求解最大似然边界T_o。最佳阈值T_o由如下方程给出

$\frac{P\left({\mathrm{\ω}}_c\right)}{P\left({\mathrm{\ω}}_n\right)}=\frac{p\left(X\right|{\mathrm{\ω}}_n)}{p\left(X\right|{\mathrm{\ω}}_c)}\·\·\·\left(6\right)$

在Gauss概率密度分布的假设下，由式(6)可得

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2)T_o^2+2({\mathrm{\μ}}_n{\mathrm{\σ}}_c^2-{\mathrm{\μ}}_c{\mathrm{\σ}}_n^2)T_o+{\mathrm{\μ}}_c^2{\mathrm{\σ}}_n^2-{\mathrm{\μ}}_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^2\ln \left[\frac{{\mathrm{\σ}}_{c}P\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}P\left({\mathrm{\ω}}_c\right)}\right]=0\·\·\·\left(7\right)$

求得T_o：

$T_o=\frac{-b\±\sqrt{b^2-4\mathrm{ac}}}{2a},$

$a={\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2,b=2({\mathrm{\μ}}_n{\mathrm{\σ}}_c^2-{\mathrm{\μ}}_c{\mathrm{\σ}}_n^2),$

$C={\mathrm{\μ}}_c^2{\mathrm{\σ}}_n^2-{\mathrm{\μ}}_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^2\ln \left[\frac{{\mathrm{\σ}}_{c}P\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}P\left({\mathrm{\ω}}_c\right)}\right]\·\·\·{\left(7\right)}^{,}$

2、双阈值的EM算法

SAR图像各像素值表示雷达回波(后向散射)的强弱，两幅不同时间的SAR图像的差值可表明各像素点上后向散射增强、不变和减弱三种情况。增强与减弱表现了地表变化的不同信息。我们将差值图像各点分为三类：增强ω_c1、不变ω_n和减弱ω_c2。在X_D中即分别为正的点集、0周围的点集和负的点集。

于是X_D的概率密度函数p(X)为

      p(X)＝p(X|ω_c1)P(ω_c1)+p(X|ω_n)P(ω_n)+p(X|ω_c2)P(ω_c2)    (8)

采用EM算法，首先对增强类ω_c1和非增强类ω_n1(ω_n和ω_c2)得到阈值T_o1，再对减弱类ω_c2和非减弱类ω_n2(ω_n和ω_c1)得到阈值T_o2，即ω_n1＝ω_n∪ω_c2，ω_n2＝ω_n∪ω_c1。将这两个结果合并得到三种分类。

如图1(a)，直方图h(X)中增大区(x≥0)内用两个阈值T_n1和T_c1来获得一个比较接近ω_n1的像素子集S_n1和另一个比较接近ω_c1像素子集S_c1。将T_n1和T_c1写成T_n1＝M_D1(1-γ)和T_n1＝M_D1(1+γ)，其中M_D1是h(X)中增强区的中值(即M_D1＝[max{X_D}/2])，本文中γ取0.5。同样，在h(X)中减弱区(X＜0)内找到比较接近ω_n2和ω_c2的像素子集S_n2和S_c2，其中M_D1＝[min{X_D}/2]。

图1(b)是以后要用到的图5ERS-2SAR对上海市观测图像实际数据中截取的延中绿地后向散射差值图像的灰度直方图。为便于比较，在图1(b)中用曲线大致标出了三类的概率密度分布。直方图进行归一化，使得直方图可以与概率密度分布相一致。

在增强区内迭代计算时，不管减弱区的像素；在减小区内迭代时，不管增强区的像素。那么得到集合

  S_n1＝{X(i，j)|0＜X(i，j)＜T_n1}，S_c1＝{X(i，j)|X(i，j)＞T_c1}

和S_n2＝{X(i，j)|T_n2＜X(i，j)＜0}，S_c2＝{X(i，j)|x(i，j)＜T_c2}，

用来计算两类ω_n1ω_c1和ω_n2，ω_c2各自相关的初始统计参数。

确定了初始参数后，用式(2～4)进行上述四类的迭代，收敛后得到四类的先验概率密度函数及其统计参数[P(ω_n1)，μ_n1，σ_n1 ²]，[P(ω_c1)，μ_c1，σ_c1 ²]，[P(ω_n2)，μ_n2，σ_n2 ²]和[P(ω_c2)，μ_c2，σ_c2 ²]。再先后求解式(7)，分别得到阈值T_o1和阈值T_o2。

3、Markov随机场(MRF)对空间结构相关的变化分类

Bruzzone和Prieto(2000)使用了两种方法来检测多时相遥感图像上的区域变化，第一种假设差值图像各像素之间是统计独立的，第二种则认为各像素是空间结构相关的。事实上，一个像素属于ω_k类，那么它相邻的像素也很有可能属于ω_k类。使用第二种方法能改善分类的性能，得到更可靠的结果。

用集合C＝{C_l，1≤l≤L}表示差值图像上所有分类结果的集合，

                C_l＝{C_l(i，j)，1≤i≤I，1≤j≤J}，C_l(i，j)∈{ω_c1，ω_c2，ω_n}其中L＝3^IJ，I，J分别为图像的高度和宽度。

根据Bayes最小误差理论，分类结果应使后验概率为最大，即

$C_k=\underset{C_l\∈C}{\arg \max }\left\{P\left(C_l\right|X_D)\right\}=\underset{C_l\∈C}{\arg \max }\left\{P\left(C_l\right)p\left(X_D\right|C_l)\right\}\·\·\·\left(9\right)$

其中X_D为差值图像，P(C₁)是类别标注的先验概率，p(X_D|C_l)是给定类别标注条件下的差值图像的像素值的联合概率密度函数。式(9)的求解需要同时估计P(C_l)和p(X_D|C_l)，计算上会过于复杂。采用Markov随机场(MRF)来建立差值图像上像素间的空间结构相关的模型，问题就可以简化。

为了用MRF来描述问题，需要引进空间邻域系统N。考虑一个坐标为(i，j)的像素，它的邻域可以用下式表达：

           N(i，j)＝{(i，j)+(v，ζ)，(v，ζ)∈N}            (10)

在本发明中，N＝{(±1，0)，(0，±1)，(1，±1)，(-1，±1)}，即一个像素的邻域由它周围的8个像素构成。

用C_l(i，j)表示像素点(i，j)的类别标注，C_lN表示像素点(i，j)的邻域像素的类别标注，C_l|(i，j)表示差值图像X_D中除了像素点(i，j)之外的其余像素的类别标注。把C_l看成一个MRF，则C_l由如下两条性质：

           P(C_l)＞0，_C_l(i，j)∈{ω_c1，ω_c2，ω_n}            (11)

           P(C_l(i，j)|C_l|(i，j))＝P(C_l(i，j)|C_lN)             (12)

MRF的二阶邻域系统和组系如图2所示。

组系是指邻域系统中的包含中心像素点的邻近像素的连通子图，这里考虑单像素和双像素的组系。如图2(b，c)所示，单像素的组系有1种形状，双像素的组系有4种形状，其数学表示分别为：L₁＝(i，j)，L₂＝{{(i，j)，(g，h)}|(g，h)∈N_i(i，j)}。

根据Hammersley-Clifford定理(Besag 1974)，MRF可以等效为Gibbs分布，因此在给定其它像素的类别标注的前提下，坐标(i，j)处的像素的类别标注的条件概率分布可表示为

P(C_l(i，j)|{C_l(g，h)，(g，h)≠(i，j)})＝P(C_l(i，j))|{C_l(g，h)，(g，h)∈N(i，j)})

$=\frac{1}{Z}\exp [-U\left(C_l(i,j)\right|\{C_l(g,h),(g,h)\∈N(i,j)\})]\·\·\·\left(13\right)$

其中U(·)是Gibbs能量函数，Z是归一化因子，

$Z=\underset{(g,h)\∈\{(i,j)\∪N(i,j)\}}{\mathrm{\Σ}}\exp [-U\left(C_l(i,j)\right|\left\{C_l(g,h)\right\})]\·\·\·\left(14\right)$

$U\left(C_l\right|X_D)=\underset{1\≤i\≤I}{\mathrm{\Σ}}\underset{1\≤j\≤J}{\mathrm{\Σ}}U\left(C_l(i,j)\right|C_{l|(i,j)},X_D)=\underset{1\≤i\≤I}{\mathrm{\Σ}}\underset{1\≤j\≤J}{\mathrm{\Σ}}U\left(C_l(i,j)\right|C_{\mathrm{lN}},X_D)\·\·\·\left(15\right)$

在差值图像的直方图为Gauss分布的情况下，

$U\left(C_l(i,j)\right|C_{\mathrm{lN}},X_D)=\frac{1}{2}\ln \left|2\mathrm{\π}{\mathrm{\σ}}_{C_l(i,j)}^2\right|+\frac{{(X(i,j)-{\mathrm{\μ}}_{C_l(i,j)})}^2}{2{\mathrm{\σ}}_{C_l(i,j)}^2}$

${+v}_1\left(C_l(i,j)\right)+\underset{(g,h)\∈N(i,j)}{\mathrm{\Σ}}{v}_2(C_l(i,j),C_l(g,h))---\left(16\right)$

其中v₁(C_l(i，j))＝α，(i，j)∈L₁。这里α为单像素组系的权值，一般为-0.7-0.5之间，实例计算中取α＝0。

$v_2(C_l(i,j),C_l(g,h))=\left\{\begin{array}{cc}-\mathrm{\β},& C_l(i,j)=C_l(g,h),\{(i,j),(g,h)\}\∈L_2\\ \mathrm{\β},& C_l(i,j)\≠C_l(g,h),\{(i,j),(g,h)\}\∈L_2\end{array}\right.\·\·\·\left(17\right)$

这里β是双像素组系的权值，一般在0.7-1.5之间，实例计算中取β＝1.0。

根据式(9)，区域变化图的生成就是给差值图像上各像素点标上类别，使得后验概率最大。在Markov模型下，就是要使式(15)的能量函数最小(Bruzzone and Prieto 2000)，重写式(15)为

$U\left(C_l\right|X_D)=\underset{1\≤i\≤I}{\mathrm{\Σ}}\underset{1\≤j\≤J}{\mathrm{\Σ}}[U_I(X(i,j)/C_l(i,j))$

$+U_C(C_l(i,j)/\{C_l(g,h),(g,h)\∈N(i,j)\left\}\right])\·\·\·\left(18\right)$

其中U_C(·)表示考虑了像素之间类别相关的能量函数，由式(16)的后两项给出。U_I(·)表示假定像素之间类别独立的能量函数，由式(16)的前两项给出，即

$U_I(X(i,j)/C_l(i,j))=\frac{1}{2}\ln \left|2\mathrm{\π}{\mathrm{\σ}}_{C_l(i,j)}^2\right|+\frac{1}{2}{(X(i,j)-{\mathrm{\μ}}_{C_l(i,j)})}^2{\left[{\mathrm{\σ}}_{C_l(i,j)}^2\right]}^{-1}\·\·\·\left(20\right)$

其中， ${\mathrm{\σ}}_{C_l(i,j)}^2\∈\{{\mathrm{\σ}}_{\mathrm{cl}}^2,{\mathrm{\σ}}_{c2}^2,{\mathrm{\σ}}_n^2\},$ ${\mathrm{\μ}}_{C_l(i,j)}\∈\{{\mathrm{\μ}}_{\mathrm{cl}},{\mathrm{\μ}}_{c2},{\mathrm{\μ}}_n\},$ 由EM算法给出。

一般地，式(18)的求解需要用模拟退火等迭代算法，本文采用ICM算法(Besag1986)，此算法计算量小，可以迅速收敛到能量函数的局部极小值。ICM算法的数学表达式可写成

$C_l=\mathrm{ICM}(\mathrm{\α\β}{\mathrm{\μ}}_{C_l(i,j)},{\mathrm{\σ}}_{C_l(i,j)}^2,X_D,C_{l0})\·\·\·\left(21\right)$

其中C_l0为C_l的初始值。

求解(18)得到C_l的步骤如下：

(a)通过EM算法，得到C_l的初始值C_l0。

(b)计算出新的C_l(i，j)，使得式(18)中的能量函数最小。

(c)重复步骤(b)直到收敛。

ICM算法收敛的条件是迭代次数超过30次或少于

的像素点发生变化。ICM算法的具体描述可参见Besag(1986)。

4、城市变化识别实例

图3、4分别为1996年6月4日和2002年4月9日中国上海地区的ERS-2 SAR遥感图像，图像的分辨率是12.5m×12.5m。计算前图像已经校准，并对初始数据进行了3×3的均值滤波以消除噪声起伏。我们从原图中分别切下580像素×980像素的两块区域，中心是上海市的陆家嘴地区。图5是二者的差值图像(也有将二者相除作为差值图像)。差值图像中灰度级大的点，代表该点上的后向散射的差值大。从差值图中可以看到上海市这5年中的变化，但是如何确定变化的分类则是很难自动确定的。

对差值图像进行EM迭代求解，得到各类在迭代过程中先验概率和统计参数在各步骤上的值。第一类迭代过程参见图6(a～c)，第二类迭代过程类似。

在完成EM求解后，分别得到散射增强区和减弱区内两个类的先验概率和概率密度分布的均值和方差，根据公式(7)获得两个阈值T_o1，T_o2用于分类。

经计算，T_o1＝0.946×10^-4，T_o2＝-1.013×10^-4(后向散射差值最大值为5.202×10^-4，最小值为-3.381×10^-4)。

图7是采用EM算法所得的区域变化的检测结果。图中绿色为后向散射增强区，红色代表后向散射减弱区，蓝色代表后向散射变化不大的区。这一分类结果比差值图像更加清晰。

为纳入空间结构相关信息，图8利用基于MRF模型的ICM算法得到结构相关区域变化的自动检测。与图7相比，同一类别的像素点更加集中，孤立的像素点更少，符合MRF模型的理论。比如，图8中奇数(1，3，5)标注的地区为平坦表面后向散射减弱的代表区域，其中1处为陆家嘴绿地，3处为浦东世纪公园，5处为延中绿地。偶数(2，4)标注的地区为多类建筑物产生双向反射和散射增强的代表区域，其中2处为浦东陆家嘴的建筑楼群，4处为上海新国际展览中心。图8所显示的上海市1-5地区的变化与上海市从1996到2002年的建设情况是一致的，说明了本文双阈值EM-MRF检测算法是有效的。将图8的结果覆盖在上海市旅游图上，如图9，是一个有趣的比较。

利用多时相SAR微波遥感图像提出的EM-MRF模型的算法，可自动检测多年间城市地区的变化。这种方法可以对差值图像进行自动分析，不需要作经验估计或人工干预。在假设像素之间统计独立的前提下，按照后向散射变化的强弱，用EM算法将差值图像分成三类：增强区、减弱区和持不变区，获得两个阈值。再进一步根据图像上相邻像素的空间结构相关性，将EM算法的结果作为初始值，利用MRF-ICM算法进行迭代计算，使差值图像的Gibbs能量到达局部最小值，从而得到关于城市变化的散射增强、减弱、不变的分类结果。用1996和2002年的上海地区的ERS-2的多时相遥感图像进行了算法验证，对比实地发展变化，证实了这种方法的有效性。本方法在获得多极化、多通道、多类遥感器多时相遥感图像后，可进一步推广到多阈值的EM-MRF算法，得到更多的信息融合的地表变化信息自动分类。

参考文献

Besag J.(1974)，Spatial interaction and the statistical analysis of lattice systems(withdiscussion).Journal of Royal Statistics Society，Series B，36(2)：192-326.

Besag J.(1986)，On the statistic analysis of dirty pictures.Journal of Royal Statistics Society，Series B，48：259-302.

Bruzzone L.and D.F.Prieto(2000)，Automatic analysis of the difference image forunsupervised change detection.IEEE Transactions on Geoscience and Remote Sensing，38(3)：1171-1182.

Dempster A.P.，N.M.Laird，and D.B.Rubin(1977)，Maximum linkelihood from incompletedata via the EM algorithm.Journal of Royal Statistics Society，39(1)：1-38.

Fukunaga K.(1990)，Introduction to Statistical Pattern Recognition，2^nd ed.London，U.K.：Academic.

Jin Y.Q.，F.Chen and L.Luo(2004)，Automatic analysis of change detection of multi-temporalERS-2SAR images by using two-thresholds EM and MRF algorithms，The Imaging ScienceJournal，52：234-241.

Kasetkasem T.and P.Varshney(2002)，An image change detection algorithm based on Markovrandom field model.IEEE Transactions on Geoscience and Remote Sensing，40(8)：1815-1823.

Merril K.R.and L.Jiajun(1998)，A comparison of four algorithms for change detection in anurban environment.Remote Sensing of Environment，63：95-100.

Moon T.K.(1996)，The expectation-maximization algorithm.Signal Processing Magazine，13(6)：47-60.

Render A.P.and H.F.Walker(1984)，Mixture densities，maximum likelihood and the EMalgorithm.SIAM Review，26(2)：195-239.

Rignot E.and J.van Zyl(1993)，Change detection techniques for ERS-1SAR data.IEEETransactions on Geoscience and Remote Sensing，31(4)：896-906.

Shahshahani B.M.and D.A.Landgrebe(1994)，The effect of unlabeled samples in reducing thesmall size problem and mitigating the Hughes phenomenon.IEEE Transactions onGeoscience and Remote Sensing，32：1087-1095.

Singh A.(1989)，Digtial change detection techniques using remotely-sensed data.InternationalJournal of Remote Sensing，10(6)：989-1003.

Claims (2)

1、一种利用星载SAR多时相图像识别地表变化的方法，其特征在于基本步骤如下：

(1)选取不同时间的两幅星载SAR微波遥感图像，计算其差值图像，该差值图像是两幅图像的相减，或者是两张图像的相除；

(2)采用EM算法将差值图像分成三类：增强区、减弱区和不变区，并获得区分的两个阈值；

(3)将EM算法的结果作为初值，利用MRF-ICM算法，进行迭代计算，使差值图像的Gibbs能量达到局部最小值，从而得到关于城市变化的散射增强、减弱、不变的分类结果。

2、根据权利要求1所述的利用星载SAR多相图像识别地表变化的方法，其特征在于所述两个阈值由下式计算获得：

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2)T_o^2+2({\mathrm{\μ}}_n{\mathrm{\σ}}_c^2-{\mathrm{\μ}}_c{\mathrm{\σ}}_n^2)T_o+{\mathrm{\μ}}_c^2{\mathrm{\σ}}_n^2-{\mathrm{\μ}}_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^2\ln \left[\frac{{\mathrm{\σ}}_{c}P\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}P\left({\mathrm{\ω}}_c\right)}\right]=0---\left(7\right)$

而σ_n、μ_n、p(w_n)由

$P^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}---\left(2\right)$

${\mathrm{\μ}}_n^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(3\right)$

${\left({\mathrm{\σ}}_n^2\right)}^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}{[X(i.j)-{\mathrm{\μ}}_n^t]}^2}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(4\right)$

迭代求得，σ_c、μ_c、p(w_c)由

$P^{t+1}\left({\mathrm{\ω}}_c\right)=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}---{\left(2\right)}^{\′}$

${\mathrm{\μ}}_c^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}---{\left(3\right)}^{\′}$

${\left({\mathrm{\σ}}_c^2\right)}^{t+1}=\frac{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}{[X(i.j)-{\mathrm{\μ}}_c^t]}^2}{\underset{X(i,j)\∈X_D}{\mathrm{\Σ}}\frac{P^t\left({\mathrm{\ω}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\ω}}_c)}{p^t\left(X(i,j)\right)}}---{\left(4\right)}^{\′}$

迭代求得，

其中ω_c为差值图像X_D的不变化类像素，ω_c为差值图像X_D的变化类像素，

p(X)为差值图像X_D的像素值的X的概率分布函数：

                        p(X)＝p(X|ω_n)P(ω_n)+p(X|ω_c)P(ω_c)    (1)

而p(w_n)、p(w_c)为w_n、w_c的先验概率分布函数，p(X|ω_k)(k＝n，c)是条件概率分布函数。

Images