[发明专利]基于自适应差异图的遥感图像变化检测方法

权利要求书

1.一种基于自适应差异图的遥感图像变化检测方法，包括如下步骤：

(1)输入两幅已配准的多时相遥感图像T1和T2，并将图像T1和T2空间位置对应像素的灰度值进行差值计算，得到差值差异图X；

(2)对第一幅图像T1以像素点i为中心取3×3大小的滑窗E_i，对第二幅图像T2同样以像素点i为中心取7×7像素大小搜索窗G_i，i为当前像素点，在搜索窗G_i内逐像素取3×3大小的滑窗，得到滑窗集合{F_j}，F_j为滑窗集合内以j为中心点的一个滑窗，定义滑窗集合{F_j}内的滑窗数目为Q个，此处Q＝25；

(3)分别将第二幅图像T2得到的Q个滑窗与第一幅图像T1的滑窗Ei进行差值计算，并将差值的结果用列向量进行表示，得到样本矩阵H，该样本矩阵H中包含Q个列向量且维数为9×Q，将第二幅图像T2的滑窗F_i与第一幅图像T1的滑窗E_i的差值列向量表示为h_ii；

(4)利用Treelet算法对样本矩阵H进行自适应聚类，得到与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R，并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值；

(5)对两幅遥感图像T1和T2中的每个像素点，重复进行步骤2至步骤4，得到自适应差异图D；

(6)用Otsu阈值法分别计算差值差异图X和自适应差异图D的阈值K₁和K₂；

(7)利用Otsu阈值K₁和K₂融合差值差异图X和自适应差异图D，得到最终差异图I；

(8)对最终差异图I利用Otsu阈值法进行分割，得到变化检测的结果图像Z。

2.根据权利要求1所述的遥感图像变化检测方法，其中步骤(4)所述的利用Treelet算法对样本矩阵H进行自适应聚类，按如下步骤进行：

(4a)计算样本矩阵H的协方差系数矩阵C：

C=C1,1C1,2···C1,QC2,1C2,2···C2,Q·····Cu,v····CQ,1CQ,2···CQ,Q]]>

其中，协方差系数矩阵C中的任意一个元素C_uv的计算公式如下

Cu,v=Σk=19(xk,u-x‾u)(xk,v-x‾v)Σk=19(xk,u-x‾u)2Σk=19(xk,v-x‾v)2]]>

其中C_u，v表示协方差矩阵C中第u行第v列的计算值，和分别表示第u个和第v个差值列向量中像素灰度差值的均值，这里1≤u≤Q，，1≤v≤Q；

(4b)计算样本矩阵H的相关系数矩阵A：

A=A1,1A1,2···A1,QA2,1A2,2···A2,Q·····Au,v····AQ,1AQ,2···AQ,Q]]>

其中，相关系数矩阵A中的任意一个元素A_u，v的计算公式如下：

Au,v=Cu,vCu,uCv,v]]>

其中A_u，v为相关系数矩阵A中第u行第v列的计算值，C_u，v为协方差矩阵C中第u行第v列的值，C_u，u为协方差矩阵C中第u行第v列的值，C_v，v为协方差矩阵C中第v行第v列的值；

(4c)定义分解层数l＝0，1，…，Q-1，当l＝0时，初始化和变量为样本矩阵H，即S⁽⁰⁾＝H，差变量D⁽⁰⁾为空集，并将和变量的下标由集合Ω表示，Ω＝{1，2，…，Q}，将差变量下标集合Φ作为空集，将正交Dirac基B⁽⁰⁾＝[φ_0，1，φ_0，2，φ_0，3，…φ_0，L]作为Q×Q维的单位矩阵；

(4d)当分解层数l≠0时，寻找相关系数矩阵A中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

(α,β)=argmaxu,v∈ΩAu,v(l-1)]]>

此时仅考虑相关系数矩阵A的上三角部分，u和v分别表示相关系数矩阵A中任意值的行和列，且只在和变量集合Ω内进行；

(4e)计算Jacobi旋转矩阵J：

其中cn＝cos(θ_l)，sn＝sin(θ_l)；

旋转角θl通过C(l)＝JTC(l-1)J和这两个公式计算得到：

θl=tan-1[Cα,α-Cβ,β±(Cα,α-Cβ,β)2+4Cα,β22Cα,β]]]>且|θl|≤π4]]>

由矩阵J分别更新去相关后的正交基B^(l)＝B^(l-1)J和协方差矩阵C^(l)＝J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，并利用协方差矩阵C^(l)按下式更新相关系数矩阵A^(l)：

Au,v(l)=Cu,v(l)Cu,u(l)Cv,v(l);]]>

(4f)定义正交基矩阵B^(l)中第α列和第β列的列向量分别为尺度函数θ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将β差变量从和变量下标集合Ω中去除，即Ω＝Ω/{β}，加入到差变量下标集合Φ中，即Φ＝{β}，此时B(l)＝[φ_l，1，...φ_l，Q-l，ψ₁，...Ψ_l]，其中φ_l，1，...φ_l，Q-l∈{φ_l}；

(4g)计算样本矩阵H沿当前正交基矩阵B^(l)方向上的归一化能量矩阵Fn^(l)：

Fn(l)=[Fns(l)s∈Ω,Fnd(l)d∈Φ]=E||H·B(l)||2E||H||2]]>

其中E||·||²为计算数学期望，为归一化能量矩阵Fn^(l)中和变量下标集合Ω所对应向量组成的矩阵，为归一化能量矩阵Fn^(l)中差变量下标集合Φ所对应向量组成的矩阵；

计算和变量下标集合Ω所对应向量组成的矩阵的能量总和以及差变量下标集合Φ所对应向量组成的矩阵的能量总和

当满足Σs∈ΩFns(l)-Σs∈ΩFns(l-1)<Σd∈ΦFnd(l)]]>时，将当前分解层数l作为最终分解层数L；

(4h)利用正交基矩阵B^(L)对样本矩阵H进行投影，得到投影矩阵：

En＝(H×B^(L)′)·(H×B^(L)′)T，

提取投影矩阵En中的前Q-L行，得到尺度对应矩阵Em，将尺度对应矩阵Em中各列向量中最大值对应的行标号作为该列向量的类别标记；

(4i)寻找与差值列向量h_ii相同类别标记的q个列向量，其中1≤q≤Q，将这q个列向量组成矩阵R，并计算该矩阵R的均值m_i，将均值m_i作为自适应差异图D中像素点i的灰度值。