集合卡爾曼濾波方法的高時空分辨率山區(qū)地表反照率反演

鄭凱旋， 林興穩(wěn)， 聞建光， 郝大磊

1. 浙江師范大學(xué) 地理與環(huán)境科學(xué)學(xué)院, 金華 321004;

2. 中國科學(xué)院空天信息創(chuàng)新研究院 遙感科學(xué)國家重點實驗室, 北京 100094;

3. 中國科學(xué)院大學(xué), 北京 100049

1 引 言

地表反照率（Albedo）定義為地表向各個方向反射的太陽短波輻射（0.3—3.0 μm）與太陽入射總輻射的比值（Liang等，2010），是表征地表反射太陽輻射能量強(qiáng)弱的物理量。地表反照率是輻射強(qiáng)迫的直接驅(qū)動因子之一，陸表覆蓋變化也是通過改變地表反照率，進(jìn)而影響輻射強(qiáng)迫，最終導(dǎo)致區(qū)域或全球氣候環(huán)境發(fā)生變化（Dickinson，1995；Trenberth 等，2009），因此，研究地表反照率的時空分布特征及其變化趨勢對研究全球或者區(qū)域氣候變化具有十分重要的意義。

遙感衛(wèi)星以其空間覆蓋廣的優(yōu)勢，為區(qū)域乃至全球尺度地表反照率的研究提供重要的數(shù)據(jù)支撐（Qu等，2015）。遙感衛(wèi)星輔以不同的地表反照率反演算法，構(gòu)建了長時間序列的地表反照率遙感產(chǎn)品（Liang等，2013）。在低分辨率尺度（公里級或百米級），地表反照率反演算法主要有基于地表二向反射分布函數(shù)（BRDF）模型算法（Schaaf 等，2002）、直接反演算法（Liu 等，2013）、靜止衛(wèi)星反演算法（孫越君 等，2018）。基于以上這些算法，生產(chǎn)了全球尺度的MODIS（Schaaf等，2002）、GLASS（Liang等，2013）和Copernicus（Lacaze等，2013）產(chǎn)品等；在高分辨率尺度（十米級），衛(wèi)星重訪周期長，短時間范圍內(nèi)難以累積足夠多角度的地表反射率，面臨著角度信息量不足難以滿足BRDF 建模需求（Gao 等，2014）的問題。目前，高分辨率地表反照率反演算法主要是利用低分辨率的BRDF 產(chǎn)品提供先驗知識，如基于MODIS BRDF 產(chǎn)品（MCD43A1）構(gòu)建的“pre-MODIS era”算法（Shuai 等，2011）、“MODIS-concurrent”算法（Shuai 等，2014）、利用Landsat 的光譜響應(yīng)函數(shù)和MODIS BRDF 數(shù)據(jù)庫構(gòu)建輻射傳輸模型的直接反演方法等（He等，2018）。上述算法在均一性較好的地表取得了較好的效果（Lin 等，2018a），但在山區(qū)地表下，地形改變了太陽—地形—傳感器之間的幾何關(guān)系，導(dǎo)致山區(qū)地表太陽輻射傳輸過程較平坦地表更為復(fù)雜（林興穩(wěn) 等，2020）。地形坡度、坡向、光照、陰影等效應(yīng)，導(dǎo)致同種地表類型在不同坡度和坡向下，地表反射率變化各異（Combal等，2000；Schaaf等，1994）以及地表BRDF 呈現(xiàn)出不對稱性（Wu 等，2019；Yin 等，2018）。目前BRDF 模型沒有考慮地形的影響，反演的地表反照率產(chǎn)品在山區(qū)地表存在一定的誤差（Lin 等，2018b）。林興穩(wěn)（2018b）將山地輻射傳輸模型和直接反演算法結(jié)合，初步構(gòu)建了適用于山區(qū)地表反照率遙感反演方法，并初步應(yīng)用于黑河地區(qū)的高分辨率環(huán)境衛(wèi)星的地表反照率的反演，但是山區(qū)云、雪多，云污染、季節(jié)性積雪等導(dǎo)致反演的地表反照率產(chǎn)品的缺失值較多，產(chǎn)品在空間分布上的存在缺失和時間序列上的存在不連續(xù)問題，嚴(yán)重制約山區(qū)地表反照率應(yīng)用。

低分辨率遙感反照率產(chǎn)品具有良好的時間連續(xù)性，但其空間分辨率過低（Holben，1986；Justice 等，1985）；高空間分辨率衛(wèi)星空間分辨率高，重訪周期長（Coops 等，2006）。許多學(xué)者將低空間分辨率遙感數(shù)據(jù)的高頻率觀測優(yōu)勢與高空間分辨率遙感數(shù)據(jù)的空間分辨率高的優(yōu)勢相結(jié)合，構(gòu)建了基于不同理論的時空融合算法（Li 等，2020）。這些算法大體上歸納為5 類；（1）基于權(quán)重函數(shù)的方法，如STARFM（Spatial and Temporal Adaptive Reflectance Fusion Model）系列模型，這類方法能夠準(zhǔn)確地捕捉物候變化，但精確性受限于景觀的特征斑塊大小，且采用簡單的線性模型，會影響山區(qū)地表參量反演（Gao等，2006）。（2）基于混合像元分解的方法，如MMT模型（Zhukov等，1999）、STDFA 模型（Wu 等，2012）等。該類模型采用混合像元分解模型，利用長時間序列的低空間分辨率影像提取地表反射率的非線性時間變化特征，進(jìn)而與高空間分辨率影像進(jìn)行融合（呂長春 等，2003）。但模型假設(shè)像元反射率不隨時間發(fā)生變化，與實際情況存在不符。（3）基于機(jī)器學(xué)習(xí)的方法，如SPSTFM 機(jī)器學(xué)習(xí)模型（Huang 和Song，2012），該方法有較好的穩(wěn)定性與可信度，但是算法過程中需要Landsat 與其對應(yīng)的MODIS 影像，導(dǎo)致算法在實際應(yīng)用中的難度增加（蔡家駿，2019）。（4）基于貝葉斯的方法，如TBDF 模型，該方法將時間相關(guān)信息融入到圖像時間序列中，將融合問題轉(zhuǎn)化為估計問題，通過最大后驗估計得到融合后的圖像，較為適合于異質(zhì)性地表（Xue等，2017）。（5）混合類方法，將基于權(quán)重函數(shù)、貝葉斯理論或者像元分解等方法相結(jié)合，得到更好的結(jié)果，如FSDAF 模型混合使用了像元分解及權(quán)重函數(shù)的方法（Zhu 等，2016）。上述方法大多未考慮山區(qū)地形，且受到景觀斑塊大小的影響，較難適用于山區(qū)地表時空連續(xù)的高空間分辨率地表反照率反演?？柭鼮V波算法（Kalman，1960），以分析誤差的最小方差為最優(yōu)標(biāo)準(zhǔn)，通過遞歸處理，可以保持過程結(jié)果，使之始終處于最優(yōu)狀態(tài)；Evensen（1994）提出了集合卡爾曼濾波資料同化的基本理論框架，克服了卡爾曼濾波僅限于處理線性問題的弱點，不再需要切線性模型和伴隨模式，解決了計算量巨大的問題。集合卡爾曼濾波EnKF（Ensemble Kalman Filter）算法通過集合預(yù)報與資料同化，在給定觀測和預(yù)報模式下去描述參量的概率分布及其發(fā)展（官元紅 等，2007），不受限于景觀的斑塊大小，故在山區(qū)也有很大的應(yīng)用空間。

本文首先基于改進(jìn)的直接反演算法獲取山區(qū)高分辨率地表反照率，利用集合卡爾曼濾波方法，構(gòu)建山區(qū)長時間序列高空間分辨率反照率產(chǎn)品，解決反演的反照率在山區(qū)地表下時空缺失值較多的情況。最后，將反演的反照率產(chǎn)品與地面站點實測數(shù)據(jù)進(jìn)行對比，評價構(gòu)建的時空連續(xù)地表反照率產(chǎn)品的精度和不確定性。

2 研究區(qū)與數(shù)據(jù)源

2.1 研究區(qū)域

黑河流域是中國西北地區(qū)第二大內(nèi)陸流域，位于河西走廊農(nóng)牧交錯帶上，橫跨3 種不同的自然環(huán)境單元，即上游地區(qū)主要為祁連山區(qū)、中游地區(qū)是平原走廊、下游地區(qū)則以低山山地為主，整個流域地區(qū)是典型的生態(tài)敏感區(qū)（趙曉冏，2012）。本研究區(qū)位于黑河流域南部，整個研究區(qū)高程范圍為889—5218 m，太陽輻射強(qiáng)烈，地表覆蓋類型包含農(nóng)田，林地，草地、濕地、裸地、冰雪覆蓋區(qū)等，主導(dǎo)地表類型為裸地及草地（圖1）。該地區(qū)因獨特的地理位置及氣候條件，是寒區(qū)和旱區(qū)生態(tài)、水文陸表過程研究的理想?yún)^(qū)域。大量科學(xué)研究與實驗已經(jīng)積累非常豐富的數(shù)據(jù)資料，已開展的科學(xué)實驗包括“非均勻下墊面地表蒸散發(fā)的多尺度觀測試驗”（Liu等，2018）、“無人機(jī)物聯(lián)網(wǎng)中繼實現(xiàn)無公網(wǎng)偏遠(yuǎn)地區(qū)環(huán)境監(jiān)測”（Zhang 等，2020）、“黑河流域生態(tài)—水文過程綜合遙感觀測聯(lián)合試驗（HiWater）”等（李新 等，2012）。其中，HiWater 實驗是黑河流域的多尺度綜合觀測試驗，構(gòu)建了長時間序列的水文氣象觀測網(wǎng)絡(luò)。該水文氣象觀測網(wǎng)包含有不同下墊面的地面觀測站點，提供多種觀測數(shù)據(jù)集，包含風(fēng)速、氣壓、四分量輻射、地表輻射溫度、光合有效輻射、土壤熱通量、土壤水分、土壤溫度等，每10 min 記錄一次并對外發(fā)布。

圖1 黑河流域研究區(qū)高程及土地利用類型圖Fig. 1 Elevation and land cover maps of Heihe Watershed area

2.2 數(shù)據(jù)源

2.2.1 地面數(shù)據(jù)

地面驗證數(shù)據(jù)采用青藏高原科學(xué)數(shù)據(jù)中心發(fā)布的2016 年和2017 年“黑河流域生態(tài)—水文過程綜合遙感觀測聯(lián)合試驗”地面輻射觀測數(shù)據(jù)集（Liu 等，2018）。由于MODIS 反照率已訂正到太陽正午時刻的值（Wu等，2017）。考慮到黑河地區(qū)的正午時刻約為下午13：00時左右，故選取12：30 pm至13：30 pm 之間觀測的地表上行短波輻射值與太陽下行短波輻射值（觀測波段范圍：0.3—3.0 μm），計算比值并求平均值作為每天的反照率觀測值，計算式為式（1）：

式中，ar表示日均當(dāng)?shù)卣鐣r刻地表短波反照率，ri表示每10 min短波上行輻射，表示每10 min短波下行輻射，n＝7。為驗證算法在不同地表的適用情況，選取了4 個不同下墊面的站點進(jìn)行算法檢驗；考慮到MODIS反照率產(chǎn)品及GF-4數(shù)據(jù)像元大小，統(tǒng)計站點周圍2 km 范圍內(nèi)的坡度分布范圍，站點信息見表1。

表1 輻射通量觀測站點信息Table 1 information of radiation flux observation site

2.2.2 衛(wèi)星數(shù)據(jù)

GF-4 衛(wèi)星是中國第一顆地球同步軌道遙感衛(wèi)星，可見光和多光譜空間分辨率優(yōu)于50 m，紅外譜段空間分辨率優(yōu)于400 m，實現(xiàn)對中國及周邊地區(qū)的觀測（孫越君 等，2018），GF-4 傳感器各波段參數(shù)見表2，本文所用GF-4 反射率數(shù)據(jù)見表3。經(jīng)過輻射定標(biāo)等預(yù)處理的GF-4 反射率數(shù)據(jù)用于坡面反照率反演。

表2 GF-4衛(wèi)星傳感器參數(shù)Table 2 Observation parameters of GF-4 satellite sensor

表3 GF-4影像數(shù)據(jù)Table 3 Information of GF-4 data

MCD43A3 C6 產(chǎn)品是日尺度16 d 合成的500 m地表反照率產(chǎn)品。是通過篩選近16 d質(zhì)量較好的像元，組成多角度觀測數(shù)據(jù)集，進(jìn)而擬合BRDF 得到地表反照率（Schaaf等，2002）；該產(chǎn)品經(jīng)過廣泛測試使用，精度已得到普遍認(rèn)可（馮智明 等，2018）。

地表高程數(shù)據(jù)采用國家青藏高原科學(xué)數(shù)據(jù)中心發(fā)布的ASTER 衛(wèi)星90 m 地表高程數(shù)據(jù)（https：//data.tpdc.ac.cn/zh-hans/data［2020-08-03］），地表覆蓋類型采用清華大學(xué)發(fā)布的2017年30 m 地面分類數(shù)據(jù)（http：//data.ess.tsinghua.edu.cn/［2020-08-03］），地表覆蓋類型共分為10 類，包括農(nóng)田、森林、草地、灌木地、濕地、水、苔原、不透水面、裸地和冰雪（圖1）。

3 方法構(gòu)建

3.1 算法流程

算法流程圖如圖2所示，首先利用山區(qū)地表反照率反演方法得到GF-4 坡面地表反射率（Lin 等，2018b）；然后通過研究區(qū)坡度、坡向、地表高程數(shù)據(jù)、土地利用數(shù)據(jù)對MCD43A3 數(shù)據(jù)進(jìn)行逐像元處理，利用高空間分辨率地表覆蓋數(shù)據(jù)提取純像元，并將MCD43A3 純像元的觀測值作為時空連續(xù)地表反照率反演的背景場；在此基礎(chǔ)上代入地面站點觀測數(shù)據(jù)和GF-4 坡面地表反照率，通過EnKF 算法得到長時間序列高分辨率反照率數(shù)據(jù)；最后進(jìn)行反演結(jié)果的精度驗證。

圖2 長時間序列山區(qū)地表反照率遙感反演流程圖Fig. 2 The flow chart of long time series land surface albedo retrieval over rugged terrain

3.2 山區(qū)地表反照率的直接反演算法

Lin 等（2018b）基于單一角度反照率反演算法，初步構(gòu)建了山區(qū)地表反照率反演算法。該方法假設(shè)高分辨率像元內(nèi)地表具有單一的坡度和坡向，考慮樹木向地性生長；在斜坡上構(gòu)建坡面坐標(biāo)系，將太陽入射和衛(wèi)星觀測幾何格網(wǎng)化，構(gòu)建相對太陽入射天頂角、相對太陽入射方位角、相對衛(wèi)星觀測天頂角和相對衛(wèi)星觀測方位角的格網(wǎng)組合。在坡面坐標(biāo)系的每個格網(wǎng)中，利用核驅(qū)動模型，構(gòu)建不同地表類型下的地表反射率和寬波段地表反照率之間的查找表（Liu 等，2013；孫長奎 等，2013），進(jìn)而在坡面坐標(biāo)系中反演坡面地表反照率。

式中，ρdd（is，φs，iv，φv）、ρhd（iv，φv）、ρdh（is，φs）、ρhh分別表示為坡面像元方向—方向反射率、半球—方向反射率、方向—半球反射率和半球—半球反射率。根據(jù)定義可將ρhd（iv，φv）、ρdh（is，φs）、ρhh展開如下：

將式（3）—（5）分別帶入到式（2）中，構(gòu)建關(guān)于ρdd（is，φs，iv，φv）的一元二次方程，利用各地表類型的MODIS純像元提供BRDF先驗知識，在坡面坐標(biāo)系下，對坡面方向—方向反射率進(jìn)行積分獲取坡面半球—方向、方向—半球和半球—半球反射率，進(jìn)而求解關(guān)于ρdd（is，φs，iv，φv）的一元二次方程，獲取坡面地表反射率。利用構(gòu)建的坡面地表反射率和寬波段地表反照率的查找表，直接反演得到地表反照率（式（6））

式中，αABT表示高分辨率坡面反照率，is，φs，iv，φv分別表示坡面坐標(biāo)系上相對的太陽天頂角、太陽方位角、觀測天頂角和觀測方位角；ρi（is，φs，iv，φv）表示第i波段的地表反射率；Ci（is，φs，iv，φv）表示坡面坐標(biāo)系下的回歸系數(shù)。

3.3 山區(qū)高分辨率地表反照率時空填補方法

本文采用改進(jìn)的集合卡爾曼濾波方法對反演的山區(qū)地表反照率遙感產(chǎn)品進(jìn)行處理，以MODIS地表反照率產(chǎn)品和地面實測數(shù)據(jù)為先驗知識，構(gòu)建時空連續(xù)的山區(qū)地表反照率產(chǎn)品。集合卡爾曼濾波方法常用來進(jìn)行資料同化，充分有效的利用已有的各種信息預(yù)測未知的狀態(tài)，關(guān)鍵步驟是計算過程中的狀態(tài)更新（朱琳，2007）；應(yīng)用在山區(qū)高分辨率地表反照率反演中，有兩點需要考慮，一是背景場的構(gòu)建，二是集合卡爾曼濾波方法構(gòu)建。

3.3.1 先驗知識背景場的構(gòu)建

MODIS 反照率數(shù)據(jù)未考慮地形影響，將其作為背景場時，須先對其進(jìn)行像元尺度的處理，增強(qiáng)其在山區(qū)的適用性。本文首先以GF-4 影像空間分辨率和坐標(biāo)系統(tǒng)為基準(zhǔn)，對MODIS 反照率產(chǎn)品進(jìn)行重采樣與投影變換；然后，在GF-4 影像上建立搜索窗口，搜索與窗口中心像元相似的像元，確定每個像元周圍的相似像元位置，主要從兩方面進(jìn)行相似像元的選?。阂皇桥c窗口中心像元具有相同地物覆蓋類型、相似高程及坡度坡向的像元視為相似像元；二是計算中心像元與鄰近像元的反照率差異值，通過閾值判斷差異是否合理，在合理區(qū)間內(nèi)的像元被視為相似像元。閾值可以通過高分辨率圖像中像元總體的標(biāo)準(zhǔn)差和估計的土地覆蓋類別數(shù)來確定（式（7））（Gao等，2006）；

式中，F(xiàn)i、Ci分別表示最鄰近時刻第i個相似像元位置對應(yīng)的高、低空間分辨率的反照率值，cov（.）表示求協(xié)方差，D（Fi）、D（Ci）分別表示Fi、Ci的方差，R的范圍［-1，1］，R值越大表示越相似。第i個相似像素與中心像素的地理距離di可由式（9）計算，窗口內(nèi)相似像素的距離范圍［1，1+20.5］（Zhu等，2010）。

結(jié)合光譜相似性和距離，可以計算出光譜和地理距離結(jié)合的綜合指數(shù)D（式（10））：

式中，D越大，表示相似像元的貢獻(xiàn)越小，利用倒數(shù)歸一化形式（式（11）），得到權(quán)重W表示為

式中，Wi的范圍為0 到1，所有相似像素的總權(quán)重為1；故可計算tp時刻搜索窗口中心像元的高分辨率反照率（式（12））：

式中，N為搜索窗口內(nèi)通過式（7）得到的近似像元個數(shù)，包含中心預(yù)測像元；C為相似像元位置對應(yīng)的重采樣后MODIS 產(chǎn)品反照率值；vi為第i個相似像元的轉(zhuǎn)換系數(shù)，是通過線性回歸分析窗口內(nèi)GF-4 反照率數(shù)據(jù)和MODIS 反照率數(shù)據(jù)的相似像元，建立回歸模型，取其斜率得到。

3.3.2 集合卡爾曼濾波算法

集合卡爾曼濾波算法根據(jù)背景場以及觀測值的特征誤差分布，對背景場和觀測值進(jìn)行分析，得到一組分析值，再用這組分析值的差異作為分析誤差的統(tǒng)計樣本進(jìn)行分析誤差協(xié)方差的估計；對這組分析值作一個短期預(yù)報后，得到一組預(yù)報值。把這組預(yù)報值的差異作為背景場誤差的統(tǒng)計樣本進(jìn)行誤差協(xié)方差的估計，以此遞歸處理整個系統(tǒng)（Li 等，2012）。當(dāng)測量及觀測誤差已知時，集合卡爾曼濾波器可以從這些數(shù)據(jù)中估計系統(tǒng)的動態(tài)變化，其動態(tài)模型一般表示為

式中，X表示動態(tài)模型狀態(tài)向量，本文中為高空間分辨率的地表反照率，Xk?1，k表示預(yù)測的高空間分辨率反照率；Mk?1，k是k?1到k時刻的地表反照率的狀態(tài)轉(zhuǎn)移矩陣，考慮前后兩日正午時刻地表反照率之間互不影響，且變化較小，此處將M設(shè)為1。Bk?1，k是k?1 到k時刻的系統(tǒng)控制矩陣，Uk表示k時刻的狀態(tài)控制量，本研究是直接獲取地表反照率，沒有控制量，所以Uk=0，Bk不需要考慮。vk-1為k?1 時刻的觀測誤差，符合正態(tài)分布。通過式（14）對k時刻的高分辨率地表反照率Xk進(jìn)行預(yù)測：

式中，Kg為卡爾曼增益矩陣；Zk表示集合卡爾曼濾波算法中的實際測量值，在本研究中即為地面站點實際測量反照率值；Hk是測量系統(tǒng)的參數(shù)，將狀態(tài)向量轉(zhuǎn)化為觀測類型的值，在本研究中，觀測值與狀態(tài)向量均為地表反照率，即已直接獲得反照率的測量結(jié)果，故H設(shè)為單位矩陣。

式中，Zk是k時刻經(jīng)過處理的地面站點實際測量反照率值。Kgk可以通過式（16）表示：

式中，R是測量噪聲協(xié)方差矩陣，HT表示H的轉(zhuǎn)置矩陣。由上一預(yù)測時刻的誤差協(xié)方差Pk?1和過程噪聲矩陣Q預(yù)測新的誤差協(xié)方差矩陣Pk-1，k，可表示如下：

式中，Pk?1，k是Xk對應(yīng)的協(xié)方差矩陣，Q是系統(tǒng)過程的噪聲協(xié)方差，MT表示M的轉(zhuǎn)置矩陣。此時得到了k狀態(tài)下最優(yōu)的估算值X。為使卡爾曼濾波器不斷的運行下去直到系統(tǒng)過程結(jié)束，還需更新k狀態(tài)下X的協(xié)方差Pk：

式中，I為單位矩陣。

在本研究中，以3.1 節(jié)得到的結(jié)果為背景場。考慮到研究區(qū)域內(nèi)并不是每個像素都具有地面測量值，故計算區(qū)域內(nèi)測量誤差的協(xié)方差矩陣R時，基于站點誤差構(gòu)建了均值為0，方差為0.01 的正態(tài)分布矩陣，作為區(qū)域內(nèi)像素的測量誤差；即根據(jù)站點像素的測量誤差設(shè)置區(qū)域的測量誤差，進(jìn)而計算區(qū)域內(nèi)測量誤差的協(xié)方差矩陣R；同理，在預(yù)測新的誤差協(xié)方差矩陣P的過程中（公式（17）），將系統(tǒng)過程噪聲設(shè)置為均值為0，方差為0.01 的矩陣，求得系統(tǒng)過程噪聲協(xié)方差矩陣Q，進(jìn)而得到新的誤差協(xié)方差矩陣P（Zhang等，2019）。

4 結(jié)果與驗證

4.1 山區(qū)時空連續(xù)地表反照率反演結(jié)果

本文算法利用日尺度MODIS 反照率產(chǎn)品構(gòu)建先驗知識背景場，利用地面觀測數(shù)據(jù)和集合卡爾曼濾波算法，結(jié)合數(shù)幅高質(zhì)量的GF-4 影像，獲得日尺度50 m 的反照率產(chǎn)品。結(jié)果顯示，本文的方法可以構(gòu)建時空連續(xù)的高分辨率地表反照率產(chǎn)品，且隨著高質(zhì)量GF-4 影像的增加，同化效果得到提升。為體現(xiàn)算法的普適性，選取四個不同地表覆蓋類型的站點進(jìn)行結(jié)果驗證。根據(jù)各站點首次獲得最優(yōu)高分辨率反照率產(chǎn)品為起始日期，下圖展示了各站點地表反照率隨時間的變化趨勢（圖3）；由于大滿及花寨子站點距離較近，在圖3（c）中用一組圖展示。

圖3 2017年研究區(qū)內(nèi)不同類型地表反照率反演結(jié)果圖Fig. 3 Shortwave albedo maps generated by EnKF for different land surface types of study area in 2017

圖3（a）為阿柔研究區(qū)2017 年第100—300 D地表反照率變化圖，其下墊面為亞高山山地草甸，地表反照率變化最為明顯，尤其是中下部分。中部顏色由褐色逐漸變?yōu)辄S色，下半部分深綠色逐漸增加，表明隨著時間變化，該區(qū)域地表反照率在下降；中部下墊面以草甸為主，此外還存在部分農(nóng)田，隨著草甸回綠及農(nóng)作物的生長，地表反照率不斷下降至0.2 左右；南部山區(qū)人跡罕至，下墊面包含森林、草地和裸地，圖3（a）中青色部分地表反照率較小，該區(qū)域為森林覆蓋，反照率較周圍草地、裸地等要小。圖中紅色部分地表高程高，下墊面以草地和裸地為主，部分地區(qū)山頂常年積雪，反照率值逼近0.6，但隨著季節(jié)變化，大部分山頂冰雪消融，地表反照率值下降；圖3（b）為張掖站點，深綠色的上半部分為蘆葦濕地，海拔較低，在該時間段生態(tài)情況較為穩(wěn)定，變化浮動較小，穩(wěn)定在0.17 左右。下半部分下墊面為裸地，常年地表反照率穩(wěn)定保持在0.24左右；圖3（c）中部為大滿氣象站點所在地，下墊面為玉米地，該區(qū)域在2017 年前100—240 d 左右地表反照率一直處于下降狀態(tài)，從0.23 左右降到0.15 附近，在圖中體現(xiàn)為中間部分綠色一直不斷加深，對應(yīng)玉米的灌漿期、拔節(jié)、封壟期地表反照率的變化；接著反照率出現(xiàn)上升，應(yīng)是因為在9月底左右農(nóng)作物處于收割階段，大量農(nóng)作物的收割導(dǎo)致該現(xiàn)象的出現(xiàn)；圖3（c）下方為花寨子站點，黃色部分的下墊面為山前荒漠，該區(qū)域在時間尺度上變化不大，反照率穩(wěn)定在0.23 左右；左下角綠色區(qū)域不斷加深，為植被區(qū)（姚彤和張強(qiáng)，2014）。圖4 展示了阿柔站點2017年第300 d地表高程與反演反照率的局部放大對比圖，從中可以直觀的發(fā)現(xiàn)山區(qū)坡度坡向變化引起的反照率變化；在圖4中，高程較低的區(qū)域下墊面以草甸為主，此外還存在部分農(nóng)田，為地表反照率值較大，隨著高程的增加，陽坡地表類型以草地位置，地表反照率較大，陰坡地表以針葉林為主，地表反照率較陽坡地表低。在地表高程更高的地方存在積雪區(qū)域，反照率較高，部分地區(qū)的地表反照率值要大于0.4。

圖4 地表反照率反演結(jié)果隨地表高程變化Fig. 4 The distribution of shortwave albedo retrieved by EnKF following the change of elevation

圖5 為阿柔、張掖、大滿、花寨子研究區(qū)2017 年時間段內(nèi)反照率及與地面實測地表反照率隨時間變化趨勢圖，從中可知，除了張掖濕地站點外，在其余3個站點下，MODIS產(chǎn)品較地面實測反照率存在低估現(xiàn)象，但在時間序列上，MCD43A3產(chǎn)品的變化趨勢與地面站點相同，能較好的反映地表的季節(jié)性變化特征。張掖站點的MCD43A3 產(chǎn)品較地面觀測值偏高，該站點的下墊面為濕地，地面站點觀測值可能包含來自水體的反射信號。MCD43A3 產(chǎn)品的空間分辨率為500 m，在該站點下兩者之間的尺度效應(yīng)明顯?；ㄕ诱军c的地面觀測數(shù)據(jù)波動性較大，通過地面實測數(shù)據(jù)的時間序列分析發(fā)現(xiàn)，花寨子站點前后相鄰日期的地表反照率變化較大，可能是站點觀測儀器穩(wěn)定性、觀測條件等因素影響，導(dǎo)致地面反照率降低，進(jìn)而影響了實際測量，但整體局勢符合實際情況。在這4個站點上，GF-4影像較少，利用GF-4衛(wèi)星直接反演的地表反照率與地面站點的一致性較好。利用本文算法得到的時空連續(xù)的地表反照率與地面站點的一致性也比較好，比MODIS 產(chǎn)品更能反映地表反照率的時空變化，準(zhǔn)確性也更高，同時也可以很好的保持時間上的連續(xù)性。比較本文算法反演的2017 年4 個站點的時空連續(xù)的地表反照率與站點實測地表反照率可知，4個站點的R2分別為0.8701、0.9865、0.9882、0.4466，大多高達(dá)0.8；RMSE分別為0.0116、0.0036、0.0036、0.0123，皆小于0.02。

圖5 2017年4個研究區(qū)反演結(jié)果與地面實測對比圖Fig. 5 Comparison between the inversion results and four observation sites measured data

4.2 預(yù)測地表反照率和直接反演地表反照率比較結(jié)果

圖6—圖8展示的是本文算法預(yù)測的GF-4地表反照率與同一天的利用當(dāng)天已有反射率數(shù)據(jù)采用直接反演算法得到的GF-4 地表反照率的比較結(jié)果。張掖、大滿及花寨子站點對比日期選取了2017 年6 月12 日，阿柔站點選擇了2016 年12 月15 日。從各站點反演結(jié)果與反照率產(chǎn)品對比圖可以看出，本文預(yù)測的地表反照率與直接反演的GF-4 地表反照率的空間分布具有很高的相似性。與直接反演算法得到的GF-4 反照率數(shù)據(jù)相比，模擬結(jié)果在不同的坡度和坡向下特征明顯，表現(xiàn)出地形特征效果。在不同下墊面下，尤其是阿柔、張掖等站點存在較大坡度，但站點反演結(jié)果與GF-4 反照率數(shù)據(jù)有較好的可比性，整體反演情況較好；說明本方法在模擬過程中能有效地利用已有數(shù)據(jù)，對崎嶇地表完成較好的模擬。反演結(jié)果與反照率產(chǎn)品對比圖顯示張掖站點的計算結(jié)果與同時期GF-4 反照率數(shù)據(jù)值集中在0.17—0.32，統(tǒng)計結(jié)果顯示兩者之間的RMSE 為0.0045，R2為0.9623；大滿及花寨子站點的計算結(jié)果與GF-4反照率的對比，RMSE為0.0096，R2為0.7372；阿柔站點的預(yù)測的地表反照率與GF-4直接反演的地表反照率的對比，RMSE為0.0163，R2為0.9163。RMSE 皆小于0.02，R2大多高達(dá)0.8，這說明了本文算法預(yù)測的地表反照率精度較好，可以使用于長時間序列的地表反照率的構(gòu)建。

圖6 張掖站點2017年6月12日GF-4反照率產(chǎn)品與反演結(jié)果對比Fig. 6 Comparison between the GF-4 albedo and inversion results of Zhangye site on June 12， 2017

圖7 大滿及花寨子站點2017年6月12日GF-4反照率產(chǎn)品與反演結(jié)果對比Fig. 7 Comparison between the GF-4 albedo and inversion results of DaMan and HuaZhaizi site on December 12， 2017

圖8 阿柔站點2016年12月15日GF-4反照率產(chǎn)品與反演結(jié)果對比Fig. 8 Comparison between the GF-4 albedo and inversion results of ARou site on December 25， 2016

5 結(jié) 論

本文參考前人工作提出一種山區(qū)地表時間序列反照率的估算方法，以地面站點觀測數(shù)據(jù)和時空連續(xù)的MODIS反照率產(chǎn)品構(gòu)建先驗知識背景場，結(jié)合少量的GF-4 數(shù)據(jù)通過集合卡爾曼率波算法得到黑河流域長時間序列高時空分辨率的反照率產(chǎn)品。通過與地面站點地表反照率測量數(shù)據(jù)對比，在不同地表坡度和地表類型下，反演的地表反照率與實測值時空一致性較好，RMSE 皆小于0.02，R2高達(dá)0.8。與直接反演的GF-4 反照率數(shù)據(jù)相比，預(yù)測的地表反照率與同期的地表反照率產(chǎn)品之間的一致性較好，全景影像的統(tǒng)計結(jié)果顯示兩幅影像的RMSE 小于0.02?？紤]到GF-4 衛(wèi)星高空間分辨率的特點，通過增加GF-4 影像數(shù)量，能夠得到更精確的預(yù)測結(jié)果。本文的研究具有以下創(chuàng)新點：第一，考慮了山區(qū)地表的地形效應(yīng)，在通過鄰近像元預(yù)測反照率時，通過選取相似像元，考慮了坡度坡向、高程、地表分類等信息，增強(qiáng)了背景場在山區(qū)的適用性，同時發(fā)揮其時間連續(xù)性優(yōu)勢，以此作為背景場來反映地表隨時間變化的物候規(guī)律；第二：利用集合卡爾曼濾波算法，通過使用較為少數(shù)的高精度影像，結(jié)合MODIS反照率產(chǎn)品，利用算法生成計算出全年的高時空分辨率反照率產(chǎn)品，有效的提高了山區(qū)反照率反演精度。

本文能夠為其他衛(wèi)星數(shù)據(jù)或者生物量的反演提供一種思路，然而具體在其他衛(wèi)星上的應(yīng)用效果還需進(jìn)一步檢驗；考慮到已發(fā)布的地面實測數(shù)據(jù)及衛(wèi)星數(shù)據(jù)，僅使用較為少數(shù)的站點進(jìn)行算法驗證，后續(xù)應(yīng)增加站點數(shù)量及觀測年份，提高算法普適性。