一种卫星影像变化检测方法与流程

1.本发明涉及多时相卫星影像变化检测领域，具体涉及一种卫星影像变化检测方法。


背景技术：

2.卫星影像变化检测通过对比分析多时相、同区域卫星影像，提取其中地物存在的变化信息，对环境调查、测绘地理空间数据更新、资源勘探、自然灾害救援与治理、战场信息动态感知、重点目标监测等领域有重要实际意义，是卫星影像处理的关键步骤之一。目前，常用的多时相卫星影像变化检测方法可分为三个层次：像素级、特征级及目标级，其中目标级变化检测主要检测具有一定概念意义的对象的变化，常用基于深度学习的方法，但检测结果依赖于目标识别和图像理解的准确性，难度较大。像素级变化检测基于原始图像上进行处理，变化检测前需进行校正、配准等预处理，现阶段的多时相卫星影像变化检测方法大多基于该层次，可大致分为基于代数运算的方法和基于变换域的方法。前者方法操作简单，耗时少，但易受光照、角度、云雾等干扰影响，无法准确描述地物目标的变化情况，变化检测精度不高。后者基于变换域的方法相较前者能大大减少数据冗余，且变化检测精度提升。其中使用最普遍的是基于主成分分析的方法，但该方法适用于信号的统计分布符合高斯分布的情况，针对卫星影像而言，其地物光谱特性往往不满足高斯分布，并不能完全消除各信号分量间的相关性，尤其是高阶相关信息，一定程度上影响了变化检测的准确性。


技术实现要素：

3.发明目的：针对现有船舶视频图像增强领域中存在的缺陷，进一步提高船舶视频图像增强的效果，本发明提供一种卫星影像变化检测方法，首先对某区域不同时期的两幅待检测卫星影像进行无下采样shearlet多尺度分解，真实保留其中的边缘及细节信息，抑制噪声等干扰，降低后续变化信息提取的复杂程度；然后通过基于负熵的独立性判决准则的方法改进传统独立主成分分析求解混矩阵，进一步提高收敛速度，对分解后的分量进行处理提取相应的变化信息，再经无下采样shearlet逆变换获得变化图像；最后基于改进倒数灰度熵阈值选取准则对变化图像进行分割，得到变化检测结果，进一步提高变化检测的准确性和运算效率。
4.技术方案：本发明公开了一种基于无下采样shearlet变换的卫星影像变化检测算法，包括如下步骤：
5.步骤1，对不同时期待变化检测的卫星遥感影像进行无下采样shearlet变换操作，得到相应的高频分量和低频分量；
6.步骤2，提取卫星遥感影像变化信息；
7.步骤3，将步骤1得到的高频分量和低频分量分别转化为矩阵形式，并进行无下采样shearlet变换逆变换得到多时相遥感变化图像；
8.步骤4，获得变化检测结果；
9.步骤5，进行阈值分割，得到最终的变化检测结果。
10.步骤1包括：
11.步骤1-1，合成膨胀仿射系统μ
ab
(ψ)表示为：
12.μ
ab
(ψ)＝{ψ
j,l,c
(x)＝|det(a)|
1/2
ψ(b
laj
x-c)}
ꢀꢀꢀꢀꢀꢀ
(1)
13.式中，为尺度参数，l2指平方可积，是实数集合、是整数集合；ψ
j,l,c
是基函数，x是一个变量，为剪切参数，尺度参数和剪切参数分别控制着无下采样shearlet变换的尺度分解数目和方向分解数目；c为平移参数；各向异性膨胀矩阵a、剪切矩阵b均为2
×
2的可逆矩阵，且|detb|＝1，det表示求行列式；f表示任意满足的函数变量，如果则称μ
ab
(ψ)中的元素为合成小波；b
l
是指剪切矩阵b的l次方；aj指各向异性膨胀矩阵a的j次方；
14.无下采样shearlet变换是合成小波的一种特殊情况，其中无下采样shearlet变换频域支撑区间近似看作大小2
2j
×2j
、斜率2-j
l的梯形对；
15.步骤1-2，多时相卫星影像多尺度分解；
16.步骤1-3，高频分量多方向剖分。
17.步骤1-2包括：采用无下采样金字塔对不同时期的两幅卫星遥感影像分别进行多尺度分解，选取尺度分解层数，得到低频分量和高频分量。
18.步骤1-3包括：构建带方向和尺度变化的meyer迈耶小波窗，分别对多尺度分解后的高频分量进行方向剖分及坐标映射，得到各方向子带的无下采样shearlet变换系数，选取方向剖分级数，即任意1个高频分量均有n1个不同的方向子带；
19.将多时相卫星遥感影像无下采样shearlet变换分解结果分别转化为向量形式，以低频分量和高频分量组成相应的两个分块向量。
20.步骤2包括：设一组观测信号e＝[e1,e2,
…
,en]
t
，它是n维源信号s＝[s1,s2,
…
,sn]
t
的观测值；第i个观测信号ei由n个独立源信号线性混合形成，即：
[0021]ei
＝a
i1
s1+a
i2
s2+
…
+a
in
sn,i＝1,2,
…
,m
ꢀꢀ
(2)
[0022]
其中m为混合矩阵的行数，a
in
为混合矩阵中第i行n列的值，sn为第n个源信号的观测值；
[0023]
用矢量表示公式(2)得到：
[0024]
e＝as
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0025]
式中，a＝[a1,a2,
…
,am]
t
为混合矩阵，am表示第m个混合矩阵基向量；ai是第i个混合矩阵基向量；ica就是仅通过观测数据e估计出混合矩阵a及它的逆矩阵w，最终估计出未知源信号y：
[0026]
y＝we＝was
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)。
[0027]
通过公式(4)能够将两幅待检测卫星遥感影像数据分离，得到互相独立的分量，通过对多时相遥感图像进行核主成分分析，得到互不相关的两个分量，即彼此间没有相关信息的分量，实现变化区域向量的提取，获得含变化信息的分量和背景分量，接着将含变化信
息的分量转化为数据矩阵，然后进行非下采样shearlet逆变换得到变化信息图像。
[0028]
步骤4包括：设多时相遥感变化图像中像素(m,n)的灰度级为f(m,n)，其邻域平均灰度级为g(m,n)，h(k)为变化图像中满足总体灰度级f(m,n)＝f(m,n)+g(m,n)＝k的像素点频数；k表示截距；
[0029]
用阈值t将变化图像按灰度级f(m,n)分为变化类ωo和非变化类ωb：
[0030]
ωo＝{(m,n)|f(m,n)＝0,1,...,t}，
[0031]
ωb＝{(m,n)|f(m,n)＝t+1,t+2,...,2l-1}，
[0032]
并设
[0033]
其中x代表横坐标，y代表纵坐标；这里是以图像左上角顶点为原点，x代表图像横横坐标，y代表图像纵坐标；
[0034]
则基于直线截距直方图的多时相遥感变化图像的总倒数灰度熵e(t)为：
[0035][0036]
式中，eo和eb分别表示变化区域的倒数灰度熵和非变化区域的倒数灰度熵，uo(t)表示目标区域的灰度均值，ub(t)表示背景区域的灰度均值；
[0037]
当总倒数灰度熵e(t)取最大时，表明变化区域内像素灰度级大小差异与非变化区域内像素灰度级东西差异之和最小，即变化图像分割的效果最理想，此时对应的灰度级阈值t
*
即为最优阈值：
[0038][0039]
步骤5包括：以多时相遥感变化图像像素的二维联合信息为基础建立直线截距直方图，所述二维联合信息包括灰度级f(m,n)和邻域平均灰度级g(m,n)；
[0040]
以倒数灰度熵作为变化图像的阈值选取准则函数，搜寻最优阈值，并以最优阈值对变化图像进行阈值分割，得到最终的变化检测结果。
[0041]
有益效果：与现有技术相比，本发明的优点在于：
[0042]
(1)独立主成分分析作为一种信源盲分离技术，在去除各信号分量间的相关信息方面展现出优越的性能，针对高阶相关信息也能起到很好的去相关作用。鉴于此，基于ica的方法被引入至多时相卫星影像变化检测领域中，且受到较多的关注。但该方法对图像中存在的噪声等干扰较为敏感，且存在计算量大、收敛速度慢等缺点，需要采取相应优化措施。
[0043]
(2)多尺度几何分析方法不仅能有效地抑制卫星影像中的噪声等干扰，而且能很好地保留其中的边缘及细节信息。近年来提出的无下采样shearlet变换是在shearlet变换
基础上的一种改进，可以近乎最优地对高维几何结构进行稀疏表示，更为准确地描述图像特征。因此，若将无下采样shearlet变换用于多时相卫星影像变化检测可望进一步改善其检测的精确度。
[0044]
(3)另一方面，阈值分割是一类常用的多时相遥感变化图像分割方法，分割结果的好坏直接影响其变化检测效果的优劣。采用基于直线截距直方图倒数灰度熵的阈值分割方法，以倒数灰度熵作为图像的阈值选取准则函数，弥补了最大熵零点处无定义值及不能反映类内灰度均匀性的缺陷，进一步改善分割效果。同时在考虑变化图像像素的灰度级和邻域平均灰度级这二维信息的同时，通过建立直线截距直方图将其转化为一维信息，大幅度提高方法的运行速度。
附图说明
[0045]
下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。
[0046]
图1是本发明的流程示意图；
[0047]
图2是本发明的t0时刻卫星影像(1994年8月意大利elba岛地区卫星影像)；
[0048]
图3是本发明的t1时刻卫星影像(1994年9月意大利elba岛地区卫星影像)；
[0049]
图4是本发明的无下采样shearlet变换示意图；
[0050]
图5是本发明的无下采样shearlet变换频域剖分图；
[0051]
图6是本发明的无下采样shearlet变换频域支撑区间；
[0052]
图7是本发明的ica原理图；
[0053]
图8是本发明的多时相卫星影像变化检测结果。
具体实施方式
[0054]
本发明提供了一种卫星影像变化检测方法，包括如下步骤：
[0055]
步骤1，对不同时期待变化检测的卫星遥感影像进行无下采样shearlet变换操作，得到相应的高频分量和低频分量；
[0056]
步骤2，提取卫星遥感影像变化信息；
[0057]
步骤3，将步骤1得到的高频分量和低频分量分别转化为矩阵形式，并进行无下采样shearlet变换逆变换得到多时相遥感变化图像；
[0058]
步骤4，获得变化检测结果；
[0059]
步骤5，进行阈值分割，得到最终的变化检测结果。
[0060]
步骤1包括：
[0061]
步骤1-1，合成膨胀仿射系统μ
ab
(ψ)表示为：
[0062]
μ
ab
(ψ)＝{ψ
j,l,c
(x)＝|det(a)|
1/2
ψ(b
laj
x-c)}
ꢀꢀꢀꢀꢀꢀ
(1)
[0063]
式中，为尺度参数，l2指平方可积，是实数集合、是整数集合；ψ
j,l,c
是基函数，x是一个变量，为剪切参数，尺度参数和剪切参数分别控制着无下采样shearlet变换的尺度分解数目和方向分解数目；c为平移参数；各向异性膨胀矩阵a、剪切矩阵b均为2
×
2的可逆矩阵，且|detb|＝1，det表示求行列式；f表示任意满足
的函数变量，如果则称μ
ab
(ψ)中的元素为合成小波；b
l
是指剪切矩阵b的l次方；aj指各向异性膨胀矩阵a的j次方；
[0064]
无下采样shearlet变换是合成小波的一种特殊情况，其中无下采样shearlet变换频域支撑区间近似看作大小2
2j
×2j
、斜率2-j
l的梯形对；
[0065]
步骤1-2，多时相卫星影像多尺度分解；
[0066]
步骤1-3，高频分量多方向剖分。
[0067]
步骤1-2包括：采用无下采样金字塔对不同时期的两幅卫星遥感影像分别进行多尺度分解，选取尺度分解层数，得到低频分量和高频分量。
[0068]
步骤1-3包括：构建带方向和尺度变化的meyer迈耶小波窗，分别对多尺度分解后的高频分量进行方向剖分及坐标映射，得到各方向子带的无下采样shearlet变换系数，选取方向剖分级数，即任意1个高频分量均有n1个不同的方向子带；
[0069]
将多时相卫星遥感影像无下采样shearlet变换分解结果分别转化为向量形式，以低频分量和高频分量组成相应的两个分块向量。
[0070]
步骤2包括：设一组观测信号e＝[e1,e2,
…
,en]
t
，它是n维源信号s＝[s1,s2,
…
,sn]
t
的观测值；第i个观测信号ei由n个独立源信号线性混合形成，即：
[0071]ei
＝a
i1
s1+a
i2
s2+
…
+a
in
sn,i＝1,2,
…
,m
ꢀꢀ
(2)
[0072]
其中m为混合矩阵的行数，a
in
为混合矩阵中第i行n列的值，sn为第n个源信号的观测值；
[0073]
用矢量表示公式(2)得到：
[0074]
e＝as
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0075]
式中，a＝[a1,a2,
…
,am]
t
为混合矩阵，am表示第m个混合矩阵基向量；ai是第i个混合矩阵基向量；ica就是仅通过观测数据e估计出混合矩阵a及它的逆矩阵w，最终估计出未知源信号y：
[0076]
y＝we＝was
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)。
[0077]
通过公式(4)能够将两幅待检测卫星遥感影像数据分离，得到互相独立的分量，通过对多时相遥感图像进行核主成分分析，得到互不相关的两个分量，即彼此间没有相关信息的分量，实现变化区域向量的提取，获得含变化信息的分量和背景分量，接着将含变化信息的分量转化为数据矩阵，然后进行非下采样shearlet逆变换得到变化信息图像。
[0078]
步骤4包括：设多时相遥感变化图像中像素(m,n)的灰度级为f(m,n)，其邻域平均灰度级为g(m,n)，h(k)为变化图像中满足总体灰度级f(m,n)＝f(m,n)+g(m,n)＝k的像素点频数；k表示截距；
[0079]
用阈值t将变化图像按灰度级f(m,n)分为变化类ωo和非变化类ωb：
[0080]
ωo＝{(m,n)|f(m,n)＝0,1,...,t}，
[0081]
ωb＝{(m,n)|f(m,n)＝t+1,t+2,...,2l-1}，
[0082]
并设
[0083]
其中x代表横坐标，y代表纵坐标；这里是以图像左上角顶点为原点，x代表图像横横坐标，y代表图像纵坐标；
[0084]
则基于直线截距直方图的多时相遥感变化图像的总倒数灰度熵e(t)为：
[0085][0086]
式中，eo和eb分别表示变化区域的倒数灰度熵和非变化区域的倒数灰度熵，uo(t)表示目标区域的灰度均值，ub(t)表示背景区域的灰度均值；
[0087]
当总倒数灰度熵e(t)取最大时，表明变化区域内像素灰度级大小差异与非变化区域内像素灰度级东西差异之和最小，即变化图像分割的效果最理想，此时对应的灰度级阈值t
*
即为最优阈值：
[0088][0089]
步骤5包括：以多时相遥感变化图像像素的二维联合信息为基础建立直线截距直方图，所述二维联合信息包括灰度级f(m,n)和邻域平均灰度级g(m,n)；
[0090]
以倒数灰度熵作为变化图像的阈值选取准则函数，搜寻最优阈值，并以最优阈值对变化图像进行阈值分割，得到最终的变化检测结果。
[0091]
实施例
[0092]
实现本发明的流程示意图如图1所示，包括如下步骤：
[0093]
步骤1：无下采样shearlet变换。图4给出了无下采样shearlet变换的分解流程，包括多尺度分解和多方向分解(图中用虚线框表示)。前者采用无下采样金字塔(nsp)对图像逐层划分；后者通过改进的剪切滤波器(sf)对每一层高频分量进行方向细分。无下采样shearlet变换的多尺度分解过程通过对图像进行j次nsp分解后，得到与原图像尺寸相同的j+1个分量，包括1个低频分量(最底层)和j个高频分量。无下采样shearlet变换的多方向分解过程主要通过改进的剪切滤波器sf实现，最后经过逆傅里叶变换和二维卷积完成。由此避免了标准剪切滤波器sf中需要进行下采样操作的问题，能有效克服shearlet变换产生的伪吉布斯效应。
[0094]
合成膨胀仿射系统(维数为2)表示为：
[0095]
μ
ab
(ψ)＝{ψ
j,l,c
(x)＝|det(a)|
1/2
ψ(b
laj
x-c)}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0096]
式中，为尺度参数，为剪切参数，分别控制着无下采样shearlet变换的尺度分解数目和方向分解数目；c为平移参数；各向异性膨胀矩阵a、剪切矩阵b均为2
×
2的可逆矩阵，且|detb|＝1。若则称μ
ab
(ψ)中的元素为合成小波。
[0097]
无下采样shearlet变换是合成小波的一种特殊情况，其中无下采样shearlet变换频域支撑区间可以近似看作大小2
2j
×2j
、斜率2-j
l的梯形对，图5、图6分别给出了无下采样shearlet变换的频域剖分图与相应的频域支撑区间。
[0098]
(1)多时相卫星影像多尺度分解。采用无下采样金字塔对不同时期的两幅卫星影像(如图2、图3所示)分别进行多尺度分解，本发明选取尺度分解层数为3，由此得到1个低频分量和3个高频分量；
[0099]
(2)高频分量多方向剖分。构建带方向和尺度变化的meyer小波窗，分别对多尺度分解后的高频分量进行方向剖分及坐标映射，得到各方向子带的无下采样shearlet变换系数。本章选取方向剖分级数为3，即任意1个高频分量均有10个不同的方向子带。
[0100]
将多时相卫星影像无下采样shearlet变换分解结果分别转化为向量形式，以低频分量和高频分量组成相应的两个分块向量。
[0101]
步骤2：基于独立主成分分析的卫星影像变化信息提取。ica方法的目的是，将数据进行某种线性变换，使其变换到相互独立的方向上，各个分量之间不仅要求正交，而且需要相互独立。ica已经成功地应用于盲源信号分离问题，设一组观测信号e＝[e1,e2,
…
,en]
t
，它是n维源信号s＝[s1,s2,
…
,sn]
t
的观测值。第i个观测信号由n个独立源信号线性混合形成，即：
[0102]ei
＝a
i1
s1+a
i2
s2+
…
+a
in
sn,i＝1,2,
…
,m
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0103]
用矢量表示上式则有：
[0104]
e＝as
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0105]
式中，a＝[a1,a2,
…
,am]
t
称为混合矩阵，ai是混合矩阵基向量。ica就是仅通过观测数据e估计出混合矩阵a及它的逆矩阵w，最终估计出未知源信号y：
[0106]
y＝we＝was
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0107]
如果通过学习能够实现wa＝i(其中i是单位矩阵)，则有y＝s，从而达到了源信号的分离。ica原理如图7所示。
[0108]
ica的解并不唯一，y与s也只是在某种判据最优情况下的近似，根据判据的不同，有多种求解ica的方法，本发明通过基于负熵的独立性判决准则的方法求得解混矩阵w，可以有效地把不动点迭代所带来的优良算法特性与负熵所带来更好的统计特性结合起来，以负熵最大作为一个搜寻方向，可以实现顺序地提取独立源，充分体现了投影追踪的思想，此外，该算法采用了定向迭代的优化算法，使得相比传统ica方法具有更快的收敛速度。
[0109]
步骤3：将变化区域的低频分量和高频分量数据向量分别转化为矩阵形式，并进行无下采样shearlet变换逆变换得到多时相遥感变化图像。
[0110]
步骤4：得到变化图像后，需要采取合适的图像分割方法获得变化检测结果，其中阈值分割方法简单、易操作、效果较好。本发明以变化图像像素的灰度级-邻域平均灰度级二维联合信息为基础建立直线截距直方图，以此将阈值搜索空间转化为一维，大大提高分
割的运行效率，并以倒数灰度熵作为阈值选取准则函数，考虑类内灰度均匀性，避免零值无定义值的缺陷，改善方法的分割精度。
[0111]
设多时相遥感变化图像中像素(m,n)的灰度级为f(m,n)，其邻域平均灰度级为g(m,n)，h(k)为该变化图像中满足f(m,n)＝f(m,n)+g(m,n)＝k的像素点频数。则可以用阈值t将该变化图像按灰度级f(m,n)分为变化类ωo和非变化类ωb：
[0112]
ωo＝{(m,n)|f(m,n)＝0,1,...,t}，ωb＝{(m,n)|f(m,n)＝t+1,t+2,...,2l-1}
[0113]
并设：
[0114][0115]
则基于直线截距直方图的多时相遥感变化图像的总倒数灰度熵e(t)为：
[0116][0117]
式中，eo和eb分别表示变化、非变化区域的倒数灰度熵，分别表示变化、非变化区域的倒数灰度熵，其可以反映变化、非变化区域的类内灰度均匀性。当总倒数灰度熵e(t)取最大时，表明变化区域内像素灰度级大小差异与非变化区域内像素灰度级东西差异之和最小，即变化图像分割的效果最理想，此时对应的灰度级阈值t
*
即为最优阈值：
[0118][0119]
步骤5：以多时相遥感变化图像像素的二维联合信息为基础建立直线截距直方图，以倒数灰度熵作为变化图像的阈值选取准则函数，搜寻最优阈值，并以此对变化图像进行阈值分割，得到最终的变化检测结果，如图8所示是本发明的多时相卫星影像变化检测结果。
[0120]
本发明提供了一种卫星影像变化检测方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。