歷史平均值法用于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及與之相關的工具試圖尋找新的具有可操作性和可重復性的的遙感影像云補償方法?；诖诵枰鉀Q補償數據的選取、處理、拼接及適用于本思路的時段劃分；影像精度設定、代碼條件設定、代碼的書寫、運行和運行結果的判定；影像像元有效率的計算及統計，然后運用歷史平均值法進行補償，將遙感影像像元云補償后的結果與補償前的遙感影像進行對比，得出影像云補償的精度及優點，并試圖探索此方法是否也同樣適用于其他衛星的數據類型。

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數據中有著良好的表現，在遙感衛星影像的其他類型數據中是否有效有待進一步研究。