[发明专利]基于Treelet图像融合的遥感图像变化检测方法

权利要求书

1.一种基于Treelet图像融合的遥感图像变化检测方法，包括如下步骤：

（1）输入变化前的图像P₁和变化后的图像P₂，分别进行中值滤波和均值滤波，得到中值滤波后的变化前图像P_1m，均值滤波后的变化前图像P_1a，中值滤波后的变化后图像P_2m，均值滤波后的变化后图像P_2a；

（2）用差值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₁，用对数比值法得到中值滤波后的变化前图像P_1m和中值滤波后的变化后图像P_2m的差异图H₂，用差值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₃，用对数比值法得到均值滤波后的变化前图像P_1a和均值滤波后的变化后图像P_2a的差异图H₄；

（3）分别将上述差异图H₁，H₂，H₃，H₄转化成一维的列向量，X₁，X₂，X₃，X₄，并将该4个列向量组成一个矩阵X=[X₁X₂X₃X₄]；

（4）对矩阵X用treelet法进行降维，得到列向量P，并将该向量P转化成最终差异图D的灰度矩阵；

（5）对所得的最终差异图D的灰度矩阵进行聚类分割，得到变化检测结果并输出。

2.根据权利要求书1中所述的方法，其中步骤(4)所述的对矩阵X用treelet法进行降维处理，降为一维并将该向量转化成最终差异图D的灰度矩阵，按如下步骤进行：

4.1）计算X的协方差矩阵C：

C=C1,1C1,2...C1,4C2,1C2,2...C2,4............C4,1C4,2...C4,4]]>

其中协方差矩阵C中的每一个元素C_m，n按照如下公式计算：

Cm,n=Σk=1N(xk,m-x‾m)(xk,n-x‾n)Σk=1N(xk,m-x‾m)2Σk=1N(xk,n-x‾n)2]]>

和分别表示矩阵X第m列和第n列元素的均值,x_k，m,为矩阵X第k行第m列的值，x_k，n为矩阵X第k行第n列的值，m,n=1,2,3,4，k＝1,2,3,…N，N为差异图像素点数；

4.2）计算相关系数矩阵ρ，以此来度量两个变量之间的相似性

ρ=ρ1,1ρ1,2...ρ1,4ρ2,1ρ2,2...ρ2,4............ρ4,1ρ4,2...ρ4,4]]>

其中

ρa,b=Ca,bCa,bCa,b]]>

Ρ_a，b表示矩阵ρ的第a行、第b列个元素，C_a，b为协方差矩阵C的第a行、第b列个元素,a，b=1,2,3,4；

4.3）初始化总的聚类层数为l=0，初始化正交Dirac基B⁽⁰⁾为4x4维的单位矩阵，C⁽⁰⁾=C；

4.4)设l＝l+1，寻找相关系数矩阵ρ中最大的值，将该值的座标序号分别标记为α和β，α>β；

4.5)利用得到的位置序号α和β，计算Jacobi旋转矩阵J：

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

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

为第l-1层协方差矩阵C在α行α列的值，为第l-1层协方差矩阵C在α行β列的值，为第l-1层协方差矩阵C在β行β列的值；

4.6)利用旋转矩阵J得到去相关后的正交基B^(l)=B^(l-1)J，并更新协方差矩阵C^(l)=J^TC^(l-1)J，其中B^(l)为更新后的正交基，C^(l)为更新后的协方差矩阵，并根据协方差矩阵C^(l)更新相关系数矩阵ρ^(l);

4.7）判断l的值，如果l＝3执行步骤4.8），否则返回步骤4.4）；

4.8）将矩阵X沿正交基矩阵B⁽³⁾的转置方向进行投影，得到投影后的矩阵P_n：

P_n=X×B^(3)T；

4.9）将投影后的矩阵P_n的第一列P转化为与输入图像同等大小的矩阵，该矩阵即为融合后最终差异图D的灰度矩阵。

3.根据权利要求书1中所述的方法，其中步骤(5)所述的对所得的最终差异图D的灰度矩阵进行聚类分割，得到变化检测结果，按如下步骤进行：

5.1)将最终差异图D的灰度矩阵转化为一个列向量T=[t_i]，t_i为该向量第i个像素的灰度值，i=1,2,3,…N，N为总像素点数，将第i个像素属于第k类的隶属度记为u_ki，u_ki∈[0,1]且k＝1,2，随机初始化隶属度矩阵U^b=[u_ki]的每一个元素，设定循环计数器b=0；

5.2)按如下公式计算第k类的聚类中心c_k：

ck=Σi=1Nuki2tiΣi=1Nuki2;]]>

5.3)按如下公式计算隶属度矩阵U^（b+1）:

U^（b+1)={u_ki}，

其中，k＝1,2,u1i=11+(ti-c1ti-c2)2,]]>u2i=1(ti-c2ti-c1)2+1;]]>

5.4)如果满足max{U^(b)-U^(b+1)}<0.001则转至步骤5.5），否则设置b=b+1,转至步骤5.2)；

5.5)通过比较c₁与c₂的大小确定变化类与非变化类，如果c₁≥c₂，则以c₁为聚类中心的那一类为变化类，以c₂为聚类中心的那一类为非变化类；如果c₁<c₂，以c₂为聚类中心的那一类为变化类，以c₁为聚类中心的一类为非变化类。