CN102496154A - 一种基于马尔科夫随机场的多时相遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于马尔科夫随机场的多时相遥感图像变化检测方法，属于数字图像处理与理解领域，目的在于自动检测已配准的多时相遥感图像中的变化与非变化区域，以提高遥感图像变化检测中的自动化程度和精度。该方法包括差值图像生成步骤、EM参数估计步骤、差值图像边缘检测步骤、自适应权重计算步骤以及马尔科夫随机场标记步骤。本发明利用差值图像的相邻像素大小自动调整马尔科夫随机场的权值大小，有效提高变化检测精度。

Description

一种基于马尔科夫随机场的多时相遥感图像变化检测方法

技术领域

本发明属于图像变化检测技术领域，具体涉及一种基于马尔科夫随机场的多时相遥感图像变化检测方法，用于自动检测已配准的多时相遥感图像中的变化与非变化区域，以提高遥感图像变化检测中的自动化程度和精度。

背景技术

图像变化检测技术是基于计算机图像处理系统，对不同时段目标或现象状态发生的变化进行识别分析，包括变化检测、变化定位、变化区分等关键技术，是目前数字图像处理与理解领域的前沿分支。它广泛应用于多个领域，例如土地使用监测、森林监测、农业勘察、城市研究等。

马尔科夫随机场(Markov Random Field，MRF)方法是主流的遥感图像变化检测方法之一。马尔科夫随机场(Markov Random Field，MRF)方法在[L.Bruzzone and P.D.Fernandez，Automatic analysis of the difference image forunsupervised change detection.IEEE transactions on Geoscience and RemoteSensing[J].2000，38(3)：1171-1182]中提出。MRF考虑包含于差值图像中的空间上下文信息，增加了最终变化检测结果的精确度。然而，传统的MRF模型，假定邻域内像素标记具有一致性。因此，邻域像素的权重被设定为相同常数。也就是说，传统MRF模型支持邻域内出现最多的标记。这对于平坦光滑、没有较多细节的区域是恰当的。因为在这些光滑区域中，像素标记的不一致性主要由于噪声，例如配准误差、传感器噪声等。但是，在细节区域中，传统的MRF模型经常得到错误标记。

发明内容

为了克服现有MRF变化检测方法的缺点，本发明提出了一种改进的基于马尔科夫随机场的多时相遥感图像变化检测方法，提高了检测精确度。

一种基于马尔科夫随机场的多时相遥感图像变化检测方法，包括：生成差值图像的步骤、采用EM方法估计差值图像的图像参数的步骤、提取差值图像的边缘像素的步骤、依据像素边缘状态设定差值图像内各像素权值的步骤以及依据图像参数和像素权值对差值图像作马尔科夫随机场类别标记的步骤；其特征在于，

所述设定差值图像各像素权值的步骤具体为：设定差值图像的非边缘像素的权值为预定常数值β_max，边缘像素的内各像素的权重为

其中，η为两相邻像素(i，j)与(g，h)的亮度差，f_γ(η)为自适应权重函数，F_(ij，gh)(η)为(i，j)与(g，h)之间的权重，

为

的反状态；

进一步地，所述马尔科夫随机场类别标记的步骤具体为：

求解贝叶斯准则方程式得到初始阈值T₀；

以T₀为阈值对差值图像作二值化处理：

迭代更新像素类别C_l(i，j)＝arg min{U[C_l(i，j)|X_D(i，j)]}，其中，U[C_l(i，j)|X_D(i，j)]＝U_data[X_D(i，j)|C_l(i，j)]+U_context[C_l(i，j)|{C_l(g，h)}]，

$U_{\mathrm{data}}\left[X_D(i,j)\right|C_l(i,j)]=\frac{1}{2}\ln |{2\mathrm{\π\σ}}_{C_l}^2|+{[X_D(i,j)-{\mathrm{\μ}}_{C_l}]}^2{\left[{\mathrm{\σ}}_{C_l}^2\right]}^{-1},$

$U_{\mathrm{context}}\left[C_l(i,j)\right|C_l(g,h),(g,h)\∈N(i,j)]=\underset{(g,h)\∈N(i,j)}{\mathrm{\Σ}}{F}_{(\mathrm{ij},\mathrm{gh})}\left(\mathrm{\η}\right){\mathrm{\δ}}_k[C_l(i,j),C_l(g,h)],$

${\mathrm{\δ}}_k[C_l(i,j),C_l(g,h)]=\left\{\begin{array}{cc}-1& C_l(i,j)=C_l(g,h)\\ 0& C_l(i,j)\≠C_l(g,h)\end{array}\right.,$

X_D(i，j)为像素(i，j)的亮度，N(i，j)为像素(i，j)的邻域，

和

分别图像参数中的均值和方差。

进一步地，迭代次数为30～100。

进一步地，采用8-邻域线过程提取差值图像的边缘像素。

进一步地，所述预定常数值10≤β_max≤100。

进一步地，所述自适应权重函数f_γ(η)满足：

$\left\{\begin{array}{c}{f}_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)=f_{\mathrm{\&gamma;}}(-\mathrm{\&eta;})\\ f_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)>0(\&ForAll;-a\&le;\mathrm{\&eta;}\&le;a,a>0)\\ {f^{\&prime;}}_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)<0(\&ForAll;\mathrm{\&eta;}>0)\\ {\lim }_{\mathrm{\&eta;}\&RightArrow;0}{f}_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)=1\\ {\lim }_{\mathrm{\&eta;}\&RightArrow;\&infin;}{f}_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)\&le;-\&infin;\end{array}\right.,$

a为f_γ(η)与X轴交点的X轴坐标，f′_γ(η)表示对f_γ(η)求导。

进一步地，所述自适应权重函数为f_1γ(η)＝exp(-η²/γ)、f_2γ(η)＝1/(1+η²/γ)²、f_3γ(η)＝1/(1+η²/γ)、f_4γ(η)＝1/(1+|η|/γ)、f_5γ(η)＝1-|η|/γ、f_6γ(η)＝2exp[-(η/γ)²]-1、f_7γ(η)＝2/[1+(η/γ)²]-1或f_8γ(η)＝2/(1+|η|/γ)-1。

本发明提供一种改进的基于MRF的多时相遥感图像变化检测方法，通过提取差值图像的边缘点，并且根据边缘点与邻域像素的亮度差值自适应调整MRF先验能量的权重，克服了传统MRF变化检测方法细节保持能力差的缺陷，提高了遥感图像变化检测中的自动化程度和精度。

附图说明

图1为本发明步骤流程图；

图2为原始遥感图像时相1示意图；

图3为原始遥感图像时相2示意图；

图4为参考变化区域图像；

图5为检测到的边缘点图像；

图6为最终变化检测结果图像。

具体实施方式

一种基于马尔科夫随机场的多时相遥感图像变化检测方法，用于自动检测已配准的多时相遥感图像中的变化与非变化区域，以提高遥感图像变化检测中的自动化程度和精度。其具体步骤为：

(1)差值图像生成步骤

选取两幅同一地理位置、不同时间取得的已配准的遥感图像，设它们具有同样大小M×N，分别记为X₁(i，j)和X₂(i，j)，1≤i≤I，1≤j≤J，将对应的像素分别相减得到差值图像记为X_D(i，j)。即表示为：

X_D(i，j)＝X₁(i，j)-X₂(i，j)

(2)EM参数估计步骤

假定差值图像中的变化和非变化类像素亮度服从高斯分布，则总体亮度直方图为混合高斯分布。使用EM参数估计得到两类各自的均值和方差记为：Θ＝{μ_c，σ_c，μ_n，σ_n}，下标n表示非变化类，下标c表示变化类。

非变化类的迭代步骤为：

$P^{t+1}\left({\mathrm{\&omega;}}_n\right)=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_n)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}$

${\mathrm{\&mu;}}_n^{t+1}=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_n)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_n)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\&sigma;}}_n^2\right)}^{t+1}=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_n)}{p^t\left(X(i,j)\right)}{[X(i,j)-{\mathrm{\&mu;}}_n^t]}^2}{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_n\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_n)}{p^t\left(X(i,j)\right)}}$

变化类的迭代步骤为：

$P^{t+1}\left({\mathrm{\&omega;}}_c\right)=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_c)}{p^t\left(X(i,j)\right)}}{\mathrm{IJ}}$

${\mathrm{\&mu;}}_c^{t+1}=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_c)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_c)}{p^t\left(X(i,j)\right)}}$

${\left({\mathrm{\&sigma;}}_c^2\right)}^{t+1}=\frac{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_c)}{p^t\left(X(i,j)\right)}{[X(i,j)-{\mathrm{\&mu;}}_c^t]}^2}{\underset{X(i,j)\&Element;X_D}{\mathrm{\&Sigma;}}\frac{P^t\left({\mathrm{\&omega;}}_c\right)p^t\left(X(i,j)\right|{\mathrm{\&omega;}}_c)}{p^t\left(X(i,j)\right)}}$

其中t＝0，1，2…为迭代次数。迭代公式中的

$p\left(x\right|{\mathrm{\&omega;}}_i)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_i^2}}\exp [-\frac{{(x-{\mathrm{\&mu;}}_i)}^2}{{2\mathrm{\&sigma;}}_i^2}]$ i∈{n，c}

p(x)＝p(x|ω_n)P(ω_n)+p(x|ω_c)P(ω_c)

其中，ω_n和ω_c分别表示非变化类和变化类，p(x|ω_i)和p(x)分别为似然概率和全概率。

初始像素类别集合的阈值选为：

T_n＝M_D(1-α)

T_c＝M_D(1+α)

其中M_D＝[max(X_D)-min(X_D)]/2为像素中值，α∈(0，1)。由阈值T_n、T_c，可将初始变化、非变化类集合分别选择为：

S_c＝{X(i，j)|X(i，j)＞T_c}

S_n＝{X(i，j)|X(i，j)＜T_n}

因此，迭代公式中的先验概率、均值和方差初始值为：

$P^0\left({\mathrm{\&omega;}}_c\right)=\frac{\left|\right|S_c\left|\right|}{I\&CenterDot;J},$ $P^0\left({\mathrm{\&omega;}}_n\right)=\frac{\left|\right|S_n\left|\right|}{I\&CenterDot;J}$

${\mathrm{\&mu;}}_n^0=\frac{\underset{X(i,j)\&Element;S_n}{\mathrm{\&Sigma;}}X(i,j)}{\left|\right|S_n\left|\right|},$ ${\mathrm{\&mu;}}_c^0=\frac{\underset{X(i,j)\&Element;S_c}{\mathrm{\&Sigma;}}X(i,j)}{\left|\right|S_c\left|\right|}$

${\left({\mathrm{\&sigma;}}_n^2\right)}^0=\frac{\underset{X(i,j)\&Element;S_n}{\mathrm{\&Sigma;}}{[X(i,j)-{\mathrm{\&mu;}}_n^0]}^2}{\left|\right|S_n\left|\right|},$ ${\left({\mathrm{\&sigma;}}_c^2\right)}^0=\frac{\underset{X(i,j)\&Element;S_c}{\mathrm{\&Sigma;}}{[X(i,j)-{\mathrm{\&mu;}}_c^0]}^2}{\left|\right|S_c\left|\right|}$

其中||·||表示集合中元素的个数。

(3)差值图像边缘检测步骤。

边缘提取可采用Canny算子、Sobel算子、LOG算子等方法，本发明优选8-邻域线过程从差值图像中提取出边缘像素。使用以下记号：

S^P＝{i|1≤i≤m}

S^E＝{(i，i′)|i∈S^P，i′∈S^P，i′∈N_i}

分别表示亮度点集合和边缘网格，m为差值图像像素总数，N_i为i的邻域。

先验能量定义为：

$U(f^P,f^E)=\underset{i\&Element;S^P}{\mathrm{\&Sigma;}}\underset{i^{\&prime;}\&Element;N_i}{\mathrm{\&Sigma;}}g(f_i^P,f_{i^{\&prime;}}^P,f_{i,i^{\&prime;}}^E)$

$=\underset{i\&Element;S^P}{\mathrm{\&Sigma;}}\underset{i^{\&prime;}\&Element;N_i}{\mathrm{\&Sigma;}}[{(f_i^P-f_{i^{\&prime;}}^P)}^2(1-f_{i,i^{\&prime;}}^E)+{\mathrm{\&alpha;f}}_{i,i^{\&prime;}}^E]$

其中

为势函数，

和分别表示像素i和i’的估计值，

表示像素i和i’的边缘状态。惩罚参数α控制多余边缘数量，本发明设为α＝(μ_c-μ_n)²。

后验能量定义为：

$U(f^P,f^E|d)=\underset{i\&Element;S^P}{\mathrm{\&Sigma;}}{(f_i^P-d_i)}^2$

$+\underset{i\&Element;S^P}{\mathrm{\&Sigma;}}\underset{i^{\&prime;}\&Element;N_i}{\mathrm{\&Sigma;}}\text{\&lambda;g}(f_i^P,f_{i^{\&prime;}}^P,f_{i,i}^E)$

其中d为差值图像观测值，d_i为差值图像像素i的观测值，λ为平滑参数，λ可设为1≤λ≤100。边缘场f^E通过对于(f^P)^*阈值化得到：

$f_{i,i}^E=\left\{\begin{array}{cc}1& \mathrm{if}{(f_i^P-f_{i^{\&prime;}}^P)}^2>\mathrm{\&alpha;}\\ 0& \mathrm{otherwise}\end{array}\right.$

(4)自适应权重计算步骤。

对于非边缘点，其权值不需要调整，而直接设置为高的常数值β_max，可设为10≤β_max≤100。对于边缘点，根据其自身与邻域像素的灰度差值自适应调整其权值，即为β_maxf_γ(η)。其中f_γ(η)为自适应权重函数，η为差值图像相邻像素与当前考虑像素的亮度差值。即：

$F_{(\mathrm{ij},\mathrm{gh})}\left(\mathrm{\&eta;}\right)={\mathrm{\&beta;}}_{\max }{f}_{\mathrm{\&gamma;}}{f}_{(i,j)}^E+{\mathrm{\&beta;}}_{\max }{\stackrel{\&OverBar;}{f}}_{(i,j)}^E$

$={\mathrm{\&beta;}}_{\max }[f_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)f_{(i,j)}^E+{\stackrel{\&OverBar;}{f}}_{(i,j)}^E]$

其中η为相邻两像素(i，j)与(g，h)的亮度差，f_γ(η)为自适应权重函数，γ为f_γ(η)的一个参数，F_(ij，gh)(η)为(i，j)与(g，h)之间的权重。

为点(i，j)的边缘状态，

为

的反状态。

定义：以γ(＞0)为参数的自适应权重函数(AWF)f_γ满足：

(i)f_γ(η)＝f_γ(-η)

$\left(\mathrm{ii}\right)f_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)>0(\&ForAll;-a\&le;\mathrm{\&eta;}\&le;a,a>0)$

$\left(\mathrm{iii}\right){f^{\&prime;}}_{\mathrm{\&gamma;}}\left(\mathrm{\&eta;}\right)<0(\&ForAll;\mathrm{\&eta;}>0)$

(iv)lim_η→0f_γ(η)＝1

(v)lim_η→∞f_γ(η)≤-∞

其中η为当前像素与邻域像素的差值，a为f_γ(η)与x轴的交点横坐标，可由γ控制，可取1≤γ≤1000，本发明设γ＝1。

列举8个AWF为：

f1γ(η)＝exp(-η2/γ)	f5γ(η)＝1-|η|/γ
		f2γ(η)＝1/(1+η2/γ)2	f6γ(η)＝2exp[-(η/γ)2]-1
f3γ(η)＝1/(1+η2/γ)	f7γ(η)＝2/[1+(η/γ)2]-1

f4γ(η)＝1/(1+|η|/γ)

f8γ(η)＝2/(1+|η|/γ)-1

(5)马尔科夫随机场标记步骤。根据上一步的自适应权重，可以将修改后的马尔科夫随机场模型应用于差值图像。迭代最小化马尔科夫随机场能量函数表达式可得到最终变化检测结果。

其具体步骤为：

(5.1)根据EM算法得到的参数Θ＝{μ_c，σ_c，μ_n，σ_n}，两类的先验概率P(ω_n)和P(ω_c)求解下述贝叶斯准则方程式得到初始阈值T₀：

$({\mathrm{\&sigma;}}_n^2-{\mathrm{\&sigma;}}_c^2)T_0^2+2({\mathrm{\&mu;}}_n{\mathrm{\&sigma;}}_c^2-{\mathrm{\&mu;}}_c{\mathrm{\&sigma;}}_n^2)T_0$

$+{\mathrm{\&mu;}}_c^2{\mathrm{\&sigma;}}_n^2-{\mathrm{\&mu;}}_n^2{\mathrm{\&sigma;}}_c^2-2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_c^2\ln \frac{{\mathrm{\&sigma;}}_{n}P\left({\mathrm{\&omega;}}_c\right)}{{\mathrm{\&sigma;}}_{c}P\left({\mathrm{\&omega;}}_n\right)}=0$

(5.2)根据初始阈值T₀将差值图像初始划分为两类，即：

(5.3)迭代更新。对差异图像中的每个像元，计算该像元的每一类所对应的局部后验能量：

U[C_l(i，j)|X_D(i，j)]＝U_data[X_D(i，j)|C_l(i，j)]

    +U_context[C_l(i，j)|{C_l(g，h)}]

其中马尔科夫随机场似然能量为：

$U_{\mathrm{data}}\left[X_D(i,j)\right|C_l(i,j)]=\frac{1}{2}\ln |{2\mathrm{\&pi;\&sigma;}}_{C_l(i,j)}^2|$

$+\frac{1}{2}{[X_D(i,j)-{\mathrm{\&mu;}}_{C_l(i,j)}]}^2{\left[{\mathrm{\&sigma;}}_{C_l(i,j)}^2\right]}^{-1}$

其中X_D(i，j)为差值图像中像素(i，j)的亮度值，C_l(i，j)为像素(i，j)的类别标记。

和

分别为通过EM参数估计得到的图像参数中的均值和方差(包括变化类和非变化类)。

马尔科夫随机场先验能量为：

$U_{\mathrm{context}}\left[C_l(i,j)\right|C_l(g,h),(g,h)\&Element;N(i,j)]$

$=\underset{(g,h)\&Element;N(i,j)}{\mathrm{\&Sigma;}}{F}_{(\mathrm{ij},\mathrm{gh})}\left(\mathrm{\&eta;}\right){\mathrm{\&delta;}}_k[C_l(i,j),C_l(g,h)]$

其中δ_k为delta函数，表示为：

${\mathrm{\&delta;}}_k[C_l(i,j),C_l(g,h)]=\left\{\begin{array}{cc}-1& C_l(i,j)=C_l(g,h)\\ 0& C_l(i,j)\&NotEqual;C_l(g,h)\end{array}\right.$

然后找出使后验能量最小的类别标记C′_l(i，j)替换当前类别标记C_l(i，j)：

C′_l(i，j)＝arg min{U[C_l(i，j)|X(i，j)]}

(5.4)如果迭代次数小于预定次数，重复(5.3)步，否则输出最后标记结果。迭代次数取30至100。

Claims (7)

1.一种基于马尔科夫随机场的多时相遥感图像变化检测方法，包括：

生成差值图像的步骤、采用EM方法估计差值图像的图像参数的步骤、提取差值图像的边缘像素的步骤、依据像素边缘状态设定差值图像内各像素权值的步骤以及依据图像参数和像素权值对差值图像作马尔科夫随机场类别标记的步骤；其特征在于，

所述设定差值图像各像素权值的步骤具体为：设定差值图像的非边缘像素的权值为预定常数值β_max，边缘像素的内各像素的权重为 

其中，η为两相邻像素(i，j)与(g，h)的亮度差，f_γ(η)为自适应权重函数，F_(ij，gh)(η)为(i，j)与(g，h)之间的权重， 

为 

的反状态。

2.根据权利要求1所述的多时相遥感图像变化检测方法，其特征在于，所述马尔科夫随机场类别标记的步骤具体为：

求解贝叶斯准则方程式得到初始阈值T₀；

以T₀为阈值对差值图像作二值化处理：

迭代更新像素类别C_l(i，j)＝argmin{U[C_l(i，j)|X_D(i，j)]}，其中，U[C_l(i，j)|X_D(i，j)]＝U_data[X_D(i，j)|C_l(i，j)]+U_context[C_l(i，j)|{C_l(g，h)}]，

X_D(i，j)为像素(i，j)的亮度，N(i，j)为像素(i，j)的邻域， 

和 分别图像参数中的均值和方差。

3.如权利要求2所述的多时相遥感图像变化检测方法，其特征在于，迭代次数为30～100。

4.如权利要求1所述的多时相遥感图像变化检测方法，其特征在于，采用8-邻域线过程提取差值图像的边缘像素。

5.如权利要求1所述的多时相遥感图像变化检测方法，其特征在于，所述预定常数值10≤β_max≤100。

6.如权利要求1所述的多时相遥感图像变化检测方法，其特征在于，所述自适应权重函数f_γ(η)满足：

a为f_γ(η)与X轴交点的X轴坐标，f′_γ(η)表示对f_γ(η)求导。

7.如权利要求1所述的多时相遥感图像变化检测方法，其特征在于，所述自适应权重函数为f_1γ(η)＝exp(-η²/γ)、f_2γ(η)＝1/(1+η²/γ)²、f_3γ(η)＝1/(1+η²/γ)、f_4γ(η)＝1/(1+|η|/γ)、f_5γ(η)＝1-|η|/γ、f_6γ(η)＝2exp[-(η/γ)²]-1、f_7γ(η)＝2/[1+(η/γ)²]-1或f_8γ(η)＝2/(1+|η|/γ)-1。 

Images