[发明专利]基于treelet融合和水平集分割的遥感图像变化检测

摘要

本发明公开了一种基于treelet融合和水平集分割的遥感图像变化检测方法，主要解决现有变化检测方法存在较多伪变化信息的问题。其实现过程是：输入两时相遥感图像，对每幅图像分别进行均值漂移滤波，得到两时相滤波后图像并分别对其进行3次不同层数下的二维平稳小波分解，对相同分解层数对应方向子带的小波系数矩阵做差；采用sobel算子对得到的水平、垂直方向小波系数差矩阵进行增强并进行二维小波逆变换重构；采用treelet算法融合不同分解层数的重构图像得到最终的差异图，对该差异图进行水平集分割得到变化检测结果。本发明能够有效提高变化检测结果的精度，同时较好的保持变化区域的边缘特征，可用于对自然灾害的分析、土地资源监测等领域。

主权项

1.一种基于treelet融合和水平集分割的遥感图像变化检测方法，包括如下步骤：（1）对输入的两幅已配准的大小均为mxn的多时相遥感图像I¹和I²分别进行均值漂移滤波，得到滤波后图像X¹和X²；（2）对滤波后图像X¹和X²分别进行3次二维平稳小波分解，且每一幅滤波后图像3次分解的最高分解层数分别为s=1，...3，一次二维平稳小波分解得到四个小波系数矩阵，即一个低频系数矩阵和三个分别表示水平、垂直、对角方向的高频小波系数矩阵；（3）对滤波后图像X¹和X²相同分解层数下对应方向子带的小波系数矩阵做差，得到每一个分解层s的低频小波系数差矩阵和三个分别表示水平、垂直、对角方向的高频小波系数差矩阵；（4）对步骤（3）中的水平方向小波系数差矩阵和垂直方向小波系数差矩阵利用sobel算子进行增强，保持低频小波系数差矩阵和对角方向小波系数差矩阵不变；（5）使用（3）中低频小波系数差矩阵、对角方向小波系数差矩阵和（4）中增强后的水平、垂直方向小波系数差矩阵，对每一个分解层s进行二维逆平稳小波变换，得到每个分解层s的重构图像RIs；（6）对每个分解层s的重构图像RIs使用treelet变换进行融合，得到融合后的差异图D；（7）对融合后的差异图D进行水平集分割，得到变化检测结果图Z；其中步骤（6）所述的采用treelet算法对重构后的图像RIs进行融合，按如下步骤进行：（6a）将重构后的图像RIs均变换为mxn大小的列向量RI1′，RI2′，RI3′，组成初始样本X=[RI1’,RI2’,RI3′]；（6b）定义treelet变换的逐层聚类层数l＝0,1，...L-1，L为初始样本X中列向量的个数，L=3，在第0层，每个变量采用初始样本X中的列向量表示，初始化基矩阵B⁰为LxL的单位矩阵及和变量下标集合Ω={1,2，...L}，计算样本X的初始协方差矩阵C⁰和初始相关系数矩阵M⁰，计算公式如下：C_ij＝E[(s_p-Es_p)(s_q-Es_q)^T]Mij=CijCiiCjj]]>式中C_ij表示初始协方差矩阵C⁰第i行第j列的计算值，s_p，s_q表示样本X中两个不同的列向量，E表示求数学期望，T表示转置操作，M_ij表示初始相关系数矩阵M⁰第i行第j列的计算值；（6c）当逐层聚类层数l≠0时，寻找相关系数矩阵M^l中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：(α,β)=argmaxi,j∈ΩMijl-1]]>这里i＜j，分别表示相关系数矩阵M^l中任意值的行和列,且只在和变量下标集合Ω内进行；（6d）计算Jacobi旋转矩阵J：其中c_n=cos(θ₁)，s_n=sin(θ₁)；旋转角θ_l由C^l＝J^TC^l-1J及计算得到，J^T表示旋转矩阵J的转置：θl=tan-1[Cααl-Cββl±(Cαβl-Cββl)2+4(Cαβl)22Cαβl]]]>且|θl|≤π4;]]>式中，表示协方差矩阵C^l中行、列位置序号均为α的元素值，表示协方差矩阵C^l中行、列位置序号均为β的元素值，表示协方差矩阵C^l中行位置序号为α、列位置序号均为β的元素值；（6e）由矩阵J更新去相关后的正交基B^l＝B^l-1J，协方差矩阵C^l＝J^TC^l-1J，其中B^l-1为更新前的正交基，B^l为更新后的正交基，C^l-1为更新前的协方差矩阵，C^l为更新后的协方差矩阵，并利用C^l更新相关系数矩阵M^l-1为M^l，M^l-1为更新前的相关系数矩阵，M^l为更新后的相关系数矩阵：Mijl=CijlCiilCjjl]]>其中为更新后相关系数矩阵M^l中第i行j列的更新值，为更新后协方差矩阵C^l中第i行j列的值，为更新后协方差矩阵C^l中第i行i列的值，为更新后协方差矩阵C^l中第j行j列的值；（6f）定义更新后正交基矩阵B^l中位置序号为α和位置序号为β的两个列向量分别为尺度函数φ_l和细节函数ψ_l，定义当前l层的尺度向量集合{φ_l}是尺度函数φ_l和上一层的尺度向量集合{φ_l-1}的合集，将差变量的位置序号β从和变量的下标集合Ω中去除，即Ω=Ω\{β}；（6g）重复步骤（6c）至（6f）直至聚类的最高层l＝L-1，可得最终的正交基矩阵B=[φ_L-1,ψ₁，...ψ_L-1]，其中φ_L-1∈{ψ_L-1}，ψ₁为第一层treelet变换得到的细节函数，ψ_L-1为最高层treelet变换得到的细节函数；（6h）将初始样本X沿正交基矩阵B转置的方向进行投影，即P=XxB^T，并将投影向量P变回mxn大小的图像，得到融合后的差异图D。