基于改進(jìn)PT-hybrid算法的地表潛熱通量變化分析

周耀坤，邢萬(wàn)秋

（河海大學(xué) 水文水資源學(xué)院，南京 210098）

0 引 言

【研究意義】地表潛熱通量（LE）是指從陸地表面到大氣中水分蒸發(fā)和植被蒸騰的熱通量，對(duì)于理解地球表面能量收支和地表-大氣相互作用具有重要意義[1-2]。LE的準(zhǔn)確模擬對(duì)于了解大范圍內(nèi)各生態(tài)系統(tǒng)的穩(wěn)定性起著至關(guān)重要的作用。【研究進(jìn)展】目前，通過(guò)在全球范圍內(nèi)建立的渦流協(xié)方差通量塔和衛(wèi)星觀測(cè)站，大量的LE模型和算法被開(kāi)發(fā)以適應(yīng)不同生態(tài)系統(tǒng)的要求[3-5]。Priestley-Taylor (PT)算法是一種較為簡(jiǎn)化的算法，該算法避免了使用表面阻抗，算法中參數(shù)產(chǎn)生的誤差減少，從而提高了LE估算的準(zhǔn)確性[6]。由于PT 算法本身并沒(méi)有考慮下墊面植被類型的影響，Yao等[7]基于Merry 再分析數(shù)據(jù)和Modis 衛(wèi)星數(shù)據(jù)對(duì)PT算法進(jìn)行了改進(jìn)，該研究在原有的PT 算法中引入了下墊面不同的植被功能類型，并給每一種植被類型率定了一套參數(shù)。研究表明，考慮下墊面特征的PT 算法能有效提高LE模擬的準(zhǔn)確性。【切入點(diǎn)】然而，目前大多數(shù)的算法都只考慮了單一像元，即一個(gè)柵格內(nèi)往往只考慮一種植被類型，但在陸地生態(tài)系統(tǒng)較為復(fù)雜的區(qū)域，考慮單一像元的算法會(huì)產(chǎn)生較大的誤差[8]。【擬解決的關(guān)鍵問(wèn)題】為了避免這種誤差，本研究結(jié)合一種動(dòng)態(tài)植被模型—LPJ-GUESS 模型，并引入亞像元法，對(duì)PT 算法進(jìn)行改進(jìn), 利用改進(jìn)的算法對(duì)全球陸地LE進(jìn)行反演，并分析全球植被與陸地LE的時(shí)空分布特征。

1 數(shù)據(jù)來(lái)源

1.1 氣象數(shù)據(jù)和遙感數(shù)據(jù)

本文所選取的氣象數(shù)據(jù)為歐洲天氣預(yù)報(bào)中心 （ECMWF）提供的空間分辨率為0.5°×0.5°的全球再分析氣象數(shù)據(jù)集，時(shí)間跨度為 1982—2014年（http://apps.ecmwf.int/datasets）。

本文選取的遙感數(shù)據(jù)主要為美國(guó)航空航天局（NASA）提供的全球農(nóng)業(yè)監(jiān)測(cè)系統(tǒng)（GIMMS）全球0.5°×0.5°歸一化植被指數(shù)（NDVI）1982—2014年數(shù)據(jù)集（https://glam1.gsfc.nasa.gov/）。本文選取的氣象數(shù)據(jù)和遙感數(shù)據(jù)都具有長(zhǎng)時(shí)間序列的優(yōu)勢(shì)，滿足本文的研究需求。

1.2 渦度通量數(shù)據(jù)

渦度通量數(shù)據(jù)主要用于驗(yàn)證和評(píng)價(jià)算法改進(jìn)前后的性能。本研究采取FLUXNET2015 數(shù)據(jù)集，該數(shù)據(jù)集包括從全球多個(gè)區(qū)域通量站點(diǎn)收集的數(shù)據(jù)（https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/），本次研究選取了其中49 個(gè)通量站點(diǎn)，如圖1 所示。數(shù)據(jù)主要包括2012—2014年的LE年值序列，用于驗(yàn)證改進(jìn)前后的算法。

1.3 其他數(shù)據(jù)

LPJ-GUESS所使用的CO2質(zhì)量濃度數(shù)據(jù)來(lái)源于1980—2000年的冰芯觀測(cè)[9]，土壤類型數(shù)據(jù)來(lái)源于聯(lián)合國(guó)糧食和農(nóng)業(yè)組織土壤數(shù)據(jù)集。

2 研究方法

2.1 LPJ-GUESS 模型介紹

LPJ-GUESS 模型是一種針對(duì)區(qū)域和全球應(yīng)用進(jìn)行優(yōu)化的動(dòng)態(tài)植被模型[10]。模型在全球動(dòng)態(tài)植被模型LPJ-DGVM 的基礎(chǔ)上發(fā)展而來(lái)，對(duì)物種個(gè)體觀測(cè)水平的模擬進(jìn)行了細(xì)化，該模型能夠有效模擬氣候變化下植被的瞬時(shí)響應(yīng)，常用于描述各種植被地理分布格局的動(dòng)態(tài)變化、生理過(guò)程、生物物理和生物地球化學(xué)過(guò)程等[11]。模型主要輸入數(shù)據(jù)有月平均氣候（氣溫、降水和日照百分率）、年平均CO2質(zhì)量濃度以及土壤類型數(shù)據(jù)。本研究基于該模型模擬地表分布最廣泛的5 種植被類型，落葉闊葉林（DBF）、落葉針葉林（DNF）、常綠闊葉林（EBF）、常綠針葉林（ENF）、稀樹(shù)草原（GRA），分析地表植被功能類型變化對(duì)于陸地蒸發(fā)能力的影響。

圖1 研究區(qū)49 個(gè)渦度通量塔的空間分布 Fig.1 The spatial distribution of 49 eddy covariance flux tower sites used in this study

2.2 PT-hybrid 算法以及算法改進(jìn)介紹

本研究是在最近的PT-hybrid 算法基礎(chǔ)上，考慮了下墊面植被動(dòng)態(tài)變化以及加入亞像元方法，對(duì)現(xiàn)有模型進(jìn)行改進(jìn)。PT-hybrid 算法是由PT 模型發(fā)展而來(lái)，其計(jì)算式為：

式中：LE為潛熱通量(W/m2)；?為濕表面狀況的PT系數(shù)，值為1.26；Δ為飽和蒸汽壓曲線的斜率（kPa/℃）；γ為濕度計(jì)常數(shù)（kPa/℃）；Rn和G分別代表了太陽(yáng)凈輻射和土壤熱通量（W/m2）；ki（i=0,…,4）是經(jīng)驗(yàn)系數(shù)。其中，G是由植被覆蓋率fc和Rn根據(jù)經(jīng)驗(yàn)方法得到的，計(jì)算式為：

式中：ag值為0.18；NDVImin和NDVImax分別為研究時(shí)段內(nèi)的NDVI的最小值和最大值，設(shè)定為常數(shù)：0.05和0.95[12]。

相比于其他物理算法而言，本研究選用PT-hybrid算法具有以下幾點(diǎn)優(yōu)勢(shì)：①該算法只需要Rn、Ta、NDVI、RH和VPD等氣象要素，簡(jiǎn)化了輸入?yún)?shù)，更適用于長(zhǎng)序列時(shí)間尺度下的LE模擬；②該算法避免使用溫度差或濕度差，并且不考慮空氣動(dòng)力和表面阻抗，通過(guò)削弱輸入?yún)?shù)帶來(lái)的誤差大幅度減小了算法結(jié)果誤差；③該算法針對(duì)不同的下墊面植被類型率定了不同的參數(shù)，如表1 所示，從而提高了LE模擬的準(zhǔn)確度。

表1 基于全球植被功能類型優(yōu)化所得系數(shù) Table 1 Coefficients based on global vegetation functional type optimization

圖2 亞像元法解釋圖 Fig.2 Sketch map of sub-pixel method

然而在該算法中，輸入的地表植被類型是靜態(tài)且恒定的。考慮到下墊面植被會(huì)因氣候變化和人類活動(dòng)而發(fā)生變化，本研究引入動(dòng)態(tài)變化的植被類型來(lái)代替土地覆蓋類型產(chǎn)品作為算法驅(qū)動(dòng)數(shù)據(jù)。另外，引入亞像元法，避免了原算法運(yùn)用單一像元在陸地生態(tài)系統(tǒng)較為復(fù)雜的地區(qū)而產(chǎn)生較大的誤差，如圖2 所示。

圖2（a）為以往算法輸入的植被覆蓋類型，一個(gè)柵格只考慮了一種植被類型，而大多數(shù)區(qū)域下墊面植被類型應(yīng)呈現(xiàn)圖2（b）的分布，因此，引入亞像元法將柵格細(xì)分，根據(jù)LPJ-GUESS 模型的輸出—不同植被功能類型在每一個(gè)柵格的葉面積指數(shù)（LAI），用不同植被類型的葉面積指數(shù)占每一個(gè)柵格總?cè)~面積指數(shù)的比重來(lái)表征不同植被功能類型在每一個(gè)柵格的占比，通過(guò)考慮同一柵格不同植被占比變化來(lái)減少算法誤差。本研究引入?yún)?shù)fv來(lái)表征不同植被功能類型的占比。公式為：

式中：fv是由LPJ-GUESS 模型模擬得到的不同PFTs的LAI計(jì)算得到。ki(i=0,…,4)是經(jīng)驗(yàn)系數(shù)，如表1所示。

2.3 Mann-Kendall 顯著性檢驗(yàn)

Mann-Kendall 檢驗(yàn)是一種檢驗(yàn)水文過(guò)程趨勢(shì)和其他相關(guān)物理變量的非參數(shù)方法，計(jì)算式為：

式中：x為獨(dú)立變量；n為數(shù)據(jù)總量，并且滿足：

統(tǒng)計(jì)量S接近正態(tài)分布，而統(tǒng)計(jì)量Z是一個(gè)標(biāo)準(zhǔn)正態(tài)變量，符合式（7）：

式中：i為群組中的個(gè)數(shù)；ti為對(duì)應(yīng)的數(shù)據(jù)。統(tǒng)計(jì)量Z被用來(lái)判斷是否存在顯著趨勢(shì)，正值代表上升趨勢(shì)，負(fù)值代表下降趨勢(shì)。

3 結(jié)果與分析

3.1 模型驗(yàn)證與評(píng)價(jià)

本文采用3 種不同的統(tǒng)計(jì)指標(biāo)來(lái)評(píng)價(jià)改進(jìn)后模型的精度，包括R2、RMSE和Bias。圖3（a）為用ERA 再分析氣象數(shù)據(jù)驅(qū)動(dòng)的PT-hybrid 算法模擬的2012—2014年LE值與49 個(gè)通量站LE實(shí)測(cè)值的對(duì)比結(jié)果。圖3（b）則是用同樣的氣象數(shù)據(jù)驅(qū)動(dòng)改進(jìn)后的算法得到的2012—2014年全球LE與實(shí)測(cè)LE的對(duì)比。驗(yàn)證結(jié)果表明，改進(jìn)后的算法能更準(zhǔn)確地模擬全球LE，R2從0.63 增加到0.74，RMSE從17.1 W/m2降低到11.7 W/m2，Bias從2.5 W/m2變?yōu)?3.2 W/m2。

圖3 基于49 個(gè)通量站點(diǎn)的年平均LE 驗(yàn)證 Fig.3 Comparison of annual LE observed from 49 flux towers and the corresponding estimated LE

為了進(jìn)一步探究算法改進(jìn)前后對(duì)于長(zhǎng)時(shí)間序列模擬的效果，本文選取時(shí)間序列較長(zhǎng)的6 個(gè)代表性站點(diǎn)（AU-Rig，CA-TP3，IT-Tor，DE-Akm，BE-Lon，US-Prr）的1982—2014年均LE值進(jìn)行重點(diǎn)分析。如圖4 所示，除了個(gè)別年份外，改進(jìn)的算法與原算法模擬的LE都呈現(xiàn)一致的趨勢(shì)性。總體而言，改進(jìn)的算法模擬結(jié)果要高于原算法的模擬結(jié)果。圖5 展示了上述6 個(gè)站點(diǎn)下墊面不同植被占比對(duì)應(yīng)的年際變化趨勢(shì)，由于氣候變化和人類活動(dòng)的影響，區(qū)域植被類型在時(shí)間尺度上會(huì)發(fā)生變化。因此，本文提出在算法中引入亞像元法，即在同一柵格里考慮多種動(dòng)態(tài)變化的植被類型，從而優(yōu)化了算法的模擬。

圖4 2 種算法在代表性站點(diǎn)的長(zhǎng)時(shí)間序列模擬結(jié)果對(duì)比 Fig.4 Comparison of results of the two algorithms in the long time series simulations at representative sites

圖5 代表性站點(diǎn)在研究時(shí)段內(nèi)下墊面不同植被功能類型的LAI 占比變化Fig.5 Example of a time series for LAI fraction of different PFTs simulated by LPJ-GUESS

3.2 全球LE 的時(shí)空變化特征分析

本文將研究區(qū)分為4 個(gè)時(shí)間段，分別為1982—1990、1991—2000、2001—2010年以及2011—2014年，重點(diǎn)從4 個(gè)不同年代的角度分析全球潛熱通量的變化。圖6 展示了4 個(gè)年代的全球LE分布以及MK趨勢(shì)性檢驗(yàn)得到的Z值分布。結(jié)果表明，LE顯著增加的地區(qū)主要分布在北緯70°～80°之間，尤其是亞歐大陸的東北部，以及非洲中部地區(qū)。而在美國(guó)南部，亞歐大陸中部，南美洲和澳洲的大部分地區(qū)，LE呈下降趨勢(shì)。從植被變化角度來(lái)看，美國(guó)南部地區(qū)和南美洲中部地區(qū)LE出現(xiàn)下降趨勢(shì)是由于EBF 的減少和GRA 的增加導(dǎo)致的；而在澳洲，由于近些年澳洲山火頻發(fā)[13]，導(dǎo)致森林面積大幅減少，造成LE的急劇降低；反觀亞歐大陸北部地區(qū)，由于全球氣候變暖，造成北極冰川融化，導(dǎo)致北極地區(qū)林線進(jìn)一步北移[14]，而這一趨勢(shì)在未來(lái)也會(huì)持續(xù)下去，因此亞歐大陸北部LE也會(huì)持續(xù)呈上升趨勢(shì)。

圖6 不同年代的年平均LE 分布圖及LE 的MK 趨勢(shì)檢驗(yàn)Z 值分布 Fig.6 distributions of annual average LE in different decades and Z-value distribution of LE trend

圖7 不同PFTs 的LAI 以及2 種算法模擬LE 的全球緯度剖面 Fig.7 Global latitudinal profiles of LAI of various PFTs simulated by LPJ-GUESS and LE simulated by two algorithms

本文從不同緯度分析了植被和LE的關(guān)系。圖7給出了改進(jìn)前后 2 種算法模擬得到的LE與LPJ-GUESS 模型模擬得到的各種植被類型在緯度上的分布。由圖7 可知，總?cè)~面積指數(shù)LAI和LE在緯度上呈現(xiàn)很強(qiáng)的一致性，分別在南北緯中緯度地區(qū)和赤道附近出現(xiàn)了峰值，而DBF、DNF、ENF 主要分布在南北緯40°～60°附近，EBF 則主要分布在赤道附近，與此同時(shí)，GRA 的分布與LE呈相反的趨勢(shì)，與先前的研究結(jié)論一致[15-16]，由此可以表明，LE與不同PFTs 在全球的分布密切相關(guān)。

4 討 論

本文引入LPJ-GUESS 模型，模擬得到全球植被的動(dòng)態(tài)變化過(guò)程，并對(duì)PT 算法進(jìn)行了改進(jìn)。模擬得到的1982—2014年全球LE值為30～95 W/m2，與Yao等[17]，C Jiménez 等[18]對(duì)于全球LE的模擬的總體范圍基本保持一致，但基于全球通量塔的地面實(shí)測(cè)LE驗(yàn)證，本文改進(jìn)的算法擁有更好的模擬效果，原因在于本文在原算法的基礎(chǔ)上做出2 點(diǎn)改進(jìn)：①采用LPJ-GUESS 模型來(lái)獲取全球植被的動(dòng)態(tài)變化過(guò)程作為PT 算法的輸入條件，替代靜態(tài)的下墊面植被類型；②研究采用亞像元法，對(duì)柵格植被進(jìn)行細(xì)化，提高了柵格植被的豐度，從而提高算法模擬的精度。

但由于不同植被內(nèi)在屬性不同，導(dǎo)致了其對(duì)于LE的影響程度也不同，本文只考慮了不同植被功能類型對(duì)于LE的定性影響，并未從定量的角度進(jìn)行探討，因此對(duì)于植被與潛熱通量關(guān)系的研究還有很大空間。此外，不同植被在不同季節(jié)的熱通量變化也有很大區(qū)別，對(duì)于年內(nèi)LE與植被的季節(jié)性研究也應(yīng)納入今后的考慮范圍。

5 結(jié) 論

1）通過(guò)與全球49 個(gè)通量塔實(shí)測(cè)數(shù)據(jù)的對(duì)比，驗(yàn)證了改進(jìn)算法的精度，改進(jìn)前后R2從0.63 增加到0.74，同時(shí)RMSE從17.1 W/m2降低到11.7 W/m2，從而優(yōu)化了算法的模擬效果。

2）采用MK 顯著性檢驗(yàn)方法分析了1982—2014年全球LE時(shí)空變化，其中亞歐大陸東北部以及非洲中部地區(qū)呈顯著的增加趨勢(shì)，而在亞歐大陸中部，南美洲和澳洲大分布地區(qū)，則呈顯著的下降趨勢(shì)。

3）從緯度角度分析不同植被分布與LE的關(guān)系，發(fā)現(xiàn)總的葉面積指數(shù)和LE在緯度分布上呈較強(qiáng)的一致性，同時(shí)出現(xiàn)3 個(gè)峰值，分別在南北緯40°～60°以及赤道附近。