[发明专利]基于多组滤波和降维的遥感图像变化检测方法

权利要求书

1.一种基于多组滤波和降维的遥感图像变化检测方法，包括如下步骤：

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

(2)计算差值差异图Xd变化类和非变化类的分类阈值；

(3)根据分类阈值计算出图像T1和图像T2中非变化类像素的灰度直方图，将图像T1中的非变化类像素和图像T2中的非变化类像素进行灰度校正，得到图像T2灰度校正后的图像T3；将图像T1和图像T3空间位置对应的像素灰度值进行差值计算，得到一幅灰度校正的差值差异图D；

(4)将Gabor滤波器函数参数中5个方向角度的任意一个与6个频率的任意一个进行两两组合，可得到30组滤波器函数，利用该组滤波器函数对灰度校正的差值差异图D进行滤波得到30组滤波后的差异图像组，记为滤波图像组矩阵E₁；

(5)将灰度校正的差值差异图D加入到滤波图像组矩阵E₁中，得到降维样本矩阵E；

(6)利用Treelet算法计算降维样本矩阵E的投影矩阵P，得到沿投影方向降维后的最终差异图像F；

(7)对得到的最终差异图像F利用K均值算法进行聚类，得到最终的变化检测结果图Z。

2.根据权利要求1所述的遥感图像变化检测方法，其中步骤(2)所述计算差值差异图Xd变化类和非变化类的分类阈值，按如下步骤计算：

(2a)初始化分类阈值T为差值差异图Xd中最大像素值和最小像素值的均值；

(2b)利用初始化分类阈值T将差值差异图Xd分成两类，分别计算变化类像素值的均值m_c和非变化类像素值的均值m_d，将m_c和m_d的均值记为T_n，当初始化分类阈值T与T_n的绝对差值大于收敛值T_s时，更新分类阈值T为T_n，T_s＜＜1；

(2c)重复步骤(2b)直至初始化分类阈值T与均值T_n的绝对差值满足收敛值T_s时，得到最终的分类阈值T_m。

3.根据权利要求1所述的遥感图像变化检测方法，其中步骤(3)所述将图像T1中的非变化类像素和图像T2中的非变化类像素进行灰度校正，按如下步骤计算：

(3a)根据差值差异图Xd的分类阈值T_m将两幅已经配准的多时相遥感图像T1和T2分成两类，得到两幅图像各自的变化类和非变化类，计算出该图像T1中的非变化类像素的灰度直方图p₁和该图像T2中非变化类像素的灰度直方图p₂；

(3b)以灰度直方图p₁为参考，通过直方图匹配使灰度直方图p₂尽可能地接近灰度直方图p₁，得到图像T2灰度校正后的图像T3；

(3c)将多时相遥感图像T1和校正后的图像T3空间位置对应的像素灰度值进行差值计算，得到灰度校正的差值差异图D。

4.根据权利要求1所述的遥感图像变化检测方法，其中步骤(5)所述利用Treelet算法计算降维样本矩阵E的投影矩阵P，得到沿投影方向降维后的最终差异图像F的过程，按如下步骤进行：

(6a)计算降维样本矩阵E的协方差系数矩阵C：

C=C1,1C1,2...C1,LC2,1C2,2...C2,L............CL,1CL,2...CL,L]]>

其中，协方差系数矩阵C中的元素C_i，j的计算公式如下

Ci,j=Σk=1M×N(vk,i-vi‾)(vk,j-vj‾)Σk=1M×N(vk,i-vi‾)2Σk=1M×N(vk,j-vj‾)2]]>

其中C_i，j表示协方差矩阵C第i行第j列的计算值，和分别表示第i个和第j个滤波后图像列向量的均值，L为降维样本矩阵E的维数，此处L＝31，且满足1≤i，j≤L；

(6b)计算滤波后图像组矩阵E的相关系数矩阵A：

A=A1,1A1,2...A1,LA2,1A2,2...A2,L............AL,1AL,2...AL,L]]>

其中，相关系数矩阵A中的元素A_i，j的计算公式如下：

Ai,j=Ci,jCi,iCj,j]]>

其中A_i，j为相关系数矩阵A中第i行第j列的计算值，C_i，j为协方差矩阵C中第i行第j列的值，C_i，i为协方差矩阵C中第i行第i列的值，C_j，j为协方差矩阵C中第j行第j列的值；

(6c)定义分解的层数l＝0，1，…，L-1，当l＝0时，初始化和变量为滤波后图像组矩阵E，即S⁽⁰⁾＝E，差变量为D⁽⁰⁾为空集，并将和变量的集合Ω可用下标集表示为Ω＝{1，2，…，L-1}，将差变量集合Ф作为空集，将正交Dirac基作为L×L维的单位矩阵；

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

(α,β)=argmaxi,j∈ΩAi,j(l-1)]]>

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

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

其中c_n＝cos(θ_l)，s_n＝sin(θ_l)；旋转角θ_l由C^(l)＝J^TC^(l-1)J及计算得：

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

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

Ai,j(l)=Ci,j(l)Ci,i(l)Cj,j(l);]]>

(6f)定义正交基矩阵B^(l)中第α列和第β列的列向量分别为尺度函数和细节函数ψ_l，定义当前l层的尺度向量集合是尺度函数和上一层的尺度向量集合的合集，将β差变量从和变量集合Ω中去除，即Ω＝Ω/{β}，加入到差变量集合Ф中，即Ф＝{β}；

(6g)重复步骤(6d)至(6f)直至分解到最高层l＝L-1时，得到投影矩阵其中

(6h)按如下公式将滤波后图像组矩阵E沿投影矩阵P的方向降维，即

E_n＝E×P＝[E_sE_d]

其中E_n表示降维后的图像组矩阵，E_s为沿主要能量方向投影的滤波后图像矩阵，E_d表示沿次要能量方向投影的滤波后图像组矩阵，将E_s列向量变为与原始差值差异图X相同的大小，即得到降维后的最终差异图像F。