2000—2021 年青藏高原地區(qū)地表凈輻射的時空變化*

崔玉祥 胡斯勒圖 李同文 姬大彬 張 顥 施建成

1（中國科學(xué)院空天信息創(chuàng)新研究院 北京 100094）

2（中山大學(xué)測繪科學(xué)與技術(shù)學(xué)院 珠海 519080）

3（中國科學(xué)院國家空間科學(xué)中心 北京 100190）

0 引言

地表凈輻射是指地球表面的入射和出射輻射之間的差值，包括長波和短波輻射，是衡量地表輻射收支總量的一個重要參數(shù)。從地表能量平衡的角度看，地表凈輻射可以轉(zhuǎn)化為熱通量，包括土壤熱通量、感熱通量和潛熱通量等，地表凈輻射是影響整個地氣能量交換和再分配的重要參數(shù)[1]。同時，作為研究地氣間相互作用的關(guān)鍵因素，地表凈輻射是地表蒸散、大氣熱量和水循環(huán)的主要驅(qū)動力，其時空變化直接影響著水分和能量交換[2-5]，用于氣候監(jiān)測、天氣預(yù)報和農(nóng)業(yè)氣象等領(lǐng)域。

青藏高原地處亞洲大陸中心地帶，平均海拔超過4 km，是地球上海拔最高且地形復(fù)雜的高原之一，素有地球第三極之稱[6]。青藏高原與大氣圈、生物圈和水圈等圈層耦合作用，對高原及其周邊地區(qū)氣候格局、亞洲季風(fēng)形成，乃至北半球大氣環(huán)流、全球的氣候變化和極端天氣形成有巨大影響[7-9]。因此，分析研究青藏高原地區(qū)地表凈輻射的時空變化具有重大意義。

對青藏高原地區(qū)地表凈輻射的時空變化特征已有很多研究，并取得了一定成果[10-13]。Ciren 等[14]采用青藏高原羊八井地區(qū)2009—2010 年期間輻射觀測數(shù)據(jù)，研究發(fā)現(xiàn)羊八井地區(qū)地表凈輻射最高值出現(xiàn)在6 月，最低值出現(xiàn)在12 月。Li 等[15]基于獅泉河地氣交換綜合觀測站和獅泉河氣象站輻射觀測數(shù)據(jù)，研究了阿里地區(qū)凈輻射通量的日循環(huán)、日變化、季節(jié)變化和年際變化特征。結(jié)果表明，阿里地區(qū)凈輻射平均日變化峰值出現(xiàn)在14:30 LT—15:00 LT，凈輻射月均值最大值出現(xiàn)在8 月，且近20 年來凈輻射的平均日變化峰值和月均值最大值出現(xiàn)時間延后現(xiàn)象。Zhao 等[16]利用青藏高原羊八井觀測站2020—2021 年的太陽短波輻射數(shù)據(jù)和地氣長波輻射數(shù)據(jù)，發(fā)現(xiàn)羊八井地區(qū)地表凈輻射呈現(xiàn)顯著的單峰月際變化特征，月平均值6 月最高，12 月最低。以上研究為青藏高原地區(qū)地表凈輻射的研究奠定了重要基礎(chǔ)，但是大部分地表凈輻射變化規(guī)律僅從時間或者空間一個維度分析，沒有在時空維度綜合考慮。另外，大多數(shù)研究只在站點位置或者某一局部區(qū)域進行研究，涉及青藏高原全域的較少。衛(wèi)星遙感數(shù)據(jù)相較于地基站點數(shù)據(jù)，擁有大面積、長時序等顯著優(yōu)勢，特別有利于大范圍長時間序列的地表凈輻射估算問題的研究[17,18]。

本文利用CERES 獲取青藏高原地區(qū)2000 年3 月1 日至2022 年2 月28 日的逐日地表輻射各分量數(shù)據(jù)，計算地表凈輻射數(shù)據(jù)，并結(jié)合經(jīng)驗正交函數(shù)分析法（Empirical Orthogonal Function，EOF）、Theil-Sen Median 趨勢分析和Mann-Kendall 分析法，深入探究青藏高原凈輻射通量的時空變化特征，以期為青藏高原地區(qū)地氣間的能量交換和再分配探究以及長期氣候變化研究等提供科學(xué)依據(jù)。

1 研究區(qū)域

青藏高原是世界上海拔最高、面積最大的高原，位于中國西南部。圖1 表明研究區(qū)** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/東部地勢較低，西部地勢較高，平均海拔超過4000 m。總體而言，該地區(qū)地勢高，大氣稀薄，太陽輻射能量更容易到達地面，是全球輻射最強的地區(qū)之一。本文在基于站點數(shù)據(jù)驗證CERES 衛(wèi)星數(shù)據(jù)有效性基礎(chǔ)上，結(jié)合多種分析方法探究青藏高原地區(qū)地表凈輻射的時空變化特征。

圖1 青藏高原地區(qū)高程及地基站點分布Fig. 1 Elevation map of Qinghai-Tibet Plateau

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源及預(yù)處理

研究采用的數(shù)據(jù)為2000 年3 月1 日至2022 年2 月28 日CERES 的SYN1 deg-Level 3 數(shù)據(jù)集*** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/逐日地表輻射各分量數(shù)據(jù)--地表上行短波輻射、地表下行短波輻射、地表上行長波輻射、地表下行長波輻射，由公式1 計算得到地表凈輻射，空間分辨率為1°×1°，數(shù)據(jù)單位為W·m-2。

使用的地基站點數(shù)據(jù)來自于中國氣象數(shù)據(jù)網(wǎng)站**** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/，基于如表1 所示的青藏高原地區(qū)所有地基站點（拉薩、西寧、那曲、昌都、玉樹），獲取2015 年12 月31 日18:00 LT 到2016 年12 月31 日17:00 LT 逐小時地表下行短波輻射，用于評估CERES 輻射參數(shù)的數(shù)據(jù)質(zhì)量，地表凈輻射

表1 青藏高原地區(qū)拉薩、西寧、那曲、昌都和玉樹5 個地基站點參數(shù)Table 1 Parameters of five ground stations in Qinghai-Tibet Plateau region,including Lhasa, Xining, Naqu, Changdu and Yushu

式中，Rsd為地表下行短波輻射；Rld為地表下行長波輻射；Rsu為地表上行短波輻射；Rlu為地表上行長波輻射，單位為W·m-2。

2.2 有效性驗證

衛(wèi)星數(shù)據(jù)比地基站點觀測數(shù)據(jù)空間覆蓋能力強，更利于反映地表凈輻射的時空分布和分配過程，但是地基站點的觀測具有更高的精度。因此，利用5 個地基站點數(shù)據(jù)對衛(wèi)星數(shù)據(jù)進行精度驗證，相關(guān)數(shù)據(jù)采用2015—2016 年共8784 h 地表下行輻射數(shù)據(jù)，以地基站點值為參考值，分別計算決定系數(shù)R2和平均絕對誤差E（單位W·m-2），即

其中，N為樣本數(shù)量，X和X分別為地基站點觀測真值及其平均值；Y和Y分別為CERES 衛(wèi)星數(shù)值及其平均值。圖2 給出了CERES 衛(wèi)星數(shù)據(jù)有效性驗證結(jié)果，五個站點的R2分別為0.92，0.87，0.88，0.85 和0.87，同時，E分別為44.41，55.51，45.03，56.80 和53.64 W·m-2。上述結(jié)果表明，CERES 數(shù)據(jù)相對穩(wěn)定且有效性較高，能夠滿足研究需求。

圖2 拉薩（a）、西寧（b）、那曲（c）、昌都（d）和玉樹（e）站點CERES地表下行短波輻射和地基站點數(shù)據(jù)驗證Fig. 2 Validation of CERES surface downward shortwave radiation data with ground-based station data at sites in Lhasa （a）, Xining （b）, Naqu （c）, Changdu （d） and Yushu （e）

2.3 研究方法

2.3.1 趨勢分析法

Theil-Sen Median 趨勢分析[19,20]，是一種穩(wěn)健性強且不依賴參數(shù)的趨勢計算方法，該算法擺脫了測量值服從正態(tài)分布的缺點，并且在缺失值和異常值影響下依然有效，因而在長時序數(shù)據(jù)趨勢分析中得到了廣泛應(yīng)用。具體計算過程如下：

式中，Median 為中位數(shù)運算，Rni和Rnj分別代表第i和j年的地表凈輻射，i和j分別為第i和j年，且j＞i；S表征地表凈輻射的變化趨勢。當S ＞0，表明該研究時段內(nèi)地表凈輻射有加強趨勢；當S=0，表明該研究時段內(nèi)地表凈輻射無變化趨勢；當S ＜0，則表明該研究時段內(nèi)地表凈輻射有減弱趨勢。

2.3.2 Mann-Kendall 檢驗

Mann-Kendall 顯著性檢驗法[21,22]，是一種不受異常值影響且操作簡單的檢驗長時序變化趨勢顯著性的方法，已經(jīng)在水文、氣象等領(lǐng)域取得成功應(yīng)用。本文采用該方法來檢驗地表凈輻射變化的顯著性。具體算法如下：

式中，1 ≤i ＜j≤n，n為時間序列長度；fji為Rnj-Rni的差值函數(shù)；Svar為S的方差，Z反映地表凈輻射變化趨勢的顯著性。當|Z|≥2.58時，表示地表凈輻射變化趨勢已通過顯著性水平α=0.01的顯著性檢驗；當1.96＜|Z|＜2.58時，則說明地表凈輻射變化趨勢通過了顯著性水平α=0.05的顯著性檢驗；當|Z|≤1.96表示地表凈輻射變化趨勢未能通過顯著性水平α=0.05的顯著性檢驗。

Theil-Sen Median 趨勢分析對異常值不敏感，可用于檢測數(shù)據(jù)中是否存在某種趨勢，而Mann-Kendall 顯著性檢驗法可以驗證被檢驗的時間序列的趨勢性是否通過某置信度下的顯著性檢驗。結(jié)合Theil-Sen Median 趨勢分析與Mann-Kendall 顯著性檢測（見表2[22]），即可更有效反應(yīng)地表凈輻射的變化趨勢。

表2 地表凈輻射變化特征劃分標準Table 2 Criteria for classifying changes in surface net radiation

2.3.3 EOF 分析

EOF 分析[23-25]是一種常用的數(shù)據(jù)分析方法，常用于對海氣耦合系統(tǒng)、氣象、氣候等領(lǐng)域中的大規(guī)模環(huán)境數(shù)據(jù)進行降維、分解和重構(gòu)。其不僅可以降低數(shù)據(jù)維度，并且能盡可能地保留原始數(shù)據(jù)的空間和時間變化特征。特征向量所對應(yīng)的空間樣本可以展現(xiàn)其空間分布特點；而主成分則表現(xiàn)時間變化，揭示了隨時間變化的不同空間模態(tài)的重要性權(quán)重。

對青藏高原地區(qū)2000 年3 月到2022 年2 月共264 個月逐像素（共299 個像素點）地表凈輻射處理成地表凈輻射距平值Xm×n，基于EOF 分析求解特征根及相應(yīng)特征向量，計算各空間模態(tài)的方差貢獻率并進行North 檢驗，具體計算公式如下：

式中，n=264，m=299，Xm×n為地表凈輻射距平值；Cm×m為協(xié)方差矩陣；λ1，λ2,...,λm為協(xié)方差矩陣Cm×m對應(yīng)的特征根；Vm×m為特征向量；Em×m為m×m的對角矩陣，Pm×n為時間系數(shù)，A為方差貢獻率，R為North 檢驗容許的誤差范圍。

3 結(jié)果與討論

3.1 年平均地表凈輻射的分布特征

圖3 給出了青藏高原地區(qū)逐年平均地表凈輻射空間分布，從空間分布來看，青藏高原年平均地表凈輻射的空間分布特征首先表現(xiàn)為緯度水平條帶特征，即年平均地表凈輻射沿緯度從南到北逐漸降低，最高值出現(xiàn)在喜馬拉雅山脈和橫斷山脈一線。其次，從經(jīng)度方向來看，年平均地表凈輻射從西向東大致呈現(xiàn)先減弱，后增強，再減弱，再增強的高—低—高—低—高模式，這主要與青藏高原的地勢有關(guān)[26]，青藏高原西部有著平均海拔6000 m 以上的喜馬拉雅山脈，唐古拉山脈橫臥青藏高原中部，而橫斷山脈盤亙于青藏高原東部。

從時間變化來看，年平均地表凈輻射2000—2016 年大都在100 W·m-2以上，2017—2021 也大都在90 W·m-2以上；青藏高原中部地區(qū)2000—2016 年處于90～100 W·m-2之間，2017—2021 處于80～90 W·m-2之間，最低值為塔里木盆地，常年處于50～60 W·m-2之間。其中，值得注意的是，青藏高原大部分地區(qū)在2000—2016 年年平均地表凈輻射基本保持平穩(wěn)狀態(tài)，而在2016—2017 年間發(fā)生突變，2017—2021 年年平均地表凈輻射也基本維持不變。

3.2 結(jié)合趨勢分析法與M-K 檢驗法的地表凈輻射趨勢分析

圖4 給出了2000—2021 年青藏高原地區(qū)地表凈輻射變化分布，分析可知，2000—2021 年青藏高原地表凈輻射整體處于基本不變和不顯著減少狀態(tài)。其中，帕米爾高原、昆侖山、藏北高原、柴達木盆地、藏東南等地區(qū)地表凈輻射基本不變；岡底斯山脈、唐古拉山脈、昆侖山脈、青海湖以及青藏高原中東部均呈現(xiàn)不顯著減少趨勢；喜馬拉雅山脈附近顯著減少；藏北高原附近出現(xiàn)零星不顯著增加區(qū)域。結(jié)合圖5 地表輻射各分量變化趨勢，可以看到喜馬拉雅山脈附近地表下行輻射呈穩(wěn)定不變或者不顯著減少趨勢，而地表上行輻射卻出現(xiàn)不顯著增加，甚至顯著增加趨勢，因此，此地地表凈輻射呈顯著減少和不顯著減少狀態(tài)。此外，青藏高原地表上行長波輻射全域呈現(xiàn)基本不顯著增加趨勢，這在一定程度上表明了青藏高原大部分地區(qū)地表溫度總體呈增長趨勢，響應(yīng)了全球變暖的氣候特征，這與Gu 等[6]探究青藏高原地表輻射通量的氣候特征結(jié)論高度吻合。

圖4 2000—2021 年青藏高原地區(qū)地表凈輻射變化分布Fig. 4 Distribution of surface net radiation changes in the Qinghai-Tibet Plateau region from 2000 to 2021

圖6 給出了2000—2021 年青藏高原地區(qū)地表凈輻射季節(jié)、年均變化趨勢以及云覆蓋面積分數(shù)變化趨勢。由圖6 分析可知，從年內(nèi)變化來看（見圖6 a），青藏高原地表凈輻射有明顯季節(jié)變化特征。2022 年來青藏高原地表凈輻射均呈現(xiàn)出夏季＞春季＞秋季＞冬季的規(guī)律。原因可能在于夏季青藏高原地區(qū)太陽高度角大，太陽總輻射相應(yīng)較大，對地表凈輻射的增量遠大于因大氣中水氣含量較高以及云量相對較高對地表直接輻射而產(chǎn)生的減少量[27]。這也與Gao 等[1]關(guān)于青藏區(qū)地表凈輻射各個季節(jié)總體趨勢保持一致。

圖6 2000—2021 年青藏高原地區(qū)地表凈輻射季節(jié)趨勢（a）、地表凈輻射年均變化趨勢（b）和云覆蓋面積分數(shù)變化趨勢（c）Fig. 6 Seasonal trends （a）, annual changes （b）, and trends in cloud area fraction （c） for surface net radiation in the Qinghai-Tibet Plateau region from 2000 to 2021

從圖6（b） 所示的年際變化來看，2000—2021 年青藏高原地表凈輻射在2016 年發(fā)生顯著變化，突降約6 W·m-2。表現(xiàn)為2000—2016 年和2017—2021年，兩個時間段內(nèi)地表凈輻射各季均值變動幅度不大，但是2017 年各季節(jié)地表凈輻射均明顯低于2016 年數(shù)值。如圖6（c） 所示，查閱資料后發(fā)現(xiàn)青藏高原地區(qū)2000—2016 年云覆蓋度年均值基本維持在47.01%， 而在2017—2020 年，這一數(shù)據(jù)達到55.87%，提升約18.75%，地表凈輻射與云覆蓋分數(shù)突變一致性較強，進一步利用皮爾遜相關(guān)性檢驗發(fā)現(xiàn)：兩者相關(guān)性系數(shù)為-0.90（p＜0.001），具有顯著負相關(guān)性，猜測可能是云覆蓋度增加會阻擋更多太陽輻射，使得地表收到的太陽輻射減少，最終導(dǎo)致地表凈輻射下降。

3.3 EOF 方法地表凈輻射通量變化分析

表3 給出了前5 個特征向量的特征根、方差貢獻率以及特征根誤差范圍，從表3 可以看出，第一空間模態(tài)貢獻最大，為75.55%，第二方差貢獻急劇減少，僅為7.58%，其后，方差貢獻率均較小，前2 個特征向量的累積方差貢獻率已經(jīng)達到83.13%，并且都通過了North 顯著性檢驗，因此前2 個特征向量基本可以反映青藏高原2000—2021 年地表凈輻射的分布類型。

表3 青藏高原地區(qū)地表凈輻射時空場EOF 分解前5 個模態(tài)特征值與方差貢獻率Table 3 Characteristic values and variances contribution rate of the first five EOF modes for the spatial-temporal pattern of surface net radiation in the Qinghai-Tibet Plateau region

第1 模態(tài)特征向量的方差貢獻率為75.55%，遠高于其他模態(tài)的方差貢獻率，可以反映青藏高原地表凈輻射分布最重要的典型場，代表區(qū)域分布的平均特征。從圖7（a）所示的第一特征向量空間分布特征來分析，第1 模態(tài)在青藏高原地區(qū)的空間系數(shù)全部為正值，這表明2000—2022 年青藏高原地表凈輻射變化趨勢具有整體一致性，即青藏高原地區(qū)地表凈輻射均呈現(xiàn)出普遍較強或者普遍較弱的特征。同時，第1 模態(tài)系數(shù)大致呈以喀喇昆侖山脈為低值中心向四周擴大的趨勢，高值區(qū)位于喀喇昆侖山脈、喜馬拉雅山脈東南部、巴顏克拉山脈、橫斷山脈附近，表明這些區(qū)域地表凈輻射變化大。在圖7（b）所示的第1 模態(tài)時間系數(shù)中，110 個月為正值，154 個月為負值，第一時間系數(shù)均在零值上下周期浮動，且呈現(xiàn)準正弦振動，年周期變化明顯，這主要是受到太陽高度角的周期變化影響。峰值主要出現(xiàn)在每年夏季的7 月和8 月，谷值出現(xiàn)在每年冬季的2 月和3 月，恰好與青藏高原夏冬季節(jié)氣候特征吻合。

圖7 2000—2021 年青藏高原地區(qū)地表凈輻射第1 模態(tài)（a）特征向量分布及對應(yīng)月時間分布（b） ，第2 模態(tài)（c）特征向量分布及對應(yīng)月時間分布（d）Fig. 7 Distribution of the first mode （a） and corresponding monthly time distribution （b） of surface net radiation, distribution of the second mode （c） and corresponding monthly time distribution （d）in the Qinghai-Tibet Plateau region from 2000 to 2021

第2 模態(tài)特征向量的方差貢獻率為7.58%，也屬于青藏高原地區(qū)較為典型的地表凈輻射表現(xiàn)形式。從圖7（c）可以看出，這種分布格局大致以喀喇昆侖山脈和喜馬拉雅山為高值中心向四周減小，青藏高原中部大范圍為正值，而在興都庫什山脈、祁連山脈和喜馬拉雅山脈南麓出現(xiàn)負值中心，大體呈現(xiàn)西南-東北反向分布模式。除喜馬拉雅山脈外，特征向量值從西南向東北依次減小，反映出青藏高原地區(qū)地表凈輻射變化也是從西南向東北依次遞減。由圖7（d）分析表明，第2 模態(tài)對應(yīng)的時間系數(shù)也有一定的年周期變化，這可能是地表綜合作用的結(jié)果。其中，時間系數(shù)為負值有123 個月，地表凈輻射位相與圖7（c） 相反，空間分布呈現(xiàn)喀喇昆侖山脈和喜馬拉雅山低值、而周圍高值的情形，典型時間為2004 年8 月。時間系數(shù)為正的月份有141 個月，呈中間高四周低分布，典型時間為2002 年12 月。

4 結(jié)論

青藏高原高原氣候復(fù)雜且大氣圈、生物圈和水圈等圈層耦合作用，導(dǎo)致其地表凈輻射的時空分布和變化復(fù)雜多樣。基于2000 年3 月至2022 年2 月共22 年的CERES 衛(wèi)星數(shù)據(jù)，采用Theil-Sen Median 趨勢分析、Mann-Kendall 檢驗以及EOF 分析多種方法，分析22 年間青藏高原地表輻射能量平衡的時空變化規(guī)律，得出結(jié)論如下。

（1）青藏高原地表凈輻射通量總體呈現(xiàn)南高北低的分布特征，最高值在喜馬拉雅山脈和橫斷山脈一線，最低值在柴達木盆地。其中，青藏高原大部分地區(qū)在2000—2016 年地表凈輻射通量明顯更強，而2017—2021 年地表凈輻射通量明顯降低。

（2）Theil-Sen Median 趨勢分析和Mann-Kendall檢驗表明，青藏高原地表凈輻射通量整體處于基本不變和不顯著減少狀態(tài)，但是喜馬拉雅山脈附近地區(qū)的地表凈輻射通量顯著減少。其他年份地表凈輻射通量大體處于平穩(wěn)狀態(tài)，2016—2017 年，出現(xiàn)近6 W·m-2的突變，這可能是對云覆蓋分數(shù)增加的響應(yīng)。

（3）EOF 分解結(jié)果表明，青藏高原地表凈輻射通量場主要有兩種模態(tài)。第1 模態(tài)方差貢獻率為75.55%，可以代表地表凈輻射通量變化的總體空間特征，地表凈輻射通量變化趨勢具有高度一致性，青藏高原南部地表凈輻射通量波動程度遠高于其北部，對應(yīng)的時間系數(shù)顯示，年周期變化明顯，第一峰值出現(xiàn)在2018 年8 月，第一谷值出現(xiàn)在2001 年3 月；第2 模態(tài)大體呈現(xiàn)西南-東北反向分布模式，除喜馬拉雅山脈外，特征向量值從西南向東北依次減小，其時間系數(shù)也有一定的年周期變化，這可能是地表綜合作用的結(jié)果。