一种基于动态模式分解的交互式光学云检测方法

1.本发明涉及云检测技术领域，特别涉及一种基于动态模式分解的交互式光学云检测方法。


背景技术：

2.随着卫星遥感技术的研究开发和新一代静止轨道卫星观测系统的升级，许多学者对遥感数据的研究与应用工作越来越全面，但由于气候因素的影响，使得遥感领域的研究受到了极大的干扰，遥感数据的研究性大幅度减少。其中，云层的遮挡就是影响因素之一。
3.因此，去除云层的遮挡不仅是为之后进行遥感图像领域的研究打下坚实的基础,也是增强卫星遥感数据有效性、切实性的重要途径之一。而云检测作为云层去除的重要前处理步骤之一，其检测结果的准确性和完整性关系到之后的遥感图像研究结果。同时，由于云自身的多变性和卫星观测图像的复杂性，使得研究高效性和低成本的云检测方法具有十分重要的意义。
4.近年来，随着遥感卫星性能的大幅度提升，以光学遥感图像为研究基础的云检测方法大量涌现。可将其大致分为三类。第一类是通过分析云的纹理结构和调节阈值的传统遥感图像检测方法；第二类是通过研究云的物理属性进行手工设计特征的统计学检测方法；第三类是基于卷积神经网络的深度学习检测方法。
5.其中，第一类云检测方法对卫星观测数据的质量有较高的要求，且基本上都是通过增强卫星云图中云层的细节来提升云检测的性能，没有对云的运动特性进行考虑；第二类云检测方法需要研究学者研究云的物理属性，基于研究成果和研究者本身的研究经验进行手工设计云检测网络结构，其过程复杂性高，且对计算设备性能有一定的要求，存在检测成本高的问题；最后一类检测方法在检测效率和检测精度较前两种检测方法都有一定的提升，但是很难保证在提升检测效率的同时提升检测精度，或者提升很低，且没有考虑添加对检测结果进行修正的操作。如何解决这些问题对于云检测具有关键意义。


技术实现要素：

6.本发明的目的是提供一种简单而有效的基于动态模式分解的交互式光学云检测方法。
7.为此，本发明技术方案如下：
8.一种基于动态模式分解的交互式光学云检测方法，步骤如下：
9.s1、获得用于构建云检测模型的指导数据，其为由基于连续时间段内的n帧卫星云图经真彩色融合处理后得到的n帧真彩色合成图像，进而合成得到的视频m：
10.s2、构建dmd算法以对输入视频中的云进行检测；
11.s201、将由步骤s1得到的卫星云图视频m从m1到mn拉长成时间间隔为δ
t
的n帧均匀
的采样数据，得到时间序列数据矩阵其表达式为：
12.s202、对数据矩阵x1进行灰度化处理，得到r行、c列的首帧灰度图像矩阵；将首帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并分别赋值给矩阵x1和矩阵x；
13.s203、对数据矩阵x2进行灰度化处理，得到r行、c列的第2帧灰度图像矩阵；将第2帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并赋值给矩阵y2；将矩阵x1与矩阵y2进行重构，以更新矩阵x1为：x1＝[x
1 y2]；将矩阵x与矩阵y2进行重构，以更新矩阵x为：x＝[x y2]；
[0014]
s204、对数据矩阵x3至数据矩阵x
n-1
进行与步骤s203相同的处理，以实现不断更新矩阵x1和矩阵x；最终，矩阵x1更新为：x1＝[x
1 y
2 y
3 y
4 y5ꢀ…ꢀyn-1
]，矩阵x更新为：x＝[x y
2 y
3 y
4 y5ꢀ…ꢀyn-1
]；
[0015]
s205、对数据矩阵xn进行灰度化处理，得到r行、c列的最后帧灰度图像矩阵；将最后帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并赋值给矩阵yn；将经过步骤s204得到矩阵x1的第一列元素删除，构建得到矩阵x2，即：x2＝[y
2 y
3 y
4 y5ꢀ…ꢀyn-1
]；将矩阵x2和矩阵yn进行重构，更新矩阵x2为：x2＝[y
2 y
3 y
4 y5ꢀ…ꢀyn
]；将步骤s204得到的矩阵x1与矩阵yn进行重构，更新矩阵x1为：x1＝[x
1 y
2 y
3 y
4 y5ꢀ…ꢀyn
]；将步骤s204得到的矩阵x与矩阵yn进行重构，更新矩阵x为：x＝[x y
2 y
3 y
4 y5ꢀ…ꢀyn
]；
[0016]
s206、将由步骤s205得到的矩阵x1代入公式为：x1＝usv
*
进行奇异值分解，以求取矩阵s的秩r；式中，v
*
表示矩阵v的转置；基于矩阵s的秩r的结果，将矩阵u中第1～r列的全部元素构建为新的矩阵，并重新赋值给矩阵u；将矩阵s中第1～r列的全部元素构建为新的矩阵，并重新赋值给矩阵s；将矩阵v中第1～r列的全部元素构建为新的矩阵，并重新赋值给矩阵v；
[0017]
s207、将经过步骤s207重新赋值后的矩阵v转置后赋值给矩阵v
star
，将经过步骤s207重新赋值后的矩阵u转置后赋值给矩阵u
star
，将经过步骤s207重新赋值后的矩阵s求逆后赋值给矩阵s
star
；进而，结合步骤s205得到的矩阵x2，求得满足xn＝a
×
x
n-1
的近似矩阵a的相似矩阵fdmd：f
dmd
＝u
×
x2×v×
s，以及相似矩阵f
dmd
特征值的对角矩阵ddmd和特征向量矩阵ydmd；继而，通过公式：faii＝u
×
ydmd，求得矩阵faii；
[0018]
s208、将矩阵f
dmd
的行数赋值给afdmd，将矩阵f
dmd
的列数赋值给bfdmd，将矩阵v
star
的列数赋值给nn，构建一个行数
×
列数为afdmd
×nn
的空矩阵vand，以及一个包含对角矩阵ddmd主对角线上向量元素的矩阵edmd；
[0019]
s209、将矩阵edmd的元素进行零次幂运算后赋值至矩阵vand的第一列，将矩阵edmd的元素进行一次幂运算后赋值至矩阵vand的第二列；将矩阵edmd的元素进行二次幂运算后赋值至矩阵vand的第三列；继续重复上述操作，直到将矩阵edmd的元素进行nn次幂运算后赋值至矩阵vand的第nn列，完成矩阵vand的元素赋值；
[0020]
s210、将矩阵ydmd赋值给矩阵l，将矩阵vand赋值给r；
[0021]
s211、将矩阵l和矩阵r代入公式：p＝(l
’×
l)
×
conj.(r
’×
r)中，求得矩阵p；式中，l’表示对矩阵l求逆，r’表示对矩阵r求逆，conj.(
·
)表示对矩阵(
·
)中每个元素进行复共轭；
[0022]
s212、通过公式：q＝conj.(diag(r
×v×s’×
l))，求得矩阵q；式中，s’表示对矩阵s求逆，diag(r
×v×s’×
l)表示先获得r
×v×s’×
l的主对角线元素并重新组成列向量，conj.(
·
)表示对上述得到的列向量(
·
)中每个元素进行复共轭；基于矩阵q通过公式：
aerfa＝diag(p-1
×
q)，求得矩阵aerfa；
[0023]
s213、将步骤s205得到的矩阵x1的第一列赋值给矩阵x
列
，并通过振幅计算公式：计算得到矩阵faii的矩阵将矩阵中主对角线上元素组成列向量，并赋值给矩阵
[0024]
s214、构建行数
×
列数为afdmd
×
(nn+1)的全零矩阵vand1；将矩阵edmd的元素进行零次幂运算后赋值至矩阵vand1的第一列，将矩阵edmd的元素进行一次幂运算后赋值至矩阵vand1的第二列；将矩阵edmd的元素进行二次幂运算后赋值至矩阵vand1的第三列；继续重复上述操作，直到将矩阵edmd的元素进行(nn+1)次幂运算后赋值至矩阵vand1的第(nn+1)列，完成矩阵vand1的元素赋值；
[0025]
s215、将矩阵faii的第一列元素赋值给矩阵faiia，将矩阵aerfa的第一行第一列的元素赋值给aerfaa，将矩阵vand1的第一行元素赋值给vand
1a
；通过公式：a_hat＝faiia×
aerfaa×
vand
1a
，求得矩阵a_hat；
[0026]
s216、基于步骤s205获得的矩阵x的行数row和列数col，构建行数
×
列数为row
×
col的全零矩阵e_hat；并利用公式：e_hat(i,j)＝x(i,j)-a_hat(i,j)完成对矩阵e_hat的重新赋值：
[0027]
s217、对矩阵e_hat的第一列元素至第row列元素的每列元素进行如下相同方式处理：将指定列元素构建为矩阵，基于首帧图片的宽度img_width和高度img_height，将矩阵重构为行数
×
列数为img_height
×
img_width的矩阵，并写入jpg文件并保存；经过该步骤依次获得矩阵d1、d2、
…
、d
row
；
[0028]
s218、通过公式：d＝d1+d2+d3+
…
+d
row
，对矩阵d1、d2、
…
、d
row
进行融合，得到矩阵d；
[0029]
s219、采用otsu方法对经过上述步骤s218得到的灰度图像矩阵d进行二值化处理，所得二值图像即为云检测结果图；
[0030]
s3、对比首帧原始图像与云检测结果图，采用交互修改的方式在云检测结果图中补全漏检云。
[0031]
进一步地，在步骤s1中，图像的真彩色融合处理方法为：
[0032]
步骤1：分别获取卫星云图中波长为0.65μm图像的反射率dn和波长为0.82μm图像的反射率dn，以计算二者的比值；
[0033]
步骤2：当时，则判定该图像的植被丰富，对波长为0.82μm的图像与波长为0.65μm的图像的计数值进行互换后，将波长为0.65μm、0.47μm和0.82μm的三通道图像合成得到的彩色图像作为真彩色合成图像；当时，则判定该图像的植被不丰富，则采用原始卫星云图作为真彩色合成图像；k为调节图像显示效果的阈值参数。
[0034]
进一步地，在步骤s1中，n为16。
[0035]
进一步地，步骤s3的具体实施步骤为：
[0036]
s301、将经过步骤s2得到的云检测结果图与视频m中的首帧图像进行观察比较，判断云检测结果图中是否存在明显的云漏检问题：
[0037]
1)若存在明显的云漏检问题，则继续步骤s302～步骤s305；
[0038]
2)若不存在明显的云漏检问题，则步骤s2得到的云检测结果图即为最终的得到的云检测结果图；
[0039]
s302、将云检测结果图的分辨率改为815
×
471，并保存为png格式；将首帧图像的分辨率改为815
×
471，并保存为jpg格式；
[0040]
s303、参照检测结果图，在首帧图像中标注出云检测图中漏检云的边框，并将边框所在像素以及位于边框内的全部像素标记为前景像素；进而，根据前景像素的像素值，在首帧图像的其它像素中遍历寻找与前景像素具有相同像素值的全部像素，同样标记为前景像素；将首帧图像中剩余像素均标记为后景像素；
[0041]
s304、基于步骤s303的标记，在首帧图像的图像矩阵中，保留全部前景像素的像素值，而将后景像素的像素值重新赋值为0，实现对前景像素和后景像素的分割，进而得到图像矩阵z1；
[0042]
s305、采用cv2.bitwise_and函数，对首帧图像的图像矩阵与由步骤s220得到的二值化图像矩阵进行融合，得到图像矩阵z2；采用cv2.add函数，对图像矩阵z1与图像矩阵z2进行融合，得到新图像矩阵z3；
[0043]
s306、重复步骤s301～步骤s305，直至云检测结果图与首帧图像之间无明显的云漏检问题。
[0044]
与现有的云检测方法相比，该基于动态模式分解的交互式光学云检测方法的有益效果在于：
[0045]
1)与传统的云检测方法相比，本发明在处理图像数据中利用奇异值求解，保留图像关键特征点的同时尽量减小观测数据规模，简单高效的对卫星云图进行检测，降低了维度，避免了云检测过程中较大的运算量，总体上大大减小了计算量，避免了云检测方法对硬件的超高要求，经济成本更低；
[0046]
2)传统云检测方法未顾及云的形状和位置时刻处于变化的特点，直接从卫星观测图像中提取云层特征进行云检测，降低了云检测精度；而本技术突破了传统方法，利用动态模态分解算法实现检测运动物体的特性；在对运动的云层进行运动前后景分割检测的同时，考虑到存在一定时间段内静止不动的云无法被检测到导致云检测结果存在误差，提出可视化交互修正方案，对检测误差进行修正，提高云检测结果；
[0047]
3)本技术的方法解决了现有技术中基于纹理分析法进行云检测难以有效的拉开云层与地物在空间上的分布距离，难以选取合适的空间特征将云与地物加以区分，以及对检测结果中存在的漏判进行修正的技术问题，且该方法无须人工标注样本，无须大量的先验知识，也不需要像机器学习那样长时间的复杂训练过程，有效降低了使用门槛并提高了实用性，是一种简单而有效的云检测方法。
附图说明
[0048]
图1为本发明的基于动态模式分解的交互式光学云检测方法的流程图；
[0049]
图2为本发明的基于动态模式分解的交互式光学云检测方法中步骤s2的dmd流程图；
[0050]
图3(a)为本发明的基于动态模式分解的交互式光学云检测方法中步骤s3的交互界面图1；
[0051]
图3(b)为本发明的基于动态模式分解的交互式光学云检测方法中步骤s3的交互界面图2；
[0052]
图4(a)为不同类型卫星云图原始图像以及采用本发明的基于动态模式分解的交互式光学云检测方法的步骤s1和步骤s2处理得到的云检测结果图；
[0053]
图4(b)为一种层积云卫星云图原始图像以及采用本发明的基于动态模式分解的交互式光学云检测方法的步骤s1至步骤s3处理得到的云检测结果图；
[0054]
图4(c)为另一种层积云卫星云图原始图像以及采用本发明的基于动态模式分解的交互式光学云检测方法的步骤s1至步骤s3处理得到的云检测结果图。
具体实施方式
[0055]
下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。
[0056]
参见图1，基于现有的2021年6月15日的风云四号fy-4a成像仪l1观测的4km的中国区域的hdf数据，采用本技术的基于动态模式分解的交互式光学云检测方法对该日的部分连续时间段内观测的卫星云图进行云检测，其具体实施步骤如下：
[0057]
s1、将静止气象卫星风云四号的多通道扫描成像辐射计(简称为agri)获取的14通道的卫星云图进行处理，并合成视频，以获得构建云检测模型的指导数据；
[0058]
具体地，该步骤s1的实施步骤为：
[0059]
s101、获取连续时间段内的n(n＝16)帧卫星云图在波长分别为0.47μm、0.65μm和0.82μm的图像，并对每帧卫星云图进行处理，以获得符合人眼视觉效果的真彩色融合图像；其中，融合图像的获取方法为：
[0060]
步骤1：分别获取卫星云图中波长为0.65μm的图像和波长为0.82μm的反射率dn，以计算得到二者的比值，即rvi比值植被指数：
[0061]
步骤2：当时，则判定该图像的植被丰富，对波长为0.82μm的图像与波长为0.65μm的图像的计数值(即图像上各像素的颜色值)进行互换；当时，则判定该图像的植被不丰富，采用原始卫星云图作为真彩色合成图像；在上式中，dn为特定通道的波段反射率；k为调节图像显示效果的阈值参数，针对不同的仪器，k的取值可以不同；
[0062]
其中，步骤s101的图像校正基于：0.65μm波长对应合成图像的绿色(g)通道；0.47μm波长对应合成图像的蓝色(b)通道；0.82μm波长对应合成图像的红色(r)通道，该通道远离人眼的红色光谱响应，却对植被很敏感的特性；因此，当检测区域在0.82μm通道发现植被丰富区域时，直接交换原始真彩色合成图像的红绿通道计数值，可以获得目视上的真彩色效果；而当植被不丰富时则合成红色通道；进而最终利用不符合人眼红绿蓝视觉特性的3个反射通道最终合成符合人眼视觉效果的真彩色合成图像；
[0063]
s102、将由步骤s101获得的n(n＝16)帧真彩色融合图像，按照时序合成卫星云图视频m：卫星云图的视频格式为avi；
[0064]
s2、构建dmd算法以对输入视频中的云进行检测；
[0065]
具体地，如图2所示，该步骤s2的实施步骤为：
[0066]
s201、将由步骤s1得到的卫星云图视频m，即：从m1到m
16
拉长成时间间隔为δ
t
的16帧均匀的采样数据，得到时间序列数据矩阵，其表达式为：
[0067]
s202、对视频m的首帧彩色图像的数据矩阵x1进行灰度化处理，得到r行、c列的首帧灰度图像矩阵；将首帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并分别赋值给矩阵x1和矩阵x，此时，矩阵x1与矩阵x相同；
[0068]
s203、将视频m的第2帧彩色图像的数据矩阵x2的进行灰度化处理，得到r行、c列的第2帧灰度图像矩阵；将第2帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并赋值给矩阵y2；将矩阵x1与矩阵y2进行重构，更新矩阵x1为：x1＝[x
1 y2]；将矩阵x与矩阵y2进行重构，更新矩阵x为：x＝[x y2]；
[0069]
s204、对视频m中的第3帧彩色图像至第15帧彩色图像(共13幅彩色图像)的数据矩阵x3至x
15
进行与步骤s203相同的处理，以实现不断更新矩阵x1和矩阵x；具体地，经过步骤s204，
[0070]
矩阵x1更新为：
[0071]
x1＝[x
1 y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15
]，
[0072]
矩阵x更新为：
[0073]
x＝[x y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15
]；
[0074]
s205、将最后一帧彩色图像，即第16帧彩色图像的数据矩阵xn进行灰度化处理，构建得到r行、c列的最后帧灰度图像矩阵；将最后帧灰度图像矩阵重构为r
×
c行、1列的新矩阵，并赋值给矩阵y
16
；将经过步骤s204得到矩阵x1的第一列元素删除，构建得到矩阵x2，即：x2＝[y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15
]；将矩阵x2和矩阵y
16
进行重构，更新矩阵x2为：x2＝[y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15 y
16
]；同时，将步骤s204得到的矩阵x1，即x1＝[x
1 y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15
]与矩阵y
16
进行重构，更新矩阵x1为：x1＝[x y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15 y
16
]；同时，将步骤s204得到的矩阵x，即x＝[x y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15
]与矩阵y
16
进行重构，更新矩阵x为：x＝[x y
2 y
3 y
4 y
5 y
6 y
7 y
8 y
9 y
10 y
11 y
12 y
13 y
14 y
15 y
16
]；
[0075]
s206、将由步骤s205得到的矩阵x1代入公式为：x1＝usv
*
进行奇异值(svd)分解，以实现保留图片关键特征点的同时尽量减小观测数据规模，进而求得矩阵s的秩r＝18；式中，v
*
表示矩阵v的转置；
[0076]
基于矩阵s的秩r＝18的结果，为了后续加快运算速度，将矩阵u中第1～18列的全部元素构建为新的矩阵，并重新赋值给矩阵u，完成对矩阵u的更新；将矩阵s中第1～18列的全部元素构建为新的矩阵，并重新赋值给矩阵s，完成对矩阵s的更新；将矩阵v中第1～18列的全部元素构建为新的矩阵，并重新赋值给矩阵v，完成对矩阵v的更新；
[0077]
s207、将经过步骤s206重新赋值后的矩阵v转置后赋值给矩阵v
star
，将经过步骤s207重新赋值后的矩阵u转置后赋值给矩阵u
star
，将经过步骤s206重新赋值后的矩阵s求逆后赋值给矩阵s
star
；进而，结合步骤s205得到的矩阵x2，求得满足xn＝a
×
x
n-1
的近似矩阵a的相似矩阵fdmd：f
dmd
＝u
×
x2×v×
s，以及相似矩阵fdmd特征值的对角矩阵ddmd和特征向量
矩阵ydmd；继而，通过公式：faii＝u
×
ydmd，求得矩阵faii；
[0078]
s208、将矩阵f
dmd
的行数赋值给afdmd，将矩阵f
dmd
的列数赋值给bfdmd，将矩阵v
star
的列数赋值给nn，构建一个行数
×
列数为afdmd
×nn
的空矩阵vand，以及一个包含对角矩阵ddmd主对角线上向量元素的矩阵edmd；
[0079]
s209、将矩阵edmd的元素进行零次幂运算后赋值至矩阵vand的第一列，将矩阵edmd的元素进行一次幂运算后赋值至矩阵vand的第二列；将矩阵edmd的元素进行二次幂运算后赋值至矩阵vand的第三列；继续重复上述操作，直到将矩阵edmd的元素进行nn次幂运算后赋值至矩阵vand的第nn列，完成矩阵vand的元素赋值；
[0080]
s210、将矩阵ydmd赋值给矩阵l，将矩阵vand赋值给r；
[0081]
s211、将矩阵l和矩阵r代入公式：p＝(l
’×
l)
×
conj.(r
’×
r)中，求得矩阵p；式中，l’表示对矩阵l求逆，r’表示对矩阵r求逆，conj.(
·
)表示对矩阵(
·
)中每个元素进行复共轭；
[0082]
s212、通过公式：q＝conj.(diag(r
×v×s’×
l))，求得矩阵q；式中，s’表示对矩阵s求逆，diag(r
×v×s’×
l)表示先获得r
×v×s’×
l的主对角线元素并重新组成列向量，conj.(
·
)表示对上述得到的列向量(
·
)中每个元素进行复共轭；基于矩阵q通过公式：aerfa＝diag(p-1
×
q)，求得矩阵aerfa；
[0083]
s213、将步骤s205得到的矩阵x1的第一列赋值给矩阵x
列
，并通过振幅计算公式：计算得到矩阵faii的矩阵进而，将矩阵中主对角线上元素组成列向量，并赋值给矩阵
[0084]
s214、构建行数
×
列数为afdmd
×
(nn+1)的全零矩阵vand1；而后，将矩阵edmd的元素进行零次幂运算后赋值至矩阵vand1的第一列，将矩阵edmd的元素进行一次幂运算后赋值至矩阵vand1的第二列；将矩阵edmd的元素进行二次幂运算后赋值至矩阵vand1的第三列；继续重复上述操作，直到将矩阵edmd的元素进行(nn+1)次幂运算后赋值至矩阵vand1的第(nn+1)列，完成矩阵vand1的元素赋值；
[0085]
s215、将矩阵faii的第一列元素赋值给矩阵faiia，将矩阵aerfa的第一行第一列的元素赋值给aerfaa，将矩阵vand1的第一行元素赋值给vand
1a
；进而，通过公式：a_hat＝faiia×
aerfaa×
vand
1a
，求得用于重构图像的背景的矩阵a_hat；
[0086]
s216、基于步骤s205获得的矩阵x的行数row和列数col，构建行数
×
列数为row
×
col的全零矩阵e_hat；并利用公式：e_hat(i,j)＝x(i,j)-a_hat(i,j)完成对矩阵e_hat的重新赋值；
[0087]
具体地，上述公式的具体操作过程为：首先，用矩阵x中行数为1、列数为1的元素数值减去矩阵a_hat中行数为1、列数为1的元素数值，并将所得差值赋值给矩阵e_hat中行数为1、列数为1的元素；接着，用矩阵x中行数为1、列数为2的元素数值减去矩阵a_hat中行数为1、列数为2的元素数值，并将所得差值赋值给矩阵e_hat中行数为1、列数为2的元素；以此类推，直至用矩阵x中行数为row、列数为col的元素数值减去矩阵a_hat中行数为row、列数为col的元素数值，并将所得差值赋值给矩阵e_hat中行数为row、列数为col的元素，即最终完成对矩阵e_hat中全部元素的重新赋值，得到更新后的矩阵e_hat；
[0088]
s217、为获得视频m中16张图像的前景，对矩阵e_hat的第一列元素至第row列元素
的每列元素进行如下相同方式处理：将指定列元素构建为矩阵，基于首帧图片的宽度img_width和高度img_height，将矩阵重构为行数
×
列数为img_height
×
img_width的矩阵，并写入jpg文件并保存；
[0089]
该步骤的具体实施过程如下：
[0090]
将矩阵e_hat的第1列的元素构建为矩阵d1，并获得首帧图片的宽度img_width和高度img_height，进而将矩阵d1重构为行数
×
列数为img_height
×
img_width的矩阵，并将矩阵d1写入jpg文件并保存；将矩阵e_hat的第2列的元素构建为矩阵d2，并获得首帧图片的宽度img_width和高度img_height，进而将矩阵d2重构为行数
×
列数为img_height
×
img_width的矩阵，并将矩阵d2写入jpg文件并保存；
[0091]
以此类推，直至将矩阵e_hat的第row列的元素构建为矩阵d
row
，并获得首帧图片的宽度img_width和高度img_height，进而将矩阵d
row
重构为行数
×
列数为img_height
×
img_width的矩阵，并将矩阵d
row
写入jpg文件并保存；、
[0092]
经过该步骤s217，依次获得矩阵d1、d2、
…
、d
row
；
[0093]
s218、为提高图像的信噪比，通过公式：d＝d1+d2+d3+
…
+d
row
，将步骤s218依次获得的矩阵d1、d2、
…
、d
row
进行融合得到矩阵d；
[0094]
s219、采用otsu方法，对经过上述步骤s218得到的灰度图像矩阵d进行二值化处理，所得二值化图像，该二值化图像即为首帧图像的云检测结果图；在该步骤中，所得二值化图像的全局阈值t为0.4078；
[0095]
s3、由于步骤s2检测的是在一段连续的时间内运动的云，因此在此时间内不可避免的存在静止不动的云或者运动特征不明显的云，即仅采用上述步骤s2得到的云检测结果可能和实际结果存在一定的误差；因此，还需要通过判断云检测图与原始图像是否存在明显的差异，进而进行进一步的云检测结果图进行交互式修改步骤，以补全漏检测到的云，提升云检测结果的完整性；
[0096]
具体地，该步骤s3的具体实施步骤如下：
[0097]
s301、将经过步骤s2得到的云检测结果图与视频m中的首帧图像进行观察比较，判断云检测结果图中是否存在明显的云漏检问题：
[0098]
1)若存在明显的云漏检问题，则继续步骤s302～步骤s305；
[0099]
2)若不存在明显的云漏检问题，则步骤s2得到的云检测结果图即为最终的得到的云检测结果图；
[0100]
s302、将云检测结果图的分辨率改为815
×
471，并保存为png格式；将首帧图像的分辨率改为815
×
471，并保存为jpg格式；
[0101]
s303、参照检测结果图，在首帧图像中标注出云检测图中漏检云的边框，并将边框所在像素以及位于边框内的全部像素标记为前景像素；进而，根据前景像素的像素值，在首帧图像的其它像素中遍历寻找与前景像素具有相同像素值的全部像素，同样标记为前景像素；将首帧图像中剩余像素均标记为后景像素；
[0102]
s304、基于步骤s303的标记，在首帧图像的图像矩阵中，保留全部前景像素的像素值，而将后景像素的像素值重新赋值为0，实现对前景像素和后景像素的分割，进而得到图像矩阵z1；
[0103]
s305、采用cv2.bitwise_and函数，对首帧图像的图像矩阵与由步骤s219得到的二
值化图像矩阵进行融合，得到图像矩阵z2；采用cv2.add函数，对图像矩阵z1与图像矩阵z2进行融合，得到新图像矩阵z3；
[0104]
s306、重复步骤s301～步骤s305，直至云检测结果图与首帧图像进之间无明显的云漏检问题。
[0105]
在本实施例的具体实践中，上述步骤基于grabcut算法实现，其具体利用在python软件中将现有的grabcut算法编写为计算程序，以实现通过计算机对步骤s2得到的云检测结果图进行一次或多次交互式修改，以补全检测到的云，提升云检测结果的完整性。
[0106]
如图3(a)所示，计算机界面上共展示有分别位于左上、左下和右下的三幅图像；其中，位于左上的图像为由cloud detection first result diagram窗口显示的经过步骤s2的云检测方法得到的首帧云检测结果图，位于左下的图像为由the original image窗口展示的视频m的首帧原始图像，该两幅图像为后面进行误差修正提供参考；位于右下的图像为由the process diagram窗口显示的进行交互操作后得到的过程标注图；在该时刻下，由于交互式修改还没有开始，因此，在图3(a)中，该图像为空白状态；另外，在后续的交互式修改操作中，计算机界面的右上角还将展示由cloud detection second result diagram窗口显示的进行交互式修改后得到的云检测结果图。
[0107]
如图3(b)所示为采用本发明的基于动态模式分解的交互式光学云检测方法得到的云检测结果图；图中显示计算机界面上共展示有分别位于左上、左下。左上和右下共四幅图像；其中，位于左上的图像为由cloud detection first result diagram窗口显示的经过步骤s2的云检测方法得到的首帧云检测结果图，位于左下的图像为由the original image窗口展示的视频m的首帧原始图像，该两幅图像为后面进行误差修正提供参考；位于右下的图像为由the process diagram窗口显示的进行交互操作后得到漏检云图像；位于右上的图像为由cloud detection second result diagram窗口显示的进行交互式修改后得到的云检测结果图；从该右上图像可以看出，通过交互式修改后，在左上的云检测结果图的基础上完成了漏检云结果的修改。
[0108]
为了进一步验证本技术对卫星云图处理的广适性，基于本技术的方法，对各类卫星云图进行处理并展示如下。
[0109]
如图4(a)所示为采用本发明的基于动态模式分解的交互式光学云检测方法的步骤s1和步骤s2得到的原始图像与云检测结果图的对比图像，涉及的云图类型包括絮状云、波状云、台风云和层积云四种类型；其中，由于在对上述卫星云图的处理中，选用的多帧原始图像所在连续时间段内不存在静止不动的云或者运动特征不明显的云，因而不存在云漏检的现象，因此，图4(a)中涉及的原始图像仅经过本技术的步骤s1和步骤s2的处理即可，不需要进行步骤s3，即不存在漏检问题得到图像。
[0110]
如图4(b)和图4(c)所示均为采用本发明的基于动态模式分解的交互式光学云检测方法的步骤s1、步骤s2和步骤s3依次得到的原始图像、云检测结果图、以及经过交互修改后的云检测结果图；由于在对上述卫星云图的处理中，选用的多帧原始图像所在连续时间段内有静止不动的云或者运动特征不明显的云的存在，因而必然存在云漏检的现象，因而还需要通过步骤s3进行交互修改以在云检测结果图中补全漏检的云；具体地，在图4(b)中，采用本技术处理的云图类型为层积云，其原始图像中标注的方框内的是在步骤s3中发现的运动特征不明显的漏检云，因此通过步骤s3对漏检云进行了补充；在图4(c)中，采用本技术
处理的云图类型为层积云，其原始图像中标注的方框内的是在步骤s3中发现的静止不动的漏检云，因此通过步骤s3对漏检云进行了补充。
[0111]
综上，采用本技术的云检测方法能够实现对不同类型的云和处于不同地面环境的云都可以检测出来，从直观的视觉效果分析，其检测的云准确度和完整度很高。