一种批量遥感影像预处理方法与流程

本发明属于遥感及影像处理技术领域，具体涉及一种批量遥感影像预处理的方法。

背景技术：

近年来，遥感影像在成像过程中受到天空中的云的影响，影响地物的判读与解译。对于遥感影像信息提取而言，含云影像将大大降低遥感影像的利用率，特别是给遥感影像预处理工作带来极大的不便，加大了信息提取处理的难度，降低了遥感影像预处理的效率。

从上世纪九十年代以来，国内不少遥感领域的科研工作者先后进行了这方面的研究。采用阈值法、小波纹理分析法、面向对象法、神经网络法开展了一系列关于云检测的工作。

如2010年采用阈值法，根据云和地物的反射率随波长的变化而变化实现FORMOSAT-2,VENDS,LANDSAT and SENTINEL-2影像云地分离，但是不能保证阈值对每景影像上的云都能达到很好的识别效果。2007年采用无抽样小波变换对全色遥感影像进行薄云去除，针对单景遥感影像效果较好，对于批量遥感数据而言，效率不高，适应性较弱。2011年采用小波SCM方法提取纹理特征识别云，取得了一定的成果，在云雪混淆的时候适应性较好，但受到传感器类型不同，同时由于影像上各种地物繁杂，空间分布并不集中，特别在影像地物类别较多的情况下，此法识别出的云的效果略差。2012年采用面向对象的Fmask算法检测云，算法精度较高，但云和冰混淆的部分难以区分，特别是在批量遥感影像处理中，效率不高。采用支持向量机实现对遥感影像厚云的检测与去除，并采用统计学补偿的方法进行修复，精度相对较高，但是在薄云的检测方面相对较弱。2015年采用BP神经网络法进行云检测，因其较强的学习能力，对样本的理解与识别而致使云判效果较好，但是运算量过大，不适合批量遥感影像处理。上述工作都不同程度地取得了一些成就，但工作方法相对单一，效果时好时坏。对云边界识别精度不准不仅导致了大量的无效影像被检索，造成巨大浪费。因此对批量遥感影像预处理云的检测工作已成为当务之急。

技术实现要素：

为解决上述问题，本发明公开了一种批量遥感影像预处理方法，包括以下步骤：

步骤1，批量获取遥感影像并建立批量遥感影像数据库，通过统计遥感影像的光谱特征和纹理特征获取含云影像；

步骤2，检测含云影像上的云域范围及其含云量，去除含云量超过含云阈值的含云影像，得到剩余含云影像；

步骤3，去除步骤2中剩余含云影像上的云，得到无云影像；

步骤4，利用得到的无云影像，更新批量遥感影像数据库。

进一步地，步骤1中所述的通过统计遥感影像的光谱特征获取含云影像是指：

步骤11，在批量遥感影像数据库中任选一遥感影像作为当前遥感影像，获取该当前遥感影像在红波段的灰度频率直方图；

步骤12，若所述的灰度频率直方图中出现两个峰值频率，随着灰度值的增大依次为第一峰值频率和第二峰值频率，且两个峰值之间存在一个频率相对平缓且低于两个峰值的缓冲区，同时第二峰值频率明显高于第一峰值频率，则将该第二峰值频率对应的灰度值标记为灰度异常值且将该当前遥感影像标记为含云影像；否则，将该当前遥感影像标记为待定影像；

步骤13，重复步骤11-12，直至批量遥感影像数据库中所有的遥感影像都被标记为止。

进一步地，步骤1中所述的通过统计遥感影像的纹理特征获取含云影像是指：

步骤14，从待定影像中任取一个作为当前待定影像，对当前待定影像进行直方图均衡化处理，分别统计当前待定影像在直方图均衡化处理前后的均值和二阶矩；

步骤15，若满足且-0.03≤ASM_前-ASM_后≤0.03，则将该当前待定影像标记为含云影像，否则标记为无云影像；

其中，Mean_后为待定影像均衡化处理后的均值，Mean_前为待定影像均衡化处理前的均值，ASM_前为待定影像均衡化处理前的二阶矩，ASM_后为待定影像均衡化处理后的二阶矩；

步骤16，重复步骤14、步骤15，直至所有的待定影像都被标记为止。

进一步地，步骤2中所述的检测含云影像上的云域范围及该云域范围的含云量，去除含云量超过含云阈值的含云影像是指：

步骤21，通过建立决策树云检测规则依次检测步骤2中获取的所有含云影像中的云域范围，所述的决策树检测规则如下：

其中，a为云在红波段的灰度异常值，x为灰度值；

步骤22，分别统计含云影像和云域范围的像素个数，通过下式计算云域范围的含云量：

其中，C为含云量，N_c为云域范围的像素个数，N_I为含云影像的像素个数；

步骤23，设定含云阈值D,若含云影像中的含云量C大于D时，去除该含云影像。

进一步地，步骤21中所述的云在红波段的灰度异常值a的计算方法为：

步骤211，对含云影像的灰度直方图中的频率做一阶向前差分；

步骤212，设定阈值L，若前后频率差分值＞L时，表明该频率段不在缓冲区；否则，该频率段为疑似缓冲区；

取最长的一段疑似缓冲区为频率缓冲区；

步骤213，通过冒泡法得到含云影像的灰度直方图上的可能异常值，计算频率缓冲区的频率均值w，若该可能异常值大于2w，则该可能异常值即为云在红波段的灰度异常值a。

进一步地，步骤3中所述的去除步骤2中剩余含云影像上的云，得到无云影像是指：

步骤31，依次选取含云影像作为含云影像A，从批量遥感影像数据库中选取与该含云影像A的经纬度相同但不同时相的所有无云影像,从所述无云影像中选择质量好且分辨率与该含云影像A相同的无云影像B；

步骤32，分别裁剪出所述含云影像A与无云影像B的云域区域C1、C2，用C2区域替换掉C1区域，镶嵌得到中间影像A1；

步骤33，对中间影像A1采用均值滤波处理得到新的遥感影像A2，即无云影像；

进一步地，步骤33中所述的对中间影像A1采用均值滤波处理得到新的遥感影像A2是指：

步骤331，选取中间影像A1中拼接线处的像素作为中心像素，并找到该中心像素8邻域的8个像素，求出这8个像素和中心像素的均值作为新的中心像素；

步骤332，以该新的中心像素为中心，分别向内外依次各缓冲5个像素；

步骤333，对中间影像A1中拼接线处所有像素都进行步骤331、步骤332后，则得到新的遥感图像A2。

与现有技术相比，本发明具有以下技术效果：

1.本发明快速高效，精度较高，提供了一种针对大批量遥感影像预处理的方法，极大程度上避免了数据过度浪费，具有较强的实用价值和推广意义；

2.本发明在影像预处理的过程中，确定了影像上的含云量，自动圈定了影像上的云域范围并去除，对影像信息提取获取数据给出了准确的依据；

3.本发明采用了工作效率较高的决策树方法，不仅效率高，而且人为干预少，减轻了人工云量判读的工作量，缩短了云检测和去除的工作时间，可用于批量选取遥感影像及遥感影像信息提取研究。

附图说明

图1为本发明总体流程图；

图2为本实施例中云在红波段的灰度频率直方图；

图3为原始初级云检测规则；

图4为资源三号卫星影像常见的含云影像图，其中(a)为含点云的影像图，(b)为含卷云的影像图，(c)为含厚云的影像图，(d)为含薄云的影像图；

图5为本实施例经过云检测后的影像图，其中(a)为经过云检测后含点云的影像图，(b)为经过云检测后含卷云的影像图，(c)为经过云检测后含厚云的影像图，(d)为经过云检测后含薄云的影像图；

图6为本实施例中均值滤波处理的示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明。

如图1所示，本发明公开了一种批量遥感影像预处理方法，具体包括以下步骤：

步骤1，批量获取遥感影像并建立批量遥感影像数据库，通过统计遥感影像的光谱特征和纹理特征获取含云影像；

根据遥感影像的轨道号及传感器类型，分别建立遥感影像数据库；

不同的传感器往往导致遥感影像的分辨率不同，因此根据影像的轨道号及传感器载荷分别建立遥感影像数据库。如表1对资源三号卫星影像按照不同的传感器及轨道号建立关系型数据库的方式存储于影像数据库中。

表1资源三号卫星影像数据库

不同的遥感影像上云的光谱特征及纹理特征表现相似。根据云在红波段反射最明显挑选出部分典型含云影像，并统计含云影像的光谱特征及纹理特征。

其中，通过统计遥感影像的光谱特征获取含云影像是指：

步骤11，在批量遥感影像数据库中任选一遥感影像作为当前遥感影像，获取该当前遥感影像在红波段的灰度频率直方图；

步骤12，若所述的灰度频率直方图中出现两个峰值频率，随着灰度值的增大依次为第一峰值频率和第二峰值频率，且两个峰值之间存在一个频率相对平缓且低于两个峰值的缓冲区，同时第二峰值频率明显高于第一峰值频率，则将该第二峰值频率对应的灰度值标记为灰度异常值且该当前遥感影像标记为含云影像；否则，将该当前遥感影像标记为待定影像；

步骤13，重复步骤11-12，直至批量遥感影像数据库中所有的遥感影像都被标记为止。

云的光谱特征可以通过灰度直方图表征，云在红波段反射是最明显的，本实施例中，如图2所示，无云影像的灰度直方图呈现近似的标准正态分布。除薄云外，点云、厚云、卷云的灰度直方图均会显示出两个明显的峰值，随着灰度值的增大依次为第一峰值频率和第二峰值频率，且两个峰值之间存在一个频率相对平缓且低于两个峰值的缓冲区，同时第二峰值频率明显高于第一峰值频率。本实施例中，如图2所示，以2014年12月7日获取，轨道号016170，Path/Row为894/127的资源三号卫星影像红波段的灰度直方图，在通过对影像上的已有的云进行灰度值统计，图2中的圆圈圈定的峰值横坐标就是云在遥感影像上的光谱异常值a，平滑部分则为两个峰值之间存在的缓冲区。

对于全色单波段影像，则可以直接获取如图2所示的灰度直方图异常值，就可以标记出含云影像。

其中，通过统计遥感影像的纹理特征获取含云影像是指：

步骤14，从待定影像中任取一个作为当前待定影像，对当前待定影像进行直方图均衡化处理，分别统计待定影像在直方图均衡化处理前后的均值和二阶矩；

步骤15，若满足且-0.03≤ASM_前-ASM_后≤0.03，则将该当前待定影像标记为含云影像，否则标记为无云影像；

其中，Mean_后为待定影像均衡化处理后的均值，Mean_前为待定影像均衡化处理前的均值，ASM_前为待定影像均衡化处理前的二阶矩，ASM_后为待定影像均衡化处理后的二阶矩；

步骤16，重复步骤14、步骤15，直至所有的待定影像都被标记为止。

利用灰度共生矩阵统计纹理特征是一种常见的方式，本发明为了突出遥感影像上的薄云，对待定影像进行直方图均衡化处理，突出影像上的云域亮度，减弱非云区域亮度。灰度共生矩阵前后变化，分别统计了均值(Mean)和二阶矩(ASM)，如表2所示：

表2含云影像灰度共生矩阵变化

从Mean角度来说，沙地由于较粗的纹理特征，均值逐步递减，而云逐渐增大，利用均衡化前后的均值做比值就可以实现区分。如式(1)：

其中，“前”、“后”表示影像均衡化前后的均值。

同时，雪在均值和二阶矩(ASM)方面前后差异较小，但从均衡化后，雪的二阶矩低于云，采用如式(2)所示的差分处理可实现云雪易混淆的地物之间的区分。

-0.03≤ASM_前-ASM_后≤0.03 式(2)

同时满足式(1)和式(2)的划定为含云影像，否则为无云。

步骤2，确定含云影像上的云域范围及该云域范围的含云量，并划定含云影像质量等级；

其中，检测含云影像上的云域范围及该云域范围的含云量，去除含云量超过含云阈值的含云影像是指：

步骤21，通过建立决策树云检测规则依次检测步骤2中获取的所有含云影像中的云域范围，所述的决策树检测规则如下：

其中，a为云在红波段的灰度异常值，x为灰度值；

其中，云在红波段的灰度异常值a的计算方法为：

步骤211，对含云影像的灰度直方图中的频率做一阶向前差分；

步骤212，设定阈值L，若前后频率差分值＞L时，表明该频率段不在缓冲区；否则，该频率段为疑似缓冲区；

取最长的一段疑似缓冲区为频率缓冲区；

步骤213，通过冒泡法得到含云影像的灰度直方图上的可能异常值，计算频率缓冲区的频率均值w，若该可能异常值大于2w，则该可能异常值即为云在红波段的灰度异常值a。

本实施例中最长的缓冲区即为图2中的平滑部分(约为350-930)，求出这段的频率均值为w。通过冒泡法得到曲线上的峰值异常值可能为103和938，只有938的频率与w的差值大于2w，且938后的灰度与均值w没有交点，因此938即为异常值。

步骤22，分别统计含云影像和云域范围的像素个数，通过下式计算云域范围的含云量：

其中，C为含云量，N_c为云域范围的像素个数，N_I为含云影像的像素个数；

步骤23，设定含云阈值D，本实施例中D取50％,若含云影像中的含云量C大于D时，去除该含云影像。

如图4、图5所示分别是原始含云影像和检测后的云域范围，其中，图5中白色的为云斑，黑色的为非云；

本实施例采用混淆矩阵可以快速评判检测的精度，如表3所示。

表3云检测后精度评价

从表3可知，Kappa系数均在0.93以上，说明具有较好的一致性，再次验证了本方案的精度。

主要根据含云量C的大小划定质量等级，主要包括优秀(C≤0％)、良好(0％<C≤30％)、及格(30％<C≤50％)、废片(C>50％)四等。对资源三号卫星影像检测了53轨共计4233景影像，检测出含云影像579景影像，影像绝大部分处于良好及及格级别，说明总体成像效果较好，大部分可以直接进行信息提取。对于含云影像，则需经过去云才可再次使用。更新影像数据库，得到表4。

表4更新后的影像数据库

步骤3，去除含云影像上的云，并得到无云影像；

步骤31，依次选取划分为优秀、良好、及格等级的含云影像作为含云影像A，从批量遥感影像数据库中选取与该含云影像A的经纬度相同但不同时相的所有无云影像,从所述无云影像中选择质量好且分辨率与该含云影像A相同的无云影像B；

步骤32，分别裁剪出所述含云影像A与无云影像B的云域区域C1、C2，用C2区域替换掉C1区域，镶嵌得到中间影像A1；

由于采用的影像都是具备地理坐标的，经过决策树云检测生成的矢量范围C同样具备地理坐标，因此裁剪出的C1和C2具备相同的地理位置，不必进行二次配准；

步骤33，对中间影像A1采用均值滤波处理得到新的遥感影像A2，即无云影像；

通过步骤32得到的中间影像A1不能直接使用，会出现镶嵌后明显的拼接线，因此需要对A1拼接线处的灰度值进行处理。本发明采用3×3窗口搜索拼接缝中心线处的像素周围的9个栅格，如(x,y)处灰度值为g(x,y),取这9个栅格的像素灰度值的均值g′(x,y)为其灰度值。同时以中心线像素为准，分别往内外各缓冲5个像素，依次采用相同的方式进行采样，取均值作为新的像素灰度值，得到新的遥感影像A2，更新遥感影像数据库。其中，中间边界是拼接线，中间边界内外的两个边界是以中间边界向内向外5个像素进行缓冲得到的缓冲区，表格中的322即是均值滤波后的像素灰度值,A1是替补后的影像，A2是均值滤波处理的图像。

步骤4，重建批量遥感影像数据库。

通过步骤3去云处理之后得到影像质量相对较好的影像，按照表4所示对影像数据库进行重新编目，对所有影像重新划分质量等级,同时更新遥感影像数据库。

本实施例中采用含有浅色矿物的遥感影像进行去云处理，遥感影像上的浅色矿物，颜色浅，密度小。含二氧化硅和三氧化二铝较高，如长石、白云母、石英等。影像上的云易与浅色矿物造成混淆，提取难度加大，而单景处理根本不适用于大批量遥感影像的生产，因此本发明先对影像上云进行批量处理，大幅度提高浅色矿物的提取精度。

从遥感影像数据库中提取浅色矿物，以硅化蚀变信息为例(主要是二氧化硅等氧化物)，首先利用光谱仪对野外采集到的矿物测定反射光谱曲线。采用波段比值运算(如其中Band₄,band₃，band₂分别是近红外波段、红波段和绿波段)，求出各变量相关矩阵系数，得到如式(5)所示的线性回归方程，即可提取出遥感影像上的浅色矿物。

其中，x，y分别是光谱特征向量和浅色矿物中二氧化硅的变量。