历史平均值法用于MODIS影像像元云补偿——以甘肃省为例

陈宝林，张斌才，吴静，李纯斌，常秀红

(1.甘肃农业大学资源与环境学院，兰州 730070；2.甘肃省基础地理信息中心，兰州 730030)

0 引言

目前在遥感技术领域，除微波传感器能部分穿透云层获取地表信息外，其他可见光传感器均未能彻底解决影像成像时会受到云或雾的干扰的问题，遥感影像往往存在信息缺失的现象，造成了遥感影像数据在空间上的不连续[1]。遥感卫星影像的成像一直都存在云覆盖及雾遮挡等问题。在遥感卫星获得的遥感影像数据中，有很大一部分是光学影像及红外影像，极易受到气候因素的影响，而云层遮挡就是其中最常见的影响之一。

云雾的存在作为不可抗拒的外在因素严重影响了卫星遥感影像成像的质量、影像的识别以及影像的应用。这给遥感影像的使用带来了极大的不便。云和雾在时间和空间分布上存在着极大的不确定性，其高度、厚度、组成种类以及太阳高度角、方位角变化导致卫星云图特征千变万化，始终是卫星遥感图像处理与应用中的一大难题[2]。

国内外很多学者利用多光谱通道资料对卫星遥感图像上的各种云/云层进行了分类研究，早期的研究皆以云雾的光谱特征为主[3]，随着计算机视觉、图像识别技术的日趋成熟，纹理特征被逐步引入到遥感影像中的云雾识别等处理中，并逐渐发挥了重要的作用，不但丰富了计算机图像处理的理论，而且提高了云雾识别与处理的效率和准确性。但要想形成一个系统、全面、成熟的云雾成像机理还需要更进一步的研究。尤其是云的存在使得遥感影像在成像过程中存在着不同程度的阴影；针对遥感影像存在阴影这一问题，很多学者提出了针对阴影的补偿方法；而现有的阴影补偿方法在普遍性及适用性上都存在着一定的局限性，因为缺少光谱信息，只有相应的R，G，B波段信息[4]，在数据预处理—镶嵌的时候会出现颜色偏差、像元0值等其他情况，如果需要产生新的地理地图，那这样的数据结果必然难以符合要求，故而需要进一步的后期处理，包括颜色均衡化，像元0值替换或像元0值补偿等其他处理方式。

从理论上讲，要想完全去除遥感影像上的阴影其可能性是微乎其微的，但遥感影像上阴影区的色彩和亮度信息为遥感影像的阴影进行补偿提供了必要的可能条件[5]。目前存在的补偿方法包括线性相关修正法、指数校正法[3]、直方图匹配法，同态滤波法[6]、归一化处理法等[7]，其不同程度地都对遥感影像阴影区域的图像增强有一定的作用，但到目前为止还没有一种被学术界公认的方法[8-9]。遥感影像阴影补偿的关键在于既要提升阴影区域的视觉效果,又不能改变非阴影区域的信息[10]。

本文以甘肃省为研究区，对研究区2017年MODIS11A1数据进行数据缺失检测，并提出基于物候节气，为时间段的历史平均值法对缺失的数据进行补偿。借助Python及与之相关的工具试图寻找新的具有可操作性和可重复性的的遥感影像云补偿方法。基于此需要解决补偿数据的选取、处理、拼接及适用于本思路的时段划分；影像精度设定、代码条件设定、代码的书写、运行和运行结果的判定；影像像元有效率的计算及统计，然后运用历史平均值法进行补偿，将遥感影像像元云补偿后的结果与补偿前的遥感影像进行对比，得出影像云补偿的精度及优点，并试图探索此方法是否也同样适用于其他卫星的数据类型。

1 研究区概况及数据源

1.1 研究区概况

本文选取甘肃省2017年MODIS11A1遥感影像数据，将历史平均值法用于MODIS影像像元云补偿。甘肃省地处黄河上游，介于E92°13′～108°46′，N32°31′～42°57′之间，东接陕西，东北与宁夏毗邻，南邻四川，西连青海、新疆，北靠内蒙古，并与蒙古国接壤，总面积4.54×105km2 [11]。甘肃省地形地貌复杂多变，地势自西南向东北倾斜，地形呈狭长状[12]；处于黄土高原、青藏高原以及内蒙古高原3大高原的过渡交汇地带；遥感影像受到云层的遮挡较为明显[13]，因此选取甘肃省来进行研究具有较强的代表性。

1.2 数据源及其预处理

本文采用的2017年MOD11A1数据是地表温度产品，来自美国国家航空航天局网站。借助于遥感土壤水分监测系统软件，以JAVA为承载展开处理，通过书写代码的方式对参与预处理的遥感影像数据进行数据的质量控制，逐月处理完2017年全年12个月的原始数据，这一过程实现了MOD11A1的2017年的原始遥感影像数据的拼接和处理。

然后运用Python书写代码对2017年的MOD11A1像元值为0的数值出现的频率(以下简称“频次”)进行提取，统计像元0值出现的频次是MODIS影像像元云补偿的基础，为接下来的影像云补偿提供了重要的依据。卫星在高空扫描地物信息的过程中因受到低层暖湿空气被迫抬升到高空冷却凝结形成的“云”的遮盖而表现为像素的像元值为0；并对2017年的MOD11A1的1 km分辨率的每日地表温度数值进行提取。将运算结果用HDFVIEW打开，在打开的界面随即出现了两个遥感影像像元值为0值的二维频次矩阵(包括一个白天像元矩阵(Day—temp-numb)和一个夜晚像元矩阵(Night—temp-numb))，如图1所示。

图1 MOD11A1数据2017年白天像元0值频次分布矩阵Fig.1 MOD11A1 data 2017 daytime pixel 0 value frequency distribution matrix

在此基础上对出现在2 580×3 080大小的一个二维矩阵中的像元值为0值的频次进行统计，统计采用分段统计、确定组距的方法进行。首先将白天的像元矩阵在HDFVIEW软件中打开，然后将其数据导入EXCEL，找出最大值和最小值，利用最大值与最小值的差大致确定如何分组；然后确定组距和组数。经统计发现白天像元矩阵(Day—temp-numb)中白天像元值出现的的最大值是313，最小值是83；夜晚像元矩阵(Night—temp-numb)中夜晚像元值出现的的最大值是363，最小值是67。本文采用的是2017年全年的数据，共有365 d，因此需要进行补偿的像元值天数应当大于等于0 d，且同时满足小于等于365 d。采用分组分段统计同一个像元在一年中有多少天像元值为0值，以期为补偿提供重要的依据。

2 研究方法

对2017年的MODIS11A1的数据进行统计，发现需要补偿的卫星遥感数据大约占到了全年的1/3～1/2；以此为切入点试图提出新的针对遥感影像像元云补偿的方法。运用Python对2017年MOD11A1的数据进行第二次处理，运用Python中的GDAL等其他工具并编写代码将第一次预处理的数据加入到Python 3.8.0 shell软件中对其进行运行；因该方法是采用倒推天数的办法对遥感影像像元云补偿，因此被称作为历史平均值法。

采用历史平均值法在补偿时间的选取上考虑引用物候中的节气。物候是指动植物与当地的生态环境协同进化而形成的生长发育节律现象[14]。依据不同动植物的生长、发育和活动的变化节律进行生产活动的时间制度称为“物候历”。在二十四节气发明之前，我们的远古祖先最初使用的是“物候历”二十四节气确立后，成为我国最早的结合天文、气象、物候知识指导农事活动的历法。据公元前2世纪的《逸周书·时训解》记载，一年二十四节气，共七十二候[15]。

本文采用的是MOD11A1的2017年的遥感影像数据，将下载的数据进行拼接之后，采用的是向前倒推求历史平均值的方法，文中随机选取2017年的7月22日(大暑)、8月7日(立秋)、12月7日(大雪)、12月22日(冬至)的白天的MODIS11A1的甘肃省的影像；对选取的4 d的遥感影像进行补偿，具体的补偿日期和补偿倒推日期(对应节气前15天)对应见表1。

表1 补偿日期和倒推日期对照表Tab.1 Comparison table of compensation date and reverse date

将原始影像导入ArcGIS并嵌套甘肃省轮廓图进行掩模处理，将裁剪后的甘肃省影像图导出。图例中的像元值代表了图中影像点的地表温度，像元值和地表温度对应关系如下：

T=0.02K-273.15，

(1)

式中：K为该点的遥感影像像元值；T为该点像元值对应的地表温度，℃；0.02是系数常数。运用像元值与地表温度关系公式结合所选取的4 d的遥感影像像元值计算出当天地表温度变化范围为：

1)7月22日对应第203天，白天(-8.29～62.61 ℃)；

2)8月7日对应第219天，白天(-6.57～63.63 ℃)；

3)12月7日对应第341天，白天(-22.37～17.87 ℃)；

4)12月22日对应第356天，白天(-25.57～18.73 ℃)。

甘肃省山地和高原约占甘肃省总土地面积的70%，多山的地形加剧了下垫面的不均一性，也决定了其夏季和冬季的遥感影像像元成像更容易受到云的影响。因此，随机选取了夏季的大暑和立秋两个节气，冬季的大雪和冬至两个节气。

本文将随机选取的大暑(7月22日)、立秋(8月7日)、大雪(12月7日)、冬至(12月22日)4个节气当天及各自之前15 d白天的遥感影像像元在拼接之后嵌套甘肃省的轮廓图进行掩模提取后的有效像元利用率进行了统计。具体方法是统计像元栅格面积大小，打开像元的属性表统计弹出对话框的Count和Sum字段，将其统计作以记录；然后打开导入的甘肃省底图属性表，添加Area字段之后点击几何计算；最后求得其面积。计算面积时单位选择km2；如此往复完成所选择的所有天数的数据统计。之后进行计算对应日期的有效像元利用率。Sum字段的计算是字段总和乘以分辨率，即Sum×0.5×0.5(分辨率是500 m×500 m)，这是为了下一步计算时不牵扯单位，因为该步骤已经将单位换算为统一单位了；然后用Sum字段乘以分辨率的数值去除以甘肃省底图的面积最后乘以100%，即为有效像元利用率。

3 结果与分析

3.1 MOD11A1数据0值统计结果

对MOD11A1数据2017年的白天和夜晚的像元0值采用划定区间，分段统计的方法进行统计并绘制频率分布直方图，结果见图2和图3。图中，横坐标为以元旦为1算起的天数，通常称为儒略日。

图2 2017年白天像元0值频次分布图Fig.2 Frequency distribution of 0 value of daytime pixels in 2017

图3 2017年夜晚像元0值频次分布图Fig.3 Frequency distribution of 0 value of night pixels in 2017

对由若干景影像组成的2 580×3 080像元这样大的一个遥感影像覆盖区域进行预处理，初步统计得到了2017年白天像元0值出现频次最多的区间是[123,133)，[133,143)和[143,153)；3个区间的频率总和为36.34%，白天像元0值出现频次总和超过200 000次以上的数值主要集中在[113,223)这个大区间上，这一数据表明了2017年有[113,223)天白天的遥感影像数据因为云的遮盖而表现为0值，缺失程度占到了2017年全年的30.96%～61.10%。这一统计结果再次证明了2017年MOD11A1白天的遥感影像因受到气象要素云的影响；从而对卫星过境时扫描的成像效果产生了较大的影响。2017年全年共有365 d，统计结果表明2017年全年有1/3～1/2的遥感影像因为存在云的遮盖而影响到了遥感影像像元的有效利用。

2017年夜晚像元0值出现频次最多的区间是[97,107)，[107,117)和[117,127)；3个区间的频率总和为36.25%，夜晚像元0值的出现频次超过200 000次以上的数值主要集中在[87，197)这个大区间上，这一数据表明了2017年有[87，197)天夜晚的遥感影像数据因为云的遮盖而表现为0值，缺失程度占到了2017年全年的23.83%～53.97%。这一统计结果证明了2017年MOD11A1夜晚的遥感影像也受到了气象要素云的影响。

将白天和夜晚的两幅像元0值频次分布直方图的每一个分段区间的最高点分别用平滑曲线连接起来可以得到两幅折线图；从折线图的大体走向可以看出2017年白天和夜晚的像元0值频次折线几乎具有相似的变化走向，大体上遵循先急剧上升后缓慢下降的总趋势。对2017年的MOD11A1遥感影像像元0值频次统计之后，又将2017 年按4个季度进行了4次统计。分别统计了2017年4个季度的MOD11A1的遥感影像的白天像元0值出现的频次，以下的统计中将每个季度缺失的天数按照处理后的结果采用等距分组，结果见表2—表5。频数表示在这一季度矩阵中的所有像元缺失天数出现的次数总和。频次占比表示对应的分组的频数在频数总和中所占的百分比。其结果基本与2017年全年白天的结果保持一致。

表2 2017年第一季度白天像元0值频次表Tab.2 Frequency table of daytime pixel 0 value in the first quarter of 2017

表3 2017年第二季度白天像元0值频次表Tab.3 Frequency table of daytime pixel 0 value in the second quarter of 2017

表4 2017年第三季度白天像元0值频次表Tab.4 Frequency table of day pixel 0 value in the third quarter of 2017

表5 2017年第四季度白天像元0值频次表Tab.5 Frequency table of daytime pixel 0 value in the fourth quarter of 2017

3.2 数据补偿结果

本文选择2017年大暑、立秋、大雪以及冬至4个节气，将白天的遥感影像与甘肃省轮廓底图进行掩模处理，得到未补偿的有效遥感影像数据，结果如图4—图7所示。

图4 2017年7月22日白天MOD11A1影像Fig.4 MOD11A1 images during the day on July 22,2017

图5 2017年8月7日白天MOD11A1影像Fig.5 MOD11A1 images during the day on August 7,2017

图6 2017年12月7日白天MOD11A1影像Fig.6 MOD11A1 images during the day on December 7,2017

图7 2017年12月22日白天MOD11A1影像Fig.7 MOD11A1 images during the day on December 22,2017

结果显示没有补偿的7月22日白天(2017年第203 d)的遥感影像像元有效利用率是21.50%；8月7日白天(2017年第219 d)的遥感影像像元有效利用率是37.10%；12月7日白天(2017年第341 d)的遥感影像像元有效利用率是66.65%；12月22日白天(2017年第356 d)的遥感影像像元有效利用率是57.34%。统计所选的4个节气中的大暑及大暑前15 d、立秋及立秋前15 d、大雪及大雪前15 d、冬至及冬至前15 d(共64 d)白天的遥感影像，64 d的像元平均有效利用率是58.42%；有效像元利用率统计见表6—表9。

表6 甘肃省2017年大暑及大暑前15 d白天像元有效率Tab.6 The daytime pixel efficiency of the 2017 Great heat and 15 days before the Heavy snow in Gansu Province

表7 甘肃省2017年立秋及立秋前15 d白天像元有效率Tab.7 The daytime pixel efficiency of the 2017 autumn begins and 15 days before the autumn begins in Gansu Province

表8 甘肃省2017年大雪及大雪前15 d白天像元有效率Tab.8 The daytime pixel efficiency of the 2017 Heavy snow and 15 days before the Heavy snow in Gansu Province

在此基础上采用本文提出的历史平均值法对遥感影像像元进行补偿，该方法的核心是采用日期倒推，然后求得所有倒推天数内某个像元的非0值的像元值的和然后除以该像元非0值的天数，得出遥感影像像元值非0值天数的像元值的平均值，从而达到遥感影像像元云补偿的目的。最后将平均之后的甘肃省的整幅影像出图得到如下的补偿后的遥感影像图。

上述4幅图是按照本文提出的思路和方法按照随机选取的的4个节气采用历史平均值法进行补偿后的遥感影像图，为了便于直观看出K值的大小即甘肃省气温的高低变化情况，对补偿后的图进行分类，采用ArcGIS中的分位数的方法将其共分为10个类别(图8—图11)。图例中不同颜色的矩形颜色由浅入深代表着K值自上而下逐渐变大，也就意味着颜色越深K值越大，即温度越高。据图分析得出2017年大暑和立秋两个节气的遥感影像补偿后的K值比较大的区域出现在甘肃省北部地区，主要行政区范围包括敦煌市、瓜州县、玉门市、金塔县以及民勤县等地。其原因主要有两点：一是7—8月份正值北半球夏季，此时太阳直射点在北半球，且副热带高压北移；此时气候干燥，气温较高。二是沿着敦煌市、瓜州县、玉门市、金塔县以及民勤县这样的一个自西向东并逐渐转向东南的弧形地带恰巧也是库姆塔格沙漠、巴丹吉林沙漠南缘、腾格里沙漠南缘所形成的弧形带状区域，夏季干燥高温的环境加剧了该地带的高温，再加之沙土比热容小短期内温度升高较快，因而使得该地的K值偏大。K值较小的区域则是沿着阿尔金山、祁连山脉、甘南高原一线的这个带状区域，主要是甘肃省紧邻青海省的一侧，其原因主要有两点：一是该地带紧邻青藏高原气候类型属于高原高山气候，常年温度相对较低；二是该地带的平均海拔已经高于4 000 m，海拔加剧了该地带气候的寒冷。

图8 2017年大暑节气影像补偿图(白天)Fig.8 Image compensation of Great heat solar terms in 2017 (daytime)

图9 2017年立秋节气影像补偿图(白天)Fig.9 Image compensation of Autumn begins solar terms in 2017 (daytime)

图10 2017年大雪节气影像补偿图(白天)Fig.10 Image compensation of Heavy snow solar terms in 2017 (daytime)

图11 2017年冬至节气影像补偿图(白天)Fig.11 Image compensation of Winter solstice solar terms in 2017 (daytime)

2017年大雪和冬至两个节气的遥感影像补偿后的K值比较大的区域出现在甘肃省中东部地区并在甘肃省的西部地区有零星分布。中东部地区主要包括陇南、甘南、庆阳以及黄河谷地等地。西部地区主要是敦煌市，因为此地紧邻库姆塔格沙漠且是疏勒河的河谷地带。K值较小的区域则是甘肃紧邻青海省的一侧；主要是沿着阿尔金山、祁连山脉这个弧形的带状区域。K值较大的区域基本与甘肃省多年平均气温≥9 ℃的等温线经过的地带一致。

4 结论与讨论

4.1 结论

1)采用历史平均值法对遥感影像像元进行云补偿是切实可行的，从补偿的效果来看由补偿前的58.43%的平均值上升到补偿后接近100%的效果。因此，该方法科学合理且具有可重复性。

2)该方法对遥感影像的像元空值进行补偿相当于给原来表现为空值的像元赋一个值，经过计算发现补偿后的像元值最大值没有大于原遥感影像像元值的最大值。再次证明了补偿是在合理范围内的有效补偿。

3)本文以甘肃省为研究区，以甘肃省2017年MOD11A1数据为实验数据，运用历史平均值法对甘肃省的遥感影像像元云补偿获得了较好的补偿效果；甘肃省地貌类型丰富，涵盖了山地、高原、平川、河谷、沙漠、戈壁。自南向北依次跨越了湿润区、半湿润区、半干旱区以及干旱区；海拔从550 m的陇南白龙江中游文县罐子沟上升到5 547 m的祁连山主峰团结峰。甘肃省地理位置特殊、地域范围较广、地形复杂多样其囊括了几乎所有影响地表温度的因子，诸如：气温、太阳辐射、海陆位置、纬度、地表湿度、海拔、下垫面类型、植被覆盖率、人口密度以及工业的发展程度等。本文选取的研究区具有一定的代表性；因此，本文的研究结果和研究结论比较有向相关区域推广的潜力。

4.2 讨论

1)遥感影像像元空值的存在严重制约了人们对地表信息的识别和提取，本文所提出的历史平均值法是以某一天遥感影像(存在缺失)为立足点，尝试采用向前倒推时间，求得所推时间段内所有非0值像元的平均值作为被补偿像元0值对应区域的填充值。因此，该方法相当于对存在空值的遥感影像像元的真实值的模拟；因而用该方法求得的遥感影像像元值与卫星过境时所得的遥感影像像元值存在一定的误差。

2)采用历史平均值法对在获得的遥感影像上表现为空值的遥感影像像元值进行补偿的方法会受到倒推天数内的遥感影像像元值存在空值的天数的多少的影响，如果影像在所选择的倒推天数内都表现为空值，那么此时所选择的倒推天数就需要重新对其定义。而此时就需要考虑继续采用历史平均值法通过增长倒推天数的方式试图模拟出遥感影像像元值表现为空值的像元。

3)历史平均值法即便已经在MODIS11A1数据中有着良好的表现，在遥感卫星影像的其他类型数据中是否有效有待进一步研究。