CN103218807A - 一种多光谱遥感影像变化检测方法 - Google Patents

Abstract

本发明涉及遥感影像处理技术领域，提出了一种多光谱遥感影像变化检测方法，主要解决当遥感地物类型较复杂时，传统方法构造的差异影像信息冗余大，伪变化信息及背景噪声干扰严重；传统基于EM算法求得的高斯混合模型（GMM）对差异影像直方图拟合效果差，检测精度低的问题。其实现过程包括下述步骤：首先，采用主成分（PCA）变换与相关系数融合法相结合的方式构造差异影像；其次，利用分支数为k的GMM对差异影像进行建模，并利用MDL-EM算法自适应估计模型各参数；最后，利用基于统计最小错误率的Bayes判别准则确定变化检测的阈值，得到目标影像的变化检测结果。

Description

一种多光谱遥感影像变化检测方法

技术领域

本发明涉及遥感影像处理技术领域，特别涉及一种多光谱遥感影像变化检测方法，适用于同一区域不同时相的多幅含有地物变化信息的多光谱遥感影像数据的变化检测。

背景技术

基于遥感影像的变化检测是根据同一区域不同时相的遥感影像，辨识观测对象或者现象的状态变化过程。它已广泛应用于资源管理与规划、环境保护等诸多领域，并为相关部门提供科学决策的依据。

当前，遥感影像的变化检测问题已经成为一个非常活跃的研究方向，如美国的国家地理空间智能研究所、意大利的Trento大学、中国科学院遥感所、武汉大学等机构和研究小组都在开展相关方面的研究和应用开发，取得了诸多成果。目前遥感影像变化检测方法主要有：代数法、变换法、分类比较法、高级模型法、GIS集成法、视觉分析法和其它方法。

其中代数法以其操作简单、易于实现成为当前变化检测中使用最广泛的方法之一，主要包括差值法、比值法等。代数法的变化检测核心问题为：差异影像的获取和阈值的选取。仅依靠单一方法获取差异影像，并不能较好地反映地面真实变化情况。通过融合的方法来构造差异影像可以克服此问题，但是其不足是该方法主要是针对单一波段遥感影像。

通过差异影像来对感兴趣区的变化和非变化像元做出判定的关键是阈值的选取，常规的手动选择阈值准确性不高，且自动化程度较低。为了克服这些缺点，L.Bruzzone等提出一种基于统计最小错误率的Bayes判别准则的变化检测算法，将差异影像视作由变化和非变化两类像元组成的高斯混合分布，然后利用（Expectation-Maximization，EM）算法估计出高斯混合分布密度函数的参数，进而确定出变化检测的阈值。但当“差异影像”含有多种地物变化类型，光谱信息较复杂时，把差异影像视作仅有两类的高斯混合分布，显然会降低其变化检测精度。另外，EM算法是在已知高斯混合模型（Gaussian Mixture Model，GMM）分支数的情况下对参数进行估计的，而在实际的应用中，这点很难做到。

发明内容

本发明针对上述已有的技术不足，提出一种多光谱遥感影像变化检测方法。该方法能最大限度地保留目标数据变化的细节信息，减少伪变化信息及背景噪声的干扰，准确估计GMM的参数（分支数，权重，均值，方差），实现了对差异影像直方图的最佳拟合，提高了目标影像的变化检测精度。

本发明解决其技术问题所采取的技术方案：一种多光谱遥感影像变化检测方法，包括下述步骤：

步骤1，数据准备：

选取需要进行变化检测的同一区域、不同时相的2幅多光谱遥感影像数据；

步骤2，构造差异影像：

对2幅多光谱遥感影像数据通过相关系数融合法构造差异影像；

差异影像的构造分以下几步进行，首先对经过步骤1得到的2幅多光谱遥感影像数据X₁，X₂，通过ENVI4.8软件分别进行PCA变换，并分别提取其第一主成分X_1pc1，X_2pc1；然后，求取X_1pc1，X_2pc1的差值就影像ΔY₁以及比值影像ΔY₂，最后对获取的差值、比值影像进行融合构造出差异影像；

对于PCA变换，采用ENVI4.8软件对目标数据进行PCA变换时以协方差矩阵的特征向量对应的矩阵作为变换矩阵；

对于n维随机变量X₁,X₂,…,X_n，若c_ij=Cov(X_i,X_j)=E{[X_i-E(X_i)][X_j-E(X_j)]都存在，i,j=1,2,…n,则n维随机变量X₁,X₂,…,X_n的协方差矩阵为：

$C=\left[\begin{array}{cccc}{c}_{11}& c_{12}& ...& c_{1n}\\ c_{21}& c_{22}& ...& c_{2n}\\ ...& ...& ...& ...\\ c_{n1}& c_{n2}& ...& c_{\mathrm{nn}}\end{array}\right]---\left(1\right)$

通过相关系数融合法构造的差异影像F在坐标(p,q)处的像元值为：

F(p,q)=λ_pq[αΔY₁(p,q)+βΔY₂(p,q)]   (2)

式中，λ_pq=ΔY₂(p,q)max(ΔY₂)；α=a|r|+b;α+β=1；r为ΔY₁与ΔY₂的相关系数；其表达式为：

$r=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(\mathrm{\Δ}{Y}_1-\stackrel{\‾}{\mathrm{\Δ}{Y}_1})(\mathrm{\Δ}{Y}_2-\stackrel{\‾}{\mathrm{\Δ}{Y}_2})/\sqrt{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(\mathrm{\Δ}{Y}_1-\stackrel{\‾}{\mathrm{\Δ}{Y}_1})}^2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(\mathrm{\Δ}{Y}_2-\stackrel{\‾}{\mathrm{\Δ}{Y}_2})}^2}---\left(3\right)$

式中，N为GMM的像元总数，

分别为ΔY₁与ΔY₂的均值，a，b为常数；

常数a，b作为权重α的调节因子，a，b∈[0,1]；通过调节a，b的值来调整差值影像和比值影像的权重α,β，取a，b∈{0.1,0.2,…,1},则a，b的值可由下式确定；

$\{a,b\}=\arg \max \left\{\frac{1}{N}\underset{p=1}{\overset{m}{\mathrm{\Σ}}}\underset{q=1}{\overset{n}{\mathrm{\Σ}}}{[(p,q)-\frac{1}{N}\underset{p=1}{\overset{m}{\mathrm{\Σ}}}\underset{q=1}{\overset{n}{\mathrm{\Σ}}}F(p,q)]}^2\right\}---\left(4\right)$

式中，m为差异影像F的行数，n为差异影像F的列数，且有N=m×n；

步骤3，建立模型：

首先，采用GMM对差异影像进行建模，然后，利用EM算法估计GMM的参数，所述参数包括权重，均值，方差；

步骤4，确定变化检测阈值：

利用基于统计最小错误率的Bayes判别准则确定变化检测的阈值，得到多光谱遥感影像的变化检测结果；

所述步骤3包括如下步骤：

首先，采用GMM对差异影像进行建模，然后，利用MDL-EM算法估计GMM的参数；

对于GMM的构造，具体实现为：设f_i为差异影像F的第i个像元灰度值；i=1,…,N，F中存在k个待识别的类，即为k个分支；每个像元之间是相互独立的，且不考虑分支间的相关性，则有：

$p\left(f_i\right)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}p(f_i,j)=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}p\left(f_i\right|j)p\left(j\right)---\left(5\right)$

f_i属于分支j的概率密度函数为：

$p\left(f_i\right|j)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j}\exp [-\frac{{(f_i-{\mathrm{\μ}}_j)}^2}{2{\mathrm{\σ}}_j^2}]---\left(6\right)$

式中：p(j)为GMM中分支j的权重，μ_j，

分别代表分支j的均值和方差；j＝1，2…k;

设：

Θ=(p(1),p(2),...,p(k);θ₁,θ₂,...,θ_k)，则对应的完全数据集的对数似然函数为:

$\ln p\left(f\right|k,\mathrm{\Θ})=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\ln \left[\underset{j=1}{\overset{k}{\mathrm{\Σ}}}p\left(f_i\right|j,\mathrm{\θ})p\left(j\right)\right]---\left(7\right)$

对于参数估计，采用MDL-EM算法来估计所建立的GMM的参数，具体步骤为：设L为差异影像灰度值集合，其一般为0,1,…,255，h(f_i)为直方图，其中f_i为其第i个像元灰度值，则利用EM算法估计GMM的参数权重，均值，方差的迭代公式为：

$p^{t+1}\left(j\right)=\frac{1}{f_{i}h\left(f_i\right)}\underset{f_i\∈L}{\mathrm{\Σ}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}---\left(8\right)$

${\mathrm{\μ}}_j^{t+1}=\frac{\underset{f_i\∈L}{\mathrm{\Σ}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}{f}_i}{\underset{f_i\∈L}{\mathrm{\Σ}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}}---\left(9\right)$

${\left({\mathrm{\σ}}_j^2\right)}^{t+1}=\frac{\underset{f_i\∈L}{\mathrm{\Σ}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}{(f_i-{\mathrm{\μ}}_j^t)}^2}{\underset{f_i\∈L}{\mathrm{\Σ}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}}---\left(10\right)$

t为迭代次数；而GMM最优模型的选择则由作为适应度函数的MDL准则来确定，其表达式为：

$R(k,\mathrm{\Θ})=-\ln P\left(f\right|k,\mathrm{\Θ})+\frac{T\ln \mathrm{NM}}{2}---\left(11\right)$

其中，

k代表GMM的分支数；M代表GMM中像元的维数；

lnP(f|k,Θ)代表差异影像对应的完全数据集的对数似然函数；TlnNM/2代表引入的针对GMM分支数k的惩罚函数。

所述步骤3中的GMM中分支数k通过如下方法确定：

在选定初始最大分支数k_max的情况下，通过引入合并最小距离分支的办法来自适应地依次减少分支数直到k_min，并依据下式，确定最优的分支数k^*和参数Θ^*；

{k^*,Θ^*}=argmin{R(k,Θ),k_min≤k≤k_max}    (12)

利用上述MDL-EM算法便能自适应地估计得到GMM的k^*和Θ^*；

上述合并最小距离分支法，具体实现为：设分支数为k时其任意两类分支l、m，其中1≤l，m≤k，其参数为：p_l,p_m,μ_l,μ_m,

合并后的新分支参数为：p_(l,m)，μ_(l,m)，

首先令μ_l=μ_m=μ_(l,m)；p_l，p_m不变，则在分支数为k时有θ_(l,m)=θ_(l)=θ_(m)∈Θ^(k)，同时注意到因为θ_(l)=θ_(m)所以可将他们视为一个分支，为了便于区分，将Θ^(k-1)中的θ_(l,m)表示成

则有

此时有p_(l,m)=p_l+p_m；再利用θ_(l,m)和

的定义，结合MDL准则有：

$R(k-1,{\mathrm{\θ}}_{{(l,m)}^{-}})-R(k,{\mathrm{\θ}}^{\left(t\right)})$

$=R(k-1,{\mathrm{\θ}}_{{(l,m)}^{-}})-R(k,{\mathrm{\θ}}_{(l,m)})+R(k,{\mathrm{\θ}}_{(l,m)})-R(k,{\mathrm{\θ}}^{\left(t\right)})$

$\≤\frac{1}{2}(1+M+\frac{M(M+1)}{2})\ln N+Q({\mathrm{\θ}}^{\left(t\right)};{\mathrm{\θ}}^{\left(t\right)})-Q({\mathrm{\θ}}_{(l,m)};{\mathrm{\θ}}^{\left(t\right)})---\left(13\right)$

$\≤\frac{1}{2}(1+M+\frac{M(M+1)}{2})\ln N$

$+Q({\mathrm{\θ}}^{\left(t\right)};{\mathrm{\θ}}^{\left(t\right)})-Q({\mathrm{\θ}}^{*};{\mathrm{\θ}}^{\left(t\right)})+Q({\mathrm{\θ}}^{*};{\mathrm{\θ}}^{\left(t\right)})-Q({\mathrm{\θ}}_{(l,m)}^{*};{\mathrm{\θ}}^{\left(t\right)})$

其中θ^*和

为通过EM算法求得的最优参数估计，Q(θ;θ^(t))为根据上一次计算得的结果来估计完全数据集Z的似然函数的期望值，其表达式为：

Q(θ;θ^(t))=E_z[lnp(f,Z|θ)|f,θ^(t)]    (14)

此时有θ^*=θ^(t)，则Q(θ^(t);θ^(t))-Q(θ^*;θ^(t))=0；由EM算法的M步可知：

则对于任意的l、m其合并分支(l,m)的参数经EM算法迭代求得的最优估计为：

p_(l,m)=p_l+p_m    (15)

${\mathrm{\μ}}_{(l,m)}=\frac{p_l{\mathrm{\μ}}_l+p_m{\mathrm{\μ}}_m}{p_l+p_m}---\left(16\right)$

${\mathrm{\σ}}_{(l,m)}^2=\frac{p_l[{\mathrm{\σ}}_l^2+{({\mathrm{\μ}}_1-{\mathrm{\μ}}_{(l,m)})}^2]+p_m[{\mathrm{\σ}}_m^2+{({\mathrm{\μ}}_m-{\mathrm{\μ}}_{(l,m)})}^2]}{p_l+p_m}---\left(17\right)$

利用上式便可将k中任意两个分支l、m合并从而将k降阶到k-1；再定义一个距离函d(l,m)：

$d(l,m)=N{p}_l\{-\frac{1}{2}[1+\ln \left(2\mathrm{\π}\right)]-\frac{1}{2}\ln \left({\mathrm{\σ}}_l^2\right)\}+N{p}_m\{-\frac{1}{2}[1+\ln \left(2\mathrm{\π}\right)]-\frac{1}{2}\ln \left({\mathrm{\σ}}_m^2\right)\}$

$-2N{p}_{(l,m)}\{-\frac{1}{2}[1+\ln \left(2\mathrm{\π}\right)]-\frac{1}{2}\ln \left({\mathrm{\σ}}_{(l,m)}^2\right)\}$

$=\frac{N{p}_l}{2}\ln \left(\frac{{\mathrm{\σ}}_{(l,m)}^2}{{\mathrm{\σ}}_l^2}\right)+\frac{N{p}_m}{2}\ln \left(\frac{{\mathrm{\σ}}_{(l,m)}^2}{{\mathrm{\σ}}_m^2}\right)---\left(18\right)$

将满足 $(l^{*},m^{*})=\arg \underset{(l,m)}{\min }d(l,m)$ 的l、m合并。

所述迭代终止条件为：

|[lnp(f|k,Θ^(t+1))-lnp(f|k,Θ^(t))]|/|lnp(f|k,Θ^(t))|<0.00001；当满足此条件时，t即为确定的迭代次数。

在所述步骤1与步骤2之间还包括数据预处理步骤：

对所选的2幅多光谱遥感影像进行辐射校正及几何校正预处理；

其中，所述数据预处理具体实现为：首先进行辐射校正，然后进行几何粗校正，最后进行几何精校正；

步骤A，对于辐射校正，将目标影像数据采用基于统计量的方法进行相对辐射校正，在假定目标数据的各个波段服从高斯分布的情况下，其表达式为：

$\frac{g_d-{\mathrm{\&mu;}}_f}{{\mathrm{\&sigma;}}_f}=\frac{g_r-{\mathrm{\&mu;}}_r}{{\mathrm{\&sigma;}}_r}---\left(19\right)$

式中，g_d为待校正影像校正后的像素灰度值，μ_f，σ_f分别为参考影像的像素均值和标准差，g_r，μ_r，σ_r分别为待校正影像的像素灰度值、均值和标准差，则(19)式可变形得到：

$g_d=\frac{g_r-{\mathrm{\&mu;}}_r}{{\mathrm{\&sigma;}}_r}\&times;{\mathrm{\&sigma;}}_f+{\mathrm{\&mu;}}_f---\left(20\right)$

步骤B，对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：1.显示基准影像和待校正影像；2.采集地面控制点；3.计算误差；4.选择多项式模型；5.采用双线性插值进行重采样输出；

其中，双线性差值法为，若求未知函数f在点P=(x,y)的值，假设已知函数f在Q₁₁=(x₁,y₁),Q₁₂=(x₁,y₂),Q₂₁=(x₂,y₁),及Q₂₂=(x₂,y₂)四个点的值；如果选择一个坐标***使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式表示为：

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy    (21)

步骤C，对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正；

其中三角剖分法为，采用逐点***法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。

本发明与现有的技术相比具有以下优点：

1.本发明通过相关系数融合的方法来构造差异影像，可以克服常用的代数法仅依靠单一方法获取差异影像，不能较好地反映复杂地物的真实变化情况的问题。

2.本发明通过引入PCA变换，以实现最大限度地消除各光谱波段间的相关性，从而有效地集中和突出不同时相间的变化信息，同时又能最大限度地保留变化的细节信息，减少伪变化信息及背景噪声的干扰，最大限度提高差异影像所包含的变化信息，可以克服采用单一波段进行融合，所造成的包含变化信息量不足的问题。

3.本发明通过MDL-EM算法估计GMM的各参数，从而实现对差异影像直方图的最佳拟合。有效克服了当遥感地物类型较复杂时，通过EM算法求得的GMM对差异影像直方图拟合效果差，检测精度低的不足。

附图说明

图1是本发明基于一种多光谱遥感影像变化检测方法的实现流程示意图。

具体实施方法

下面参照图1，对本发明的具体实施步骤做进一步详细描述：

步骤1数据准备

选取需要进行变化检测的同一区域、不同时相的2幅多光谱遥感影像数据。

步骤1b数据预处理

对所选的2幅多光谱遥感影像进行辐射校正及几何校正等预处理。

在预处理过程中，首先进行辐射校正，然后进行几何粗校正，最后进行几何精校正。

对于辐射校正，将目标影像数据采用基于统计量的方法进行相对辐射校正，在假定目标数据的各个波段服从高斯分布的情况下，其表达式为：

$\frac{g_d-{\mathrm{\&mu;}}_f}{{\mathrm{\&sigma;}}_f}=\frac{g_r-{\mathrm{\&mu;}}_r}{{\mathrm{\&sigma;}}_r}---\left(1\right)$

式中，g_d为待校正影像校正后的像素灰度值，μ_f，σ_f分别为参考影像的像素均值和标准差，g_r，μ_r，σ_r分别为待校正影像的像素灰度值、均值和标准差，则(1)式可变形得到：

$g_d=\frac{g_r-{\mathrm{\&mu;}}_r}{{\mathrm{\&sigma;}}_r}\&times;{\mathrm{\&sigma;}}_f+{\mathrm{\&mu;}}_f---\left(2\right)$

对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：（1）显示基准影像和待校正影像；（2）采集地面控制点；（3）计算误差；（4）选择多项式模型；（5）采用双线性插值进行重采样输出。

双线性差值法，若求未知函数f在点P=(x,y)的值，假设我们已知函数f在Q₁₁=(x₁,y₁),Q₁₂=(x₁,y₂),Q₂₁=(x₂,y₁),及Q₂₂=(x₂,y₂)四个点的值。如果选择一个坐标***使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy    (3)

对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

三角剖分法为，采用逐点***法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影像。

步骤2构造差异影像

将经过步骤2预处理的2幅多光谱遥感影像数据，通过相关系数融合法来构造差异影像。

差异影像的构造可分以下几步进行，首先对经过步骤2得到的2幅多光谱遥感影像数据X₁，X₂，通过ENVI4.8软件分别进行PCA变换，并提取其第一主成分X_1pc1，X_2pc1。然后，求取X_1pc1，X_2pc1的差值、比值影像ΔY₁与ΔY₂，最后对获取的差值、比值影像进行融合构造出差异影像。

对于PCA变换，采用ENVI4.8软件对目标数据进行PCA变换时以协方差矩阵的特征向量对应的矩阵作为变换矩阵。具体操作步骤为：（1）打开目标影像文件。（2）在主菜单中，选择Transforms→Principal Components→Forward PC Rotation→Compute new Statistics andRotate。在Principal Components Input File对话框中，选择影像文件。在Calculate using中选择“Covariance Matrix”（协方差矩阵）。并将提取的第一主成分X_1pc1，X_2pc1输出。

对于n维随机变量（X₁,X₂,…,X_n），若c_ij=Cov(X_i,X_j)=E{[X_i-E(X_i)][X_j-E(X_j)](i,j=1,2...n)都存在，则n维随机变量（X₁,X₂,…,X_n）的协方差矩阵为：

$C=\left[\begin{array}{cccc}{c}_{11}& c_{12}& ...& c_{1n}\\ c_{21}& c_{22}& ...& c_{2n}\\ ...& ...& ...& ...\\ c_{n1}& c_{n2}& ...& c_{\mathrm{nn}}\end{array}\right]---\left(4\right)$

通过相关系数融合法构造的差异影像F在坐标(p,q)处的像元值为：

F(p,q)=λ_pq[αΔY₁(p,q)+βΔY₂(p,q)]    (5)

式中，λ_pq=ΔY₂(p,q)/max(ΔY₂)；α=a|r|+b;α+β=1；r为ΔY₁与ΔY₂的相关系数；其表达式为：

$r=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}(\mathrm{\&Delta;}{Y}_1-\stackrel{\&OverBar;}{\mathrm{\&Delta;}{Y}_1})(\mathrm{\&Delta;}{Y}_2-\stackrel{\&OverBar;}{\mathrm{\&Delta;}{Y}_2})/\sqrt{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(\mathrm{\&Delta;}{Y}_1-\stackrel{\&OverBar;}{\mathrm{\&Delta;}{Y}_1})}^2\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{(\mathrm{\&Delta;}{Y}_2-\stackrel{\&OverBar;}{\mathrm{\&Delta;}{Y}_2})}^2}---\left(6\right)$

式中，N——GMM的像元总数，

与ΔY₂的均值。a，b——常数。

常数a，b作为权重α的调节因子，a，b∈[0,1]。通过调节a，b的值来调整差值影像和比值影像的权重α,β，取a，b∈{0.1,0.2,…,1},则a，b的值可由下式确定。

$\{a,b\}=\arg \max \left\{\frac{1}{N}\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{n}{\mathrm{\&Sigma;}}}{[(p,q)-\frac{1}{N}\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{n}{\mathrm{\&Sigma;}}}F(p,q)]}^2\right\}---\left(7\right)$

式中，m——F的行数，n——F的列数，且有N=m×n。

步骤3建立模型

首先，采用GMM对差异影像进行建模，然后，利用MDL-EM算法估计GMM的参数（分支数，权重，均值，方差）。

对于GMM的构造，具体实现为：设f_i（i=1,…,N）为差异影像F的第i个像元灰度值；F中存在k个待识别的类（分支），每个像元之间是相互独立的，且不考虑分支间的相关性，则有：

$p\left(f_i\right)=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}p(f_i,j)=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}p\left(f_i\right|j)p\left(j\right)---\left(8\right)$

f_i属于分支j的概率密度函数为：

$p\left(f_i\right|j)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_j}\exp [-\frac{{(f_i-{\mathrm{\&mu;}}_j)}^2}{2{\mathrm{\&sigma;}}_j^2}]---\left(9\right)$

式中：μ_j，

——分支j的均值和方差。

设：

Θ=(p(1),p(2),...,p(k);θ₁,θ₂,...,θ_k)。则对应的完全数据集的对数似然函数为：

$\ln p\left(f\right|k,\mathrm{\&Theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\ln \left[\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}p\left(f_i\right|j,\mathrm{\&theta;})p\left(j\right)\right]---\left(10\right)$

对于参数估计，采用MDL-EM算法来估计所建立的GMM的参数，具体步骤为：设L为差异影像灰度值集合（一般为0,1,…,255），h(f_i)为直方图，其中f_i为其第i个像元灰度值，则利用EM算法估计GMM的参数（权重，均值，方差）的迭代公式为：

$p^{t+1}\left(j\right)=\frac{1}{f_{i}h\left(f_i\right)}\underset{f_i\&Element;L}{\mathrm{\&Sigma;}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}---\left(11\right)$

${\mathrm{\&mu;}}_j^{t+1}=\frac{\underset{f_i\&Element;L}{\mathrm{\&Sigma;}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}{f}_i}{\underset{f_i\&Element;L}{\mathrm{\&Sigma;}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}}---\left(12\right)$

${\left({\mathrm{\&sigma;}}_j^2\right)}^{t+1}=\frac{\underset{f_i\&Element;L}{\mathrm{\&Sigma;}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}{(f_i-{\mathrm{\&mu;}}_j^t)}^2}{\underset{f_i\&Element;L}{\mathrm{\&Sigma;}}\frac{h\left(f_i\right)p^t\left(f_i\right|j)p^t\left(j\right)}{p^t\left(f_i\right)}}---\left(13\right)$

式中：p(j)，μ_j，

——GMM中分支j的权重、均值和方差；t——迭代次数。迭代终止条件为：|[lnp(f|k,Θ^(t+1))-lnp(f|k,Θ^(t))]|/lnp(f|k,Θ^(t))|<0.00001，当满足此条件时，t即为确定的迭代次数。而GMM最优模型的选择则由作为适应度函数的MDL准则来确定，其表达式为：

$R(k,\mathrm{\&Theta;})=-\ln P\left(f\right|k,\mathrm{\&Theta;})+\frac{T\ln \mathrm{NM}}{2}---\left(14\right)$

其中，

k——GMM的分支数；M——GMM中像元的维数；

lnP(f|k,Θ)——差异影像对应的完全数据集的对数似然函数；TlnNM/2——引入的针对GMM分支数k的惩罚函数。

在选定初始最大分支数k_max的情况下，通过引入合并最小距离分支的办法来自适应地依次减少分支数直到k_min，并依据下式，确定最优的分支数k^*和参数Θ^*。

{k^*,Θ^*}=argmin{R(k,Θ),k_min≤k≤k_max}   (15)

合并最小距离法，上述合并最小距离分支法，具体实现为：设分支数为k时其任意两类分支l，m(1≤l,m≤k)，的参数为：p_l,p_m,μ_l,μ_m,

合并后的新分支参数为：p_(l,m)，μ_(l,m)，

首先令μ_l=μ_m=μ_(l,m)；p_l，p_m不变，则在分支数为k时有θ_(l,m)=θ_(l)=θ_(m)∈Θ^(k)，同时注意到因为θ_(l)=θ_(m)所以可将他们视为一个分支，为了便于区分，将Θ^(k-1)中的θ_(l,m)表示成

则有

此时有p_(lm)=p_l+p_m；再利用θ_(lm)和

的定义，结合MDL准则有：

$R(k-1,{\mathrm{\&theta;}}_{{(l,m)}^{-}})-R(k,{\mathrm{\&theta;}}^{\left(t\right)})$

$=R(k-1,{\mathrm{\&theta;}}_{{(l,m)}^{-}})-R(k,{\mathrm{\&theta;}}_{(l,m)})+R(k,{\mathrm{\&theta;}}_{(l,m)})-R(k,{\mathrm{\&theta;}}^{\left(t\right)})$

$\&le;\frac{1}{2}(1+M+\frac{M(M+1)}{2})\ln N+Q({\mathrm{\&theta;}}^{\left(t\right)};{\mathrm{\&theta;}}^{\left(t\right)})-Q({\mathrm{\&theta;}}_{(l,m)};{\mathrm{\&theta;}}^{\left(t\right)})---\left(16\right)$

$\&le;\frac{1}{2}(1+M+\frac{M(M+1)}{2})\ln N$

$+Q({\mathrm{\&theta;}}^{\left(t\right)};{\mathrm{\&theta;}}^{\left(t\right)})-Q({\mathrm{\&theta;}}^{*};{\mathrm{\&theta;}}^{\left(t\right)})+Q({\mathrm{\&theta;}}^{*};{\mathrm{\&theta;}}^{\left(t\right)})-Q({\mathrm{\&theta;}}_{(l,m)}^{*};{\mathrm{\&theta;}}^{\left(t\right)})$

在这里θ^*和为通过EM算法求得的最优参数估计，Q(θ;θ^(t))为根据上一次计算得的结果来估计完全数据集Z的似然函数的期望值，其表达式为：

Q(θ;θ^(t))=E_z[lnp(f,Z|θ)|f,θ^(t)]    (17)

此时有θ^*=θ^(t)，则Q(θ^(t);θ^(t))-Q(θ^*;θ^(t))=0。由EM算法的M步可知：

则对于任意的l，m其合并分支(l,m)的参数经EM算法迭代求得的最优估计为：

p_(l,m)=p_l+p_m    (18)

${\mathrm{\&mu;}}_{(l,m)}=\frac{p_l{\mathrm{\&mu;}}_l+p_m{\mathrm{\&mu;}}_m}{p_l+p_m}---\left(19\right)$

${\mathrm{\&sigma;}}_{(l,m)}^2=\frac{p_l[{\mathrm{\&sigma;}}_l^2+{({\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_{(l,m)})}^2]+p_m[{\mathrm{\&sigma;}}_m^2+{({\mathrm{\&mu;}}_m-{\mathrm{\&mu;}}_{(l,m)})}^2]}{p_l+p_m}---\left(20\right)$

利用上式便可将k中任意两个分支l，m合并从而将k降阶到k-1。再定义一个距离函d(l,m)：

$d(l,m)=N{p}_l\{-\frac{1}{2}[1+\ln \left(2\mathrm{\&pi;}\right)]-\frac{1}{2}\ln \left({\mathrm{\&sigma;}}_l^2\right)\}+N{p}_m\{-\frac{1}{2}[1+\ln \left(2\mathrm{\&pi;}\right)]-\frac{1}{2}\ln \left({\mathrm{\&sigma;}}_m^2\right)\}$

$-2N{p}_{(l,m)}\{-\frac{1}{2}[1+\ln \left(2\mathrm{\&pi;}\right)]-\frac{1}{2}\ln \left({\mathrm{\&sigma;}}_{(l,m)}^2\right)\}$

$=\frac{N{p}_l}{2}\ln \left(\frac{{\mathrm{\&sigma;}}_{(l,m)}^2}{{\mathrm{\&sigma;}}_l^2}\right)+\frac{N{p}_m}{2}\ln \left(\frac{{\mathrm{\&sigma;}}_{(l,m)}^2}{{\mathrm{\&sigma;}}_m^2}\right)---\left(21\right)$

将满足 $(l^{*},m^{*})=\arg \underset{(l,m)}{\min }d(l,m)$ 的l，m合并。

步骤5确定变化检测阈值

利用基于统计最小错误率的Bayes判别准则确定变化检测的阈值从而实现对多光谱遥感影像的变化检测。具体实现为：令与类j对应的判别函数为g_j(f_i)=p(f_i|j)p(j)，如果g_j(f_i)=max{g₁(f_i),g₂(f_i),...,g_k(f_i)}则把f_i判为j。而在各分支交叠不是很严重的情况下,可把各分支按均值μ_j从小到大排列其相邻两分支j和j+1之间的阈值可以由其交点确定。

确定了阈值后，再对步骤3得到的差异影像进行二值化处理便能得到目标数据的变化检测结果。

Claims (5)

1.一种多光谱遥感影像变化检测方法，其特征在于：包括下述步骤： 

步骤1，数据准备： 

选取需要进行变化检测的同一区域、不同时相的2幅多光谱遥感影像数据； 

步骤2，构造差异影像： 

对2幅多光谱遥感影像数据通过相关系数融合法构造差异影像； 

差异影像的构造分以下几步进行，首先对经过步骤1得到的2幅多光谱遥感影像数据X₁，X₂，通过ENVI4.8软件分别进行PCA变换，并分别提取其第一主成分X_1pc1，X_2pc1；然后，求取X_1pc1，X_2pc1的差值就影像ΔY₁以及比值影像ΔY₂，最后对获取的差值、比值影像进行融合构造出差异影像； 

对于PCA变换，采用ENVI4.8软件对目标数据进行PCA变换时以协方差矩阵的特征向量对应的矩阵作为变换矩阵； 

对于n维随机变量X₁,X₂,…,X_n，若c_ij=Cov(X_i,X_j)=E{[X_i-E(X_i)][X_j-E(X_j)]都存在，i,j=1,2,…n,则n维随机变量X₁,X₂,…,X_n的协方差矩阵为： 

通过相关系数融合法构造的差异影像F在坐标(p,q)处的像元值为： 

F(p,q)=λ_pq[αΔY₁(p,q)+βΔY₂(p,q)]   (2) 

式中，λ_pq=ΔY₂(p,q)/max(ΔY₂)；α=a|r|+b;α+β=1；r为ΔY₁与ΔY₂的相关系数；其表达式为： 

式中，N为GMM的像元总数，

分别为ΔY₁与ΔY₂的均值，a，b为常数； 

常数a，b作为权重α的调节因子，a，b∈[0,1]；通过调节a，b的值来调整差值影像和比值影像的权重α,β，取a，b∈{0.1,0.2,…,1},则a，b的值可由下式确定； 

式中，m为差异影像F的行数，n为差异影像F的列数，且有N=m×n； 

步骤3，建立模型： 

首先，采用GMM对差异影像进行建模，然后，利用EM算法估计GMM的参数，所述参数包括权重，均值，方差； 

步骤4，确定变化检测阈值： 

利用基于统计最小错误率的Bayes判别准则确定变化检测的阈值，得到多光谱遥感影像的变化检测结果。

2.根据权利要求1所述的一种多光谱遥感影像变化检测方法，其特征在于：所述步骤3包括如下步骤： 

首先，采用GMM对差异影像进行建模，然后，利用MDL-EM算法估计GMM的参数； 

对于GMM的构造，具体实现为：设fi为差异影像F的第i个像元灰度值；i=1,…,N，F中存在k个待识别的类，即为k个分支；每个像元之间是相互独立的，且不考虑分支间的相关性，则有： 

f_i属于分支j的概率密度函数为： 

式中：p(j)为GMM中分支j的权重，μ_j，

分别代表分支j的均值和方差；j＝1，2…k; 

设：

Θ=(p(1),p(2),...,p(k);θ₁,θ₂,...,θ_k)，则对应的完全数据集的对数似然函数为 

对于参数估计，采用MDL-EM算法来估计所建立的GMM的参数，具体步骤为：设L为差异影像灰度值集合，其一般为0,1,…,255，h(f_i)为直方图，其中f_i为其第i个像元灰度值，则利用EM算法估计GMM的参数权重，均值，方差的迭代公式为： 

t为迭代次数；而GMM最优模型的选择则由作为适应度函数的MDL准则来确定，其表达式为： 

其中，

k代表GMM的分支数；M代表GMM中像元的维数； 

lnP(f|k,Θ)代表差异影像对应的完全数据集的对数似然函数；TlnNM/2代表引入的针对GMM分支数k的惩罚函数。 

3.根据权利要求2所述的一种多光谱遥感影像变化检测方法，其特征在于：所述步骤3中的GMM中分支数k通过如下方法确定： 

在选定初始最大分支数k_max的情况下，通过引入合并最小距离分支的办法来自适应地依次减少分支数直到k_min，并依据下式，确定最优的分支数k^*和参数Θ^*； 

{k^*,Θ^*}=argmin{R(k,Θ),k_min≤k≤k_max}   (12) 

利用上述MDL-EM算法便能自适应地估计得到GMM的k^*和Θ^*； 

上述合并最小距离分支法，具体实现为：设分支数为k时其任意两类分支l、m，其中1≤l，m≤k，其参数为：p_l,p_m,μ_l,μ_m,

合并后的新分支参数为：p(_l,m)，μ_(l,m)，

首先令μ_l=μ_m=μ_(l,m)；p_l，p_m不变，则在分支数为k时有θ_(l,m)=θ_(l)=θ_(m)∈Θ^(k)，同时注意到因为θ_(l)=θ_(m)所以可将他们视为一个分支，为了便于区分，将Θ^(k-1)中的θ_(l,m)表示成

则有

此时有p_(l,m)=p_l+p_m；再利用θ_(l,m)和

的 定义，结合MDL准则有： 

R(k-1,θ_(l,m)-)-R(k,θ^(t)) 

其中θ^*和

为通过EM算法求得的最优参数估计，Q(θ;θ^(t))为根据上一次计算得的结果来估计完全数据集Z的似然函数的期望值，其表达式为： 

此时有θ^*=θ^(t)，则Q(θ^(t);θ^(t))-Q(θ^*;θ^(t))=0；由EM算法的M步可知： 

则对于任意的l、m其合并分支(l,m)的参数经EM算法迭代求得的最优估计为： 

p_(l,m)=p_l+p_m(15) 

利用上式便可将k中任意两个分支l、m合并从而将k降阶到k-1；再定义一个距离函d(l,m)： 

将满足

的l、m合并。 

4.根据权利要求3所述的一种多光谱遥感影像变化检测方法，其特征在于：迭代终止条件为：|[lnp(f|k,Θ^(t+1))-lnp(f|k,Θ(t))]|/lnp(f|k,Θ^(t))|<0.00001；当满足此条件时， t即为确定的迭代次数。 

5.根据权利要求1所述的一种多光谱遥感影像变化检测方法，其特征在于：在所述步骤1与步骤2之间还包括数据预处理步骤： 

对所选的2幅多光谱遥感影像进行辐射校正及几何校正预处理； 

其中，所述数据预处理具体实现为：首先进行辐射校正，然后进行几何粗校正，最后进行几何精校正； 

步骤A，对于辐射校正，将目标影像数据采用基于统计量的方法进行相对辐射校正，在假定目标数据的各个波段服从高斯分布的情况下，其表达式为： 

式中，g_d为待校正影像校正后的像素灰度值，μ_f，σ_f分别为参考影像的像素均值和标准差，g_r，μ_r，σ_r分别为待校正影像的像素灰度值、均值和标准差，则(19)式可变形得到： 

步骤B，对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：1.显示基准影像和待校正影像；2.采集地面控制点；3.计算误差；4.选择多项式模型；5.采用双线性插值进行重采样输出； 

其中，双线性差值法为，若求未知函数f在点P=(x,y)的值，假设已知函数f在Q₁₁=(x₁,y₁),Q₁₂=(x₁,y₂),Q₂₁=(x₂,y₁),及Q₂₂=(x₂,y₂)四个点的值；如果选择一个坐标***使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式表示为： 

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy    (21) 

步骤C，对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正； 

其中三角剖分法为，采用逐点***法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。 

Images