張曉龍，沈 冰，黃領梅，權 全，莫淑紅，梁曉燕

(西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048)

地表凈輻射(Rn)是影響陸-氣能量交換和再分配過程的重要參數[1-2]，是氣候變化的重要驅動力，被廣泛應用于氣候監(jiān)測、天氣預報和農業(yè)氣象研究，是當前全球氣候變化研究領域的熱點之一[3]。同時，Rn也是構建各類生態(tài)模式的重要參數之一[4-5]，尤其是在生態(tài)系統(tǒng)的蒸散過程中起到非常重要的作用[6]。在全球氣候變化的背景下，研究Rn時空變化特征，對于了解過去和當前的氣候狀態(tài)以及對未來的氣候變化、蒸散發(fā)的估算、植物生長發(fā)育過程、生態(tài)系統(tǒng)生物量的形成與累積等研究具有重要意義。目前獲取地表凈輻射的方法主要是將地面實測數據按照某種規(guī)則進行空間插值，或通過經驗統(tǒng)計方法建立基于地面輻射實測數據和氣象資料估算地表輻射收支的經驗-統(tǒng)計模型，通過反演得到一定范圍內地表凈輻射[2, 7-8]。這兩種方法雖然在點尺度上估算精度較高，但由于氣象站數量有限且分布不均，地表覆被類型復雜，實時氣象參數(如氣溶膠，云的含量等)難以獲取，因此不能滿足大范圍凈輻射估算的需求[12-13]。遙感反演法應是區(qū)域由點及面較佳的技術手段。近年來，MODIS數據、Landsat系列數據等被廣泛應用于地表輻射收支的研究中[9-11]。

我國氣象系統(tǒng)僅有50個地表凈輻射觀測臺站，且均建于1993年以后。由于站點少，資料序列短，開展的相應研究也較少。青藏高原生態(tài)環(huán)境脆弱，被稱為氣候的“放大器”[12-13]，但常規(guī)氣象觀測臺站分布稀少且極不均勻，為缺資料地區(qū)空間范圍上凈輻射估算帶來極大阻礙，遠滿足不了科研和業(yè)務的需求。當前主流遙感反演產品使用地面站點較少，特別是在中國寒區(qū)、旱區(qū)[18]，數據精度還存在一定誤差。利用地面氣象站數據對遙感產品進行數據融合得到的再分析常規(guī)氣象要素的驅動數據集不僅可以滿足空間分布的要求，而且在數據精度上也可以大大提高，為進一步研究提供極大的便利。地表反照率決定了太陽輻射能量被地表反射與吸收的比例，是氣候系統(tǒng)的重要驅動因子之一[13]。地表發(fā)照率的分布與土地利用類型、植被狀況、水分狀況、積雪、地形等因素密切相關，在時間和空間上呈現高度異質性。因此，本文以黃河源區(qū)為例，利用ITPCAS氣象要素驅動數據集和GLASS地表反照率數據集，取用聯合國糧農組織(FAO)1998年推薦用于計算參考作物蒸散發(fā)Penman-Monteith公式估算缺資料地區(qū)長時間序列逐日Rn，在再分析數據集估算精度驗證的基礎上，用GIS空間分析技術與Mann-Kendall趨勢分析方法揭示其時間動態(tài)特征與空間演變規(guī)律，為進一步研究該地區(qū)氣候變化及陸地生態(tài)系統(tǒng)蒸散的影響提供技術支撐和理論依據。

1 研究區(qū)概況

研究區(qū)為龍羊峽水庫以上的黃河源區(qū)，位于32.2°～37.1°N，95.9°～103.4°E之間，屬于青藏高原東部(見圖1)。

圖1 研究區(qū)概況及氣象站點分布Fig.1 General situation of the study area and distribution of meteorological stations

黃河源區(qū)總體為西北高，東南低，阿尼瑪卿山為最高點，龍羊峽為最低處，海拔在2 508～6 253 m之間，有冰川、盆地、高山和峽谷等地貌。黃河源區(qū)屬大陸性高原氣候，干濕交替、冷熱分明，氣壓低，溫差大，輻射強，光熱條件較為豐富。受印度洋季風影響，降水集中在6～9月[5]，年平均降水234.6～839.2 mm，空間分布不均[12]。黃河源區(qū)3/4的面積分布著高寒植被，其中高寒草甸和高寒草原面積最大。黃河源區(qū)為多年凍土和季節(jié)性凍土的過渡區(qū)域，多年凍土分布下界平均氣溫在-2.5℃ ～ -3.5℃范圍內。由于氣候變暖，黃河源區(qū)近年來凍土發(fā)生顯著退化。

2 數據來源與方法

2.1 數據來源及預處理

再分析氣象數據集使用由中國科學院青藏高原研究所開發(fā)的ITPCAS(中國區(qū)域高時空分辨率地面氣象要素)驅動數據集(http://westdc.westgis.ac.cn/)[13]。該數據集以GEWEX-SRB 輻射資料、Princeton 再分析資料、TRMM 降水數據以及GLDAS數據為背景場，并融合中國氣象觀測數據制作而成。氣象要素包括降水量、溫度、比濕、風速、大氣壓力、向下短波和向下長波輻射等，時間分辨率為3 h，空間分辨率為0.1°。該數據集在黃河源區(qū)數據采樣點1 421個，可滿足黃河源區(qū)氣象要素空間分布的要求(見圖1)。本研究中所涉及的日尺度參數是通過3 h時段參數計算得到。為了分析Rn的時空變化特性，利用ArcGIS軟件CUBIC重采樣技術將0.1°數據集生成空間分辨率為1 km的數據集。氣象站數據源自國家氣象科學數據共享服務平臺(http://data.cma.cn/)。黃河源區(qū)再分析數據集數據采樣點和氣象臺站點分布見圖1。

地表反照率使用全球陸表參量GLASS數據產品，取值范圍0～1，被用于計算短波凈輻射，來源與北京師范大學(http://glass-product.bnu.edu.cn)。該產品2000-2012年反照率數據是基于MODIS數據利用AB(Angular Bin)算法開發(fā)的，其空間分辨率為1 km，時間分辨率為8天。Liu等利用地面實測數據和MODIS地表反照率數據MOD43對GLASS反照率數據產品進行對比驗證，結果顯示兩者具有很好的一致性，均方根誤差小于0.05[14]，所以該數據集可以直接使用。因數據源的時間限制，本研究計算得到2000-2012年白天地表反照率每一日的多年平均值，然后線性插值得到年內每天的反照率數據集。

DEM數據源自地理空間數據云(http://www.gscloud.cn/)中90 m 分辨率SRTM 產品。為了確保數據空間尺度的一致性，將90 m DEM數據利用ArcGIS軟件重采樣生成1 km分辨率的數據。

2.2 地表凈輻射估算方法

本研究采用FAO推薦的標準方法[15]，由凈短波輻射Rns和凈長波輻射Rnl之差得到Rn(MJ·m-2·d-1)：

Rn=Rns-Rnl

(1)

Rns=(1-α)Rs

(2)

(3)

Rso=(0.75+2z×10-5)Ra

(4)

(5)

cos(φ)cos(δ)sin(ωs)

(6)

式中：α為地表反照率；Rs和Rso為向下短波輻射和凈空向下短波輻射(MJ·m-2·d-1)；Ra為地外輻射(MJ·m-2·d-1)，n為日照時數，N為可照時數，z為海拔高度(m)；Tmax,k和Tmin,k分別為日最大、最小絕對溫度(K)；σ為Stefan-Boltzman常數(4.903×10-9MJ·m-2·d-1)；ea為實際水汽壓(kPa)；as=0.25，as+bs=0.75；Gsc為日輻射常數(0.082 MJ·m-2·min-1)，dr為日地距離，ωs為日落時角，φ為太陽時角，δ為太陽赤緯角。

需要指出的是，基于中國區(qū)域高時空分辨率地面氣象要素驅動數據集估算Rn_ITP時，Rs直接采用驅動數據集數據。而基于氣象臺數據估算Rn_m時，Rs采用式(5)計算得到。

2.3 地表凈輻射時空變化分析方法

采用Mann-Kendall趨勢分析方法對多站點和區(qū)域平均長時間序列的Rn進行時間趨勢變化檢驗和突變分析?；谥鹊腗ann-Kendall趨勢分析方法是一種非參數統(tǒng)計檢驗法，與傳統(tǒng)參數方法相比，樣本可不遵從特定分布，亦不受個別異常值干擾，可較客觀地反映樣本序列變化趨勢，該方法被廣泛應用于序列突變分析和趨勢檢驗。

利用線性回歸分析可得到每個柵格的Rn在一定時間序列的變化趨勢和強度[6]。本文通過IDL程序得到每個柵格多年變化趨勢b，進而得到黃河源區(qū)Rn年際變化率空間分布圖。計算公式為：

(7)

式中：b為變化趨勢，如果b<0則表示減少趨勢，反之b>0則表示增加趨勢；n為時間長度；i為年序號；xi為第i年的柵格值。

3 結果與分析

3.1 再分析數據集精度檢驗

黃河源區(qū)處于高寒區(qū)域，氣象站點稀少，且分布極不均勻。所以使用氣象站點進行空間插值的結果并不能較為真實地反應研究區(qū)實際情況。

本研究使用ITPCAS驅動數據集估算黃河源區(qū)Rn的時空變化。將再分析數據集代入式(1)～(6)得到的1979—2015年逐日Rn_ITP與氣象站數據代入式(1)～(6)得到的對應的逐日Rn_m。

為了定量評價Rn的估算結果精度，將7個典型氣象站的氣象數據估算的Rn_m與對應的1 km氣象要素數據集計算的Rn_ITP進行對比分析，Rn_ITP與Rn_m的散點圖見圖2。虛線為散點的擬合直線，對角線為1∶1等值線。

從圖2可以看出整體上Rn_ITP與Rn_m的相關系數R為0.98，其顯著性在P<0.01范圍內顯著，表明散點分布有一定的穩(wěn)定性，二者有顯著的相關性。大部分的點分布在1:1等值線上方，表明Rn_ITP較Rn_m略微偏大。

圖2 Rn_ITP與Rn_m的散點圖Fig.2 Scatter plots of Rn_ITPvaluesand Rn_mvalues

對逐日Rn的估算精度進行檢驗的評價指標為：均方根誤差(RMSE)、相關系數(R)和誤差，在各站點上計算結果見表1，其中區(qū)站號為國家氣象站統(tǒng)一編號，N為樣本數。從表1中可知，Rn_ITP與Rn_m的相關系數在0.98以上，且均通過顯著性檢驗(P<0.01)，各站點RMSE在1.29～1.77之間，值得注意的是，Rn_ITP均比Rn_m偏大12%～19%，主要原因可能是在氣象數據估算Rs時外輻射到達地表的比例參數選用默認值(as= 0.25和bs= 0.50)導致ITPCAS驅動數據集中Rs比氣象站基于日照時數計算的Rs偏高，也可能是氣象采樣點和氣象站并不完全重合，存在著坡度、海拔和坡向等差異。綜上，通過該數據集估算的Rn可以較好的反映區(qū)域變化規(guī)律，該再分析數據集在黃河源區(qū)有較好的適用性。

表1 再分析數據集估算Rn的結果驗證

3.2 黃河源區(qū)地表凈輻射的變化特征

3.2.1 地表凈輻射時間變化特征

精度檢驗后將ITPCAS驅動數據集和GLASS地表反照率數據集代入公式，得到1979―2015年黃河源區(qū)逐日Rn。對逐日數據進行處理，得到各月、季、年Rn數據。季節(jié)采用國內常用的氣象季節(jié)劃分法，即上年12月至當年2月為冬季、3月至5月為春季、6月至8月為秋季。各月Rn為月內每日平均得到，各季節(jié)Rn為相應月份Rn的平均值，年平均Rn為逐月數據平均得到。

圖3為1979-2015年黃河源區(qū)年、春季、夏季、秋季、冬季平均Rn變化曲線。1979-2015年黃河源區(qū)年、春季、夏季、秋季、冬季多年平均Rn分別為9.63、12.11、13.70、7.68和5.02 MJ·m-2·d-1；其最大值分別為10.00(1997年)、12.80(1980年)、14.59(1994年)、8.07(1997年)和5.38(1983年)MJ·m-2·d-1；最小值分別為9.12(2009年)、11.42(1989年)、12.63(2009年)、7.24(1988年)和4.76(2009年)MJ·m-2·d-1。

利用Mann-Kendall趨勢分析方法得到黃河源區(qū)年、四季Rn的變化趨勢和統(tǒng)計值，見表2。從表2可知，黃河源區(qū)年、四季Rn均呈下降趨勢，每10年變化幅度為-0.067 MJ·m-2·d-1；春、冬季變化幅度最大，每10年變化幅度為-0.075 MJ·m-2·d-1；秋季變化幅度最小，每10年變化幅度為-0.039 MJ·m-2·d-1；夏、秋季沒有通過顯著性檢驗。通過Mann-Kendall非參數檢驗分析研究區(qū)年、春季、夏季、秋季、冬季Rn突變點，得到結果：年Rn和冬季Rn序列存在突變點，分別為2004年和2003年左右；而其他季節(jié)Rn突變點不明顯。這可能與當地氣候變化和人類活動有關。2000年以后黃河源區(qū)呈現“變暖變濕”的趨勢，云量上升，降雨/雪增加，溫度也增加，植被NDVI也呈上升趨勢。另外三江源自然保護區(qū)于2000年成立，實施了一系列的生態(tài)保護與建設工程，青海省也實施了人工增雨、治理草原鼠害、退耕/牧還草和生態(tài)移民等 “綜合性”工程(2004年)。 這些都可能是年Rn和冬季Rn序列在2004年和2003年左右存在突變的原因。

圖3 1979-2015年研究區(qū)年、春季、夏季、秋季、冬季平均Rn變化曲線Fig.3 Curves of annual and seasonal mean Rn during 1979-2015 in study area

時間序列起始年份終止年份 樣本數量Test Z顯著性每10年變化幅度/(MJ·m-2·d-1)年Rn1979201537-2.68**-0.067春季Rn1979201537-2.31*-0.075夏季Rn1979201537-1.14-0.074秋季Rn1979201537-1.03-0.039冬季Rn1979201537-3.60***-0.075

注：***通過顯著性水平為0.001的檢驗；**通過顯著性水平為0.01的檢驗；*通過顯著性水平為0.05的檢驗。

圖4為黃河源區(qū)年內各月平均Rn的變化趨勢，可知，區(qū)域平均Rn在年內呈單峰型，1～4月上升迅速，5月達到最大值(14.31 MJ·m-2·d-1)，9月份以后Rn開始快速下降，12月份最小(3.34 MJ·m-2·d-1)。通過最大、最小值誤差線可知，1979-2015年Rn變化幅度夏季大于冬季，其中7月份Rn變化幅度最大，變化范圍為11.79～15.86 MJ·m-2·d-1。

3.2.2 地表凈輻射空間變化特征

對1979-2015年黃河源區(qū)逐柵格計算得到區(qū)域上年、春、夏、秋、冬季Rn和其變化幅度的空間分布及分類面積統(tǒng)計，結果見圖5。

從圖中可知，年和各季節(jié)Rn的空間分布相似，即不同季節(jié)對Rn的空間分布影響不大，最大值分布在扎陵、鄂陵湖及龍羊峽等大面積水域附近，最小值分布在中部的阿尼瑪卿山附近和共和盆地附近。年和各季節(jié)Rn的變化趨勢的空間分布也相似，研究區(qū)Rn以下降趨勢為主，主要分布在研究區(qū)西部山區(qū)、龍羊峽南部及共和盆地附近。研究區(qū)Rn上升趨勢主要分布在中部的阿尼瑪卿山附近和研究區(qū)東部的零星區(qū)域，年、春、夏、秋、冬季Rn上升趨勢分別占流域的8.4%、15.8%、4.6%、44.1%和0.0%。

3.2.3 海拔對地表凈輻射的影響

Rn的空間分布與海拔分布呈現一定的統(tǒng)計規(guī)律。將DEM數據重采樣到1 km柵格，然后將該DEM分成5類對多年平均的年、春、夏、秋、冬季Rn和其變化比例的進行統(tǒng)計分析，見表3。從表3中可知，Rn在4 000～4 500 m的區(qū)域內最大，而在4 500 m以上最小，這是因為4 500 m以上區(qū)域存在積雪和冰川，地表反照率高，從而Rn較小。

Rn的變化幅度中除了秋季在3 500～4 000 m區(qū)域內呈0.12%的上升趨勢外，其余皆呈現下降趨勢，年、春季和夏季Rn下降幅度最明顯的區(qū)域為4 500 m以上區(qū)域，秋季和冬季下降幅度最明顯的區(qū)域為3 000 m以下區(qū)域。

圖4 研究區(qū)Rn年內變化 Fig.4 Changes of monthly mean Rn in study area (The bars indicate the maximum and minimum monthly mean values)

DEM分類面積比例Rn/(MJ·m-2·d-1)年春季夏季秋季冬季每10年Rn變化比例/%年春季夏季秋季冬季3 000 m以下2.6%9.3611.8213.707.344.59-0.66-0.21-0.66-0.52-2.013 000～3 500 m16.0%9.3911.9513.297.474.85-0.78-0.67-0.79-0.31-1.723 500～4 000 m26.1%9.7312.3713.307.885.36-0.57-0.68-0.610.12-1.274 000～4 500 m37.4%9.8212.3213.927.965.07-0.79-0.80-0.84-0.13-1.644 500 m以上17.9%0.921.131.400.700.46-0.96-0.83-1.07-0.31-1.84

4 討 論

Rn的時空變化改變了區(qū)域地表天氣和氣候系統(tǒng)的多樣性、地表蒸散的速率及不同區(qū)域生態(tài)系統(tǒng)的結構和功能。劉新安等[1]對幾種推算地表凈輻射的方法進行了評價；陳征等[3]、葉晶等[16]利用遙感資料推算了個別地區(qū)的地表凈輻射，高揚子等[2]計算了全國699個站點的地表凈輻射，并分析了各分量的空間分布特征，對中國地表凈輻射做了較為系統(tǒng)的分析。

有研究表明，1980年以來全球Rn具有普遍的下降趨勢[17]， 國內近50年來站點平均Rn在年、季節(jié)均呈現較明顯的下降趨勢，每10年均降幅為0.74 W/m2(0.064 MJ·m-2·d-1)[2]，本研究得到與此結論一致。

影響特定區(qū)域Rn時空變化特征的因素主要有天文輻射、云、氣溶膠、溫度、土壤濕度和土壤植被等[2]。

有研究表明，云量不是造成中國區(qū)域短波輻射下降的主要因素，大氣中持續(xù)增多的氣溶膠是主導因素[18]。大氣氣溶膠主要來自于化石燃料的燃燒、工業(yè)生產等人為活動及火山爆發(fā)等自然現象。

自上世紀80年代以來，隨著經濟快速發(fā)展，化石燃料的大量燃燒造成大氣氣溶膠的快速增加，這可能是區(qū)域乃至全球Rn持續(xù)下降的主要原因。

有研究表明，青藏高原積雪變化具有確定的長期增加的趨勢[19]，從而導致的青藏高原冬季凈輻射的顯著降低，另外近年來由于環(huán)境污染和氣候變化等原因導致冬季的霧霾天氣偏低，這也可能造成冬季Rn顯著下降的原因之一。

圖5 1979-2015年研究區(qū)年和四季平均Rn以及對應的每10年Rn變化幅度的空間分布Fig.5 Distribution of mean Rn values and its change rates every 10 years during 1979-2015 in study area

5 結 論

本文以黃河源區(qū)為例，利用ITPCAS氣象要素驅動數據集和GLASS地表反照率數據集，基于FAO1998年推薦的Penman-Monteith公式估算缺資料地區(qū)長時間序列逐日Rn，在再分析數據集估算精度驗證的基礎上，用 GIS 空間分析技術與Mann-Kendall趨勢分析方法揭示其時間動態(tài)特征與空間演變規(guī)律。本文得到以下主要結論。

1) 各站點Rn_ITP與Rn_m的相關系數均在0.98以上，且均通過顯著性檢驗(P<0.01)，各站點RMSE在1.29～1.77之間，表明ITPCAS氣象要素驅動數據集在黃河源區(qū)有較好的適用性，通過該數據集估算的Rn可以較好地反映區(qū)域變化規(guī)律，該數據集可用于中國西北高寒山區(qū)缺資料地區(qū)。

2) 在1979-2015年期間年、春、夏、秋、冬季區(qū)域平均Rn的每10年下降幅度分別為-0.067、-0.075、-0.074、-0.039和-0.075 MJ·m-2·d-1，其中夏、秋季Rn趨勢沒有通過顯著性檢驗；年Rn的變化趨勢在2004年左右存在突變點，冬季Rn的變化趨勢在2003年左右存在突變點，其他季節(jié)突變點不明顯；區(qū)域平均Rn在年內呈單峰型，1～4月上升迅速，5月達到最大值(14.31 MJ·m-2·d-1)，9月份以后Rn開始快速下降，12月份最小(3.34 MJ·m-2·d-1)。

3) 年和各季節(jié)Rn值及變化趨勢的空間分布相似；年和各季節(jié)Rn值最大值分布在扎陵、鄂陵湖及龍羊峽等大面積水域附近，最小值分布在中部的阿尼瑪卿山附近和共和盆地附近；研究區(qū)91.6%的面積的年平均Rn呈下降趨勢，在季節(jié)尺度中，冬季Rn下降最明顯，幾乎100%的區(qū)域處于下降趨勢之中；秋季Rn下降最不明顯(約55.9%)。

4) 地形對Rn的空間分布產生一定的影響，Rn在4 000～4 500 m的區(qū)域內最大，而在4 500 m以上最小；Rn的變化幅度中除了秋季在3 500～4 000 m區(qū)域內呈0.12%的上升趨勢外，其余皆呈現下降趨勢；年、春季和夏季Rn下降幅度最明顯的區(qū)域為4 500 m以上區(qū)域，秋季和冬季下降幅度最明顯的區(qū)域為3 000 m以下區(qū)域。

本文對長時間序列的逐日逐柵格Rn進行計算，并進行時空趨勢變化分析，為無資料地區(qū)今后開展地表輻射研究、氣候變化預測、生態(tài)系統(tǒng)蒸散發(fā)過程和機理研究具有重要意義。由于缺少實測的Rn觀測數據，在計算精度上存在一定的不確定性。另外本文未定量分析造成黃河源區(qū)Rn下降的主要驅動因子，Rn的影響過程與機理需在今后研究工作中進一步探討。