一种GF‑4卫星序列图像云与云下阴影检测方法与流程

本发明涉及遥感图像辐射处理技术，具体的说，涉及一种针对高分辨率遥感卫星序列图像的云与云下阴影检测技术。

背景技术：

遥感图像的辐射预处理一直是遥感数据处理的主要课题之一。新的遥感卫星投入使用后，对新数据的预处理便是紧要的问题。遥感图像的预处理通常包含几何预处理于辐射预处理，其中的辐射预处理除了关键的辐射定标外，对图像云量的统计也是一项重要的步骤。根据卫星数据用户方的使用需求，云量所占比例成为选择卫星图像的一项重要指标，比如在分类与制图应用中多选用云量尽量少的图像，而在气象与减灾应用中更关注云量多的图像。目前绝大多数遥感卫星图像数据产品都包含云量信息，很多都包含专门的云覆盖标记波段，方便用户逐像素区分云与地表。遥感卫星图像的云检测方法很多，根据卫星数据的特点开发适用的云检测算法是数据预处理的重要步骤，比如波段范围以可见光与近红外为主的高分辨率多光谱遥感图像云检测多采用简单的直方图统计结合自动阈值的方法，或者基于无云地表参考数据的自动阈值方法，这类方法对于雪地与高亮干燥地表容易造成误检；包含红外波段的、定量化程度高的观测卫星数据多依赖云层的低温特性检测云，一般精度较高，数据产品中常包含专门的云标记波段。

高分四号卫星(以下简称gf-4)是中国于2015年12月发射的一颗地球同步轨道卫星，搭载空间分辨率为50米的全色、多光谱相机和400米分辨率的中波红外相机，采用面阵凝视方式成像，成像间隔快至20秒，具备高时间、高空间分辨率的优势。自2016年2月3号国防科技工业局公布首批影像以来，gf-4已经获取了中国及周边区域大量数据，在检测森林火灾、洪涝灾害等方面发挥着重要作用。gf-4卫星图像的预处理同样包括几何与辐射两个部分，几何预处理包括系统成像模型的构建、控制点匹配与几何精纠正等，目标是实现同一地图投影下成像数据的逐像素配准。辐射处理包括辐射定标、云与云下阴影的检测等，目标是使得图像像素值能精确地描述地表辐射状况。

gf-4卫星数据预处理技术研发在充分利用之前卫星数据处理成果的基础上，还要考虑gf-4卫星图像本身的特性，研发专门的处理算法。根据gf-4卫星数据的特点，云检测实现的技术途径主要有两个：一是利用云的低温特性，根据云在中波红外波段亮度低与在可见光波段亮度高的特性检测，这是目前采用较多的方法；二是利用云的运动特性从序列图像中检测云，由于gf-4卫星采用地球同步轨道，并面阵凝视成像方式，容易获取同一地理区域下大量的图像，根据云的运动特性，利用序列图像便可区分出云与云下阴影。

通过对gf-4具体数据的分析与实验，发现中波红外数据难以用于精细的云检测，原因主要有：分辨率差异过大，一个400米分辨率的中波红外像素对应一个8×8像素块、共64个50米分辨率的可见光像素；覆盖地理范围不重合且像素尺寸不同，中波红外为1204×1024像素，可见光为10240×10240像素；成像时间存在差异，中波红外与可见光图像被设计为成像时间固定间隔45秒，云的运动特性便使得两者数据中云的位置、形态存在差异；数据质量上，中波红外相机成像还存在明显的噪声问题。

针对gf-4卫星图像云检测的问题，在中波红外图像应用困难的情况下，研究利用序列图像提高云检测的精度成为一种有效的技术途径，gf-4在实际的灾害监测应用中，对灾害区域的多天定时连续观测数据通常具有较强的序列特性，这为基于多幅/序列图像的云检测提供了有利的条件。工程化生产gf-4卫星数据辐射预处理产品需要一种稳健的序列图像云与云下阴影检测算法。

技术实现要素：

本发明的目的是针对gf-4静止卫星图像预处理应用，提供一种序列图像云与云下阴影检测技术，特别是对于gf-4卫星50米空间分辨率的全色、多光谱图像的l1级数据产品中的云检测产品生产，提供一种生产标记云与云下阴影波段产品的算法流程。本技术基于成熟的遥感图像相对辐射校正算法及s-g(savitzky-golay)滤波算法，根据gf-4卫星图像的辐射预处理需求与序列图像辐射特性定制的快速云与云下阴影检测算法流程。

本发明的基本思路为：对于gf-4静止卫星获取的同一地理区域序列图像数据，对在没有系统几何模型的情况下，利用图像自动匹配获取的线性函数实现图像之间相对位置关系的逐像素配准，利用自动相对辐射校正减小序列图像之间不同成像时间导致的辐射差异，利用s-g滤波结合自动阈值在序列图像中逐像素修正地表，通过比较像素修正前后的值划分出云与云下阴影，最终输出结果是序列中每幅图像都对应一个单波段云与云下阴影标记数据。

所述的gf-4卫星序列图像，限定为50米空间分辨率的全色、多光谱图像，图像四个角点的经纬度差异不超过±0.3度，可以是gf-4卫星凝视模式下获取的序列图像，也可以是不同天获取的，按照获取时间先后排列的序列图像。

本发明的技术方案提供的gf-4卫星序列图像云与云下阴影检测方法，其特征在于包括以下实施步骤：

a数据预处理，获取序列图像相对配准线性参数；

b线性相对辐射归一，自动提取序列图像两两之间伪不变特征地物点，通过统计比较各图像伪不变特征地物点辐射差异，找出整体辐射差异大的数据进行相对辐射校正；

c序列图像云检测，根据序列图像的数量，选择基于s-g滤波的算法，或者基于统计的自动阈值法标记云与云下阴影，获取每个图像的云与云下阴影掩膜波段数据；

d修正检测结果，先根据检测的阴影像元与最近的云像元的距离修正阴影像元，再通过顺序叠加显示原图与云检测结果图找出检测结果存在问题的数据进行云区修正，如果存在多幅云检测结果差的数据，则将这些数据组成新的序列图像重新进行上述云检测处理。

上述实施步骤的特征在于：

步骤a中所述数据预处理，包括对输入的序列图像进行数据完整性检查、地理覆盖范围检查、按地表平均辐射亮度排序、以及程序运行的一些准备初始化处理。所述获取序列图像相对配准线性参数，具体过程是对同一地理区域的序列图像，利用图像四个角点的近似经纬度坐标确定图像间的大致相对位置关系，并根据定位误差确定图像分块自动匹配的检索范围，通过自动匹配，获得控制点数据并拟合出一个线性函数，序列图像之间通过不同的线性函数实现相对配准，而不对图像本身进行插值与重采样等变换。

步骤b中所述的线性相对辐射归一，用于减小序列图像两两之间因数据获取时间的不同导致的辐射差异，序列中相邻两幅图像之间采用遥感领域中的相对辐射校正技术，这里采用基于多源典型相关分析的ir-mad(re-weightedmultivariatealterationdetectiontransformation)变化实现自动提取两图像的伪不变特征地物点；所述通过统计比较各图像伪不变特征地物点辐射差异，采用的比较方法是对于序列中的一幅图像，其与序列前后的图像提取的伪不变特征地物点之间的均值差异在整个序列中的大小比较，如果该图与序列中前后图像之间的伪不变地物点均值差异都很大，则需要对该图像进行相对辐射校正；所述找出整体辐射差异大的数据进行相对辐射校正，相对辐射校正采用伪不变地物点拟合出的线性函数完成。

步骤c中所述序列图像云检测，根据序列图像的数量采用不同的检测方法，当图像数量大于等于10，选择基于s-g滤波的算法，先对序列图像进行逐像素的序列滤波，再根据每个图像滤波前后像素值的比较区分云与云下阴影；当图像数量小于10，则直接统计序列图像逐像素的均值与中值，根据每个图像像素值与统计的均值与中值的差异区域云与云下阴影；该步骤获取每个图像的云与云下阴影掩膜波段数据；

步骤d中所述修正检测结果，先根据检测的阴影像元与最近的云像元的距离修正阴影像元，根据云下阴影距离云的最大距离，比如gf-4图像可取值500或者1000素数，对于每个检测到的阴影像元，如果在该半径像素范围内没有找到检测的云像元，则判定为误检，将该像元从阴影像元中剔除；通过顺序叠加显示原图与云检测结果图找出检测结果存在问题的数据进行云区修正，这里所述存在问题的数据，主要针对两种情况，一是检测结果存在较多碎片与漏洞，此时采用计算机形态学的办法去除碎片并填充漏洞，二是检测结果存在大比例的明显误检，比如将高亮地表、水面等误检为云，此时需要对改图进行重新检测，如果存在多幅云检测结果差的数据，则将这些数据组成新的序列图像重新进行云检测处理。

本发明与现有技术相比有如下特点：本发明提供了一种针对gf-4序列图像云与云下阴影检测解决方案，通过线性函数实现gf-4序列图像逐像素位置关系的快速配准，利用线性相对辐射校正减小数据间的辐射差异，利用s-g滤波结果划分出云与云下阴影。算法自动化程度高，云与云下阴影检测过程无需人机交互，用户仅需对最终的检测结果进行简单的检查，对个别数据进行重新处理。涉及的关键步骤采用成熟的算法实现，具有较高的稳定性与适用性。对于gf-4卫星数据预处理中的云与云下阴影检测产品的生产与产品精度的提升，提供了关键的技术支撑。

附图说明：

图1是gf-4卫星序列图像云检测流程图

图2是单个像素处s-g滤波与阈值云检测的示意图

图3是检测结果形态学修正示意图

具体实施方式：

本技术的思想是利用云与云下阴影在序列图像中的运动特性实现云与云下阴影的检测，其必要条件为：gf-4卫星获取的数据易构成序列，且在不同时间获取的同一区域的多幅图像之间仅存在旋转与平移的位置关系。该必要条件对于gf-4卫星图像是满足的，是由于gf-4卫星采用静止卫星轨道，相对于地球的位置是固定的，其成像几何是不变的，包括对地球可观测范围内的任意一点到卫星传感器成像点的几何关系都是固定的。静止卫星位置的固定性，保证了对同一区域的多个观测图像，在图像中心点与四个角点坐标相同的情况下，所有图像的系统成像几何模型是相同的，使得图像的系统畸变与空间分辨率都是一致的。并且gf-4卫星采用面阵凝视方式成像，成像时间是瞬时的，图像的成像过程不会引入新的几何畸变，卫星的抖动与传感器的抖动都很难对成像造成影响。最后，gf-4卫星成像定位精度高。通过指向控制，gf-4可实现对中国及周边地区进行自由观测，也可采用凝视模式对幅宽400千米的固定区域进行持续观测，其定位精度达±0.1度，使得gf-4卫星具有获得同一固定区域的序列图像的能力。自2016年2月3号国防科技工业局公布首批影像以来，gf-4卫星已经获取了中国及周边区域大量数据，其中包含大量的可构成序列的数据，包含gf-4卫星凝视模式下可获取同一天相近时刻的序列图像，也有不同天获取的同一区域图像构成序列。

采用本发明实现gf-4卫星序列图像云与云下阴影检测流程如图1所示，现结合附图对其进行描述。

处理单元111数据预处理，数据预处理针对gf-4卫星序列图像，对于gf-4数据发布的50米空间分辨率的可见光与近红外图像，有两种方式可以获得序列的图像：一是gf-4相机在凝视工作模式下获取的同一天、相近时间内获取的连续多幅图像组成的序列；二是非凝视模式下，不同天获取的、同一地区的多幅图像。在gf-4卫星的实际运营中，获取第二种非凝视模式构成的序列图像数据很容易。因为灾害监控作为gf-4的重要任务，当某一区域发生灾害时，gf-4可以在经纬度误差0.1度的范围内对固定区域进行多天、多次观测，从而构成序列图像。数据预处理算法程序对输入的序列图像进行数据完整性检查、地理覆盖范围检查、序列排序、以及程序运行的一些准备初始化处理。其中序列排序并不是按照数据获取的时间进行排序，而是根据各图像地表的平均辐射亮度排序。

根据各图像地表的平均辐射亮度排序在本发明中，对于gf-4静止卫星数据的后续处理是关键的。不同于太阳同步轨道卫星获取地球同一地理位置图像的时刻是相近的，地球同步轨道卫星相对地球位置不变，可选择一天中的任意时刻成像，白天里不同时刻成像根据太阳高度的不同，图像整体辐射亮度也是不同的，比如早上8点钟与正午12点钟成像的图像辐射差异。具体排序方法是统计序列图像中蓝光波段的直方图，排除过亮与过暗像素值后，将剩余像素的均值作为排序的依据。这里排除过亮与过暗像素值的目的是过滤掉可能的云与云下阴影，保证得到的均值尽量代表地表的辐射亮度情况。

序列图像排序后，由于本发明属于gf-4数据辐射预处理，采用的图像是未经系统几何校正的原始数据，仅根据图像定位信息直接叠加无法实现地物的逐像素配准。这里采用线性相对配准方法，具体过程是对同一地理区域的序列图像，利用图像四个角点的近似经纬度坐标确定图像间的大致相对位置关系，并根据定位误差确定图像分块自动匹配的检索范围，通过自动匹配，获得控制点数据并拟合出一个线性函数y＝ax+b，序列图像之间通过不同的线性函数参数a与b实现相对配准，而不对图像本身进行插值与重采样等变换。

通过预处理的序列图像数据，如果是非凝视模式获取的不同天数据，则需要进行处理单元112的线性相对辐射归一处理。该处理是算法流程的关键步骤，用于减小序列图像两两之间因数据获取时间的不同导致的辐射差异，序列中相邻两幅图像之间采用遥感领域中的相对辐射校正技术，对于序列中的一幅图像，其与序列前后的图像提取的伪不变特征地物点之间的均值差异在整个序列中的大小比较，如果该图与序列中前后图像之间的伪不变地物点均值差异都很大，则需要对该图像进行相对辐射校正。这里采用ir-mad变化实现自动提取两图像的伪不变特征地物点。

ir-mad变换源于nielsenetal.(1998)提出的mad变换，该算法为了遮蔽两时相图像中的变化像素，首先形成两幅图像n个通道内像素值的线性组合。用随机向量x和y分别表示目标图与参考图重叠区内筛选出的像素值。根据以下变换公式：

u＝a^tx＝a1x1+a2x2+λ+anxn

v＝b^ty＝b1y1+b2y2+λ+bnyn

其中ai与bi为mad系数，mad变换最小化u与v之间的正相关。在服从约束：var(u)＝var(v)＝1的前提下，定义mad变量：

mad＝var(u-v)＝var(u)+var(v)-2cov(u，v)＝2(1-corr(u，v))→maximum

最小化正相关系数corr(u，v)是一个标准的统计过程，即所谓的广义特征值问题。求出的mad变量各个分量相互正交，并且是线性变换的不变量。本发明之所以选择mad变换来提取不变特征点，正是由于mad变换的这种对变量x与y之间的线性关系不敏感的特性，可以很好地适应不同时间获取的gf-4图像之间存在的较大辐射差异。ir-mad变换进一步提高了mad算法的精度与稳定性。

利用ir-mad变化自动提取出的两图像的伪不变特征地物点，采用最小二乘方法拟合出一个整体的线性函数y＝ax+b，利用传统的遥感图像线性相对辐射校正将序列中辐射差异大的图像校正到相邻一幅图像的辐射水平。

根据序列图像包含的图像数目是否大于10幅，确定后续处理采用处理单元113，或者处理单元114。

处理单元113统计与自动阈值云检测，具体算法是，假若序列图像包含n幅图像，n≤10，对序列图像中的每一个相对配准的像素位置，统计n个像素的均值vmin与中值vmid，假若均值与中值数值上相差不大，比如|vmin-vmid|＜10，则n个像素全部标记为地表。假若中值与均值数值差异大，则将n个像素逐次与vmid进行比较，如果vi-vmid＞vcloud，则判定第i个像素为云，vcloud为云阈值；如果vi-vmid＜vshadow，则判定第i个像素为云下阴影，vshadow为云下阴影阈值，为负值；vcloud与vshadow的可取值±2|vmin-vmid|，或者±3|vmin-vmid|。

处理单元114s-g滤波与阈值云检测，该步处理是先对序列图像进行s-g滤波，再通过比较滤波前后数值的变化差异判定是否为云与云下阴影。

s-g滤波通过滑动窗口多项式拟合来达到对序列数据进行平滑的目的(savitzky&golay，1964)。序列数为n，对其中长度为n＝2m+1的子序列进行k(k≤n)阶多项式拟合可表示为：

s-g滤波过程是对序列中的某一点t0及其左右m邻域共n＝2m+1个点(ti，yi)，i∈[-m，m]，进行k阶(k≤n)的多项式拟合，用拟合后的滑动窗口中心的数据(t0，y0)置换原始时间序列中的数据(t0，y0)，然后向右移动窗口，使窗口中心移至序列中下一数据，重复上述过程，直到滑动窗口到达序列末尾。平滑窗口系数是通过最小二乘法方式求得。

假若序列图像包含n幅图像，n＞10，对序列图像中的每一个相对配准的像素位置，经过s-g滤波后，将n个像素vi逐次与滤波后的值vi-sg进行比较，如果vi-vi-sg＞vcloud，则判定第i个像素为云；如果vi-vi-sg＜vshadow，则判定第i个像素为云下阴影；vcloud与vshadow的可取经验值，比如20或者30，根据具体数据情况可调整阈值。序列图像中单个像素处s-g滤波与阈值云检测的示意图见图2。

对云下阴影的检测结果进行修正，根据检测的阴影像元与最近的云像元的距离修正，给出云下阴影距离云的最大距离，比如gf-4图像可取值500或者1000素数，对于每个检测到的阴影像元，如果在该半径像素范围内没有找到检测的云像元，则判定为误检，将该像元从阴影像元中剔除。

检测结果的精度建议进行人工检查，对于精度较差的结果可根据具体数据情况通过调整阈值后进行再次处理。另外由于云检测是基于单个像素的，有时会出现检测结果局部区域出现大量碎片与漏洞的情况。为了改善云检测的边界效果，采用计算机形态学的方法进行云区整饰处理，去除云边界外小于一定像素数的孤立云区，填充云边界内小于一定像素数的空洞，然后整饰云边界。效果示意图见图3。由于实际的云本身也可能是过分离散的，所以是否对云区进行整饰由用户决定。

云与云下阴影检测的结果保存为8位单波段图像，地表取值0，云下阴影取值1，云取值2。n幅序列图像对应n幅检测结果，作为gf-4初级数据产品提供给用户。

本发明的实例已经在pc平台上实现，目前已经交付用户方进行测试与使用，作为gf-4数据辐射预处理中云特征参数反演关键技术。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。