基于Landsat 8 的川西復雜地形太陽輻射估算研究

郭豫賓 , 張船紅,2

（1. 四川建筑職業技術學院現代測繪技術中心, 德陽 618000；2. 成都理工大學地球物理學院, 成都 610059）

引言

太陽輻射作為陸地生態系統中物質和能量交換過程的重要組成部分，通過與地表大氣圈、生物圈和水圈的交互作用實現地表能量交換，并促使地表輻射能量平衡，對地表自然環境和人類社會活動有著十分重要的作用[1-3]。

隨著地球科學研究的不斷深入，如何獲取區域地表太陽輻射數據，并認識其時空分布特征變得日趨重要[4]。山地因受地形、植被、坡度、坡向及地形遮蔽等因素影響，區域地表所接收的連續分布的太陽輻射量較難測量。由于山區地面觀測站點稀少、涵蓋面小且分布不均，空間分布相比平坦地表呈現更高的時空異質性，單一利用山區有限地面觀測站點通過插值來獲取區域地表太陽輻射有很大的不確定性，如何獲得區域的連續準確的太陽輻射數據有待研究[5-6]。

對于山地太陽輻射估算，前人提出了眾多方法理論。如傅抱璞[7]提出坡地太陽輻射的計算方法。而針對山地丘陵特點，唐曉萍等[8]利用重慶地區34個氣象臺站的日照時數、平均氣溫和沙坪壩輻射站太陽總輻射數據，重點分析了該地區氣溫、總輻射的時空變化特征。陳中鈺等[9]利用氣象站點逐月日照數據和逐月太陽總輻射數據，用氣候學模擬計算了四川144個氣象站逐月太陽總輻射值。鐘燕川等[10]則利用氣象站常規數據和數字高程模型（Digital Elevation Model，DEM），考慮坡度、坡向等地形因子建立氣候經驗統計模型和分布式模型，計算了四川省太陽總輻射并分析其時空分布情況。但上述工作是借助氣象數據、地形等原始資料，導致計算量大而繁瑣[4]。

近年來，伴隨遙感及計算機技術的發展，利用數字高程模型來計算地形坡度、坡向和地形遮蔽度等，為其可視化表達提供了便利，借助遙感影像結合氣象站點觀測資料和數字高程模型來模擬復雜地區地表太陽輻射成為可能。本文在借鑒已有太陽輻射模型的基礎上，充分考慮山地的海拔、日照時數、坡向、坡度、天空可見因子、大氣透過率和地表表觀反照率等因素，利用Landsat 8 遙感影像、氣象輻射站點常規觀測資料和DEM 數據，估算川西山地地表瞬時太陽輻射，研究其不同季節的空間分布差異以及影響因子與山地太陽輻射空間分布的關系，以期為提升山地地表太陽輻射研究水平提供科技支持。

1 數據源及研究區概況

1.1 數據源

文中所用數據主要包含以下幾種：（1）遙感影像。主要是Landsat 8 影像數據，行列號為131/038。影像成像質量佳，云量較少，日期分別為2016.05.12、2016.07.15、2016.10.03 和2016.12.22。用于計算區域地表反照率，植被指數。（2）數字高程模型（DEM）?？臻g分辨率30 m，用于計算區域坡度和坡向。（3）常規氣象輻射站點觀測數據。研究區內分布有馬爾康、紅原和小金三個地面氣象站點，獲取了包括氣溫日值、相對濕度日值、日照時數等觀測數據，主要用于計算研究區大氣透過率。其中，紅原站可提供正午時刻輻射強度和日總太陽輻射量，用于山地太陽輻射估算結果的精度驗證（表1）。

表1 數據源及特征

1.2 研究區概況

文中所選區域位于甘孜自治州與阿壩自治州交界，以Landsat 8 影像（行列號131/038）確定區域作為研究范圍。該區域地形地貌變化明顯，地形以山區為主，地處橫斷山構造帶，境內山峰林立，山丘河谷相間，平均海拔3000 m 以上，為川西高原一部。區域內主要有岷山、巴顏喀拉山等。

研究區主要氣候特征是河谷干暖，山地冷濕，光照豐富。山原地帶夏季溫涼，冬春寒冷，干濕季明顯，年平均氣溫5.6～8.9℃。高山峽谷地帶，氣候立體變化明顯，呈明顯的垂直性差異。區域內分布的三個氣象地面觀測站海拔分別是2664.4 m、2438 m 和3491.6 m（圖1）。

圖1 研究區位置及其高程分布

2 瞬時太陽輻射估算原理

地表太陽輻射是太陽照射到地球大氣層頂，經過大氣層的吸收、反射、散射后達地球表面的那部分太陽輻射[11]。大氣作用大大削弱了到達地面的太陽輻射，山區受緯度、海拔、天氣及日照時數的影響，地表瞬時太陽輻射估算更為困難。為準確估算山地地表太陽輻射，在前人研究的基礎上，依據復雜地形地表太陽輻射估算思路，將研究區正午時刻地表接收的瞬時太陽輻射分為太陽直接輻射、天空散射輻射和周圍地形反射輻射，分別進行計算[4,12]。最后，通過求和的方式估算出研究區特定時間的地表瞬時太陽輻射。理論公式如下：

式中：Id是入射到地表的太陽直接輻射，Ias為散射輻射，Ir為周圍地形反射輻射[13]。

2.1 太陽直接輻射

未考慮地形起伏狀況下，太陽直接輻射是太陽以平行光線的形式直接投射到地面上的輻射?？紤]山區大氣作用及地形坡面的影響，太陽直接輻射計算公式[14]如下：

式中：I0為太陽常數，當需要具體某一日期的太陽常數時，需根據日地距離系數 δ進行訂正，一般取值為1367 W/m2；τb為研究區直接輻射的大氣透過率；δ為日地距離系數；θ為太陽入射角（rad）。當 cosθ＜0時，太陽直接輻射被完全遮蔽，上述結果取0 值。

2.1.1 直射輻射大氣透過率計算

直接輻射大氣透過率受山地地形、天氣及海拔影響，直接獲取難度較大。目前普遍采用Modis 數據進行估算，本文受數據源限制，采用氣象站點觀測數據來估算研究區日大氣透過率進行代替計算。具體模型[15]為：

式中：ss和sr表示站點日照時數，取值由氣象站點觀測數據為準；τ0,dry為晴朗大氣狀況下，垂直地表方向的瞬時輻射透過率，最佳取值為0.87[13]；Pa為隨海拔變化的研究區大氣壓（Pa），Pa=101.3exp(-Ea/8200)，Ea為當地高程值，文中以區域平均海拔值代入計算；P0為標準大氣壓，P0=101.3 kPa；e為近表面水汽壓；根據氣象學公式，近表面水汽壓由氣象站點觀測的氣溫及相對濕度數據計算得到，因篇幅有限，在此不再贅述；α表示氣壓對透過率的影響，α=-6.1×10-5；mβ表示太陽天頂角為 β時的光學空氣質量，mβ=1/cosβ；太陽天頂角與太陽高度角互為余角[15-16]。

2.1.2 日地距離系數計算

日地距離又稱太陽距離，以日地平均距離為標準進行修正后得到的值。當地日地距離值與日地平均距離的比值為日地距離系數，記為δ。計算方法如下：

式中：年角γ=2π×(n-1)/365（單位：rad），n為日期在一年中的日序數，即1～365，閏年為366 天。

2.1.3 太陽入射角計算

太陽直線光線與地表坡面法線形成的夾角即坡面上的實際入射角。傾斜面太陽入射角 θ的計算方法與水平面不同，具體計算公式[17]為：

式中：S、A分別是研究區坡度、坡向，由DEM 計算得到；β為太陽天頂角，φ0為太陽方位角。本文計算正午時刻太陽高度角是依據研究區中心地理緯度、太陽赤緯及時角進行計算。太陽方位角則由太陽高度角、地理緯度和太陽赤緯計算得到，計算已較為成熟，這里不再贅述。

2.2 太陽散射輻射

太陽散射輻射是太陽輻射通過大氣時，受到大氣中氣體、塵埃、氣溶膠等的作用，從天空的各個角度到達地表的一部分太陽輻射，又稱天空散射輻射。計算模型[14]見下式：

式中：I0，β，δ上文已介紹；τd為散射輻射透過率。天空晴朗無云時可認為是一個均勻介質，它與直射輻射透過率 τb為 線性關系[18]，即：

復雜地形天空散射輻射受天空可視因子Vd影響，該因子是從某特定位置查看天空時，對整個可見或被遮擋天空的柵格表達，與地形坡度、坡向密切相關。山區天空散射輻射是在平坦地表天空散射輻射的基礎上，利用天空可視因子Vd進行修正。天空可視因子Vd估算公式如下：

式中：S、A分別為研究區坡度、坡向；Hφ為水平張角，即天頂向下每個方向相對于水平面的夾角[17]。研究區內山形復雜，Hφ取值由水平方向至天頂方向，即π/2；方位角 φ取16 方位角值代入計算，單位為rad。

2.3 周圍地形反射輻射

周圍地形反射輻射與研究區坡度及地表反照率相關。本文采用簡化近似計算，具體方法[19]為：

式中：ρ表示坡面地表平均反照率。山區坡面平均反照率ρ是在Landsat 8 影像進行大氣校正的基礎上，利用ENVI 軟件FLAASH 模塊計算相應波段反照率，再根據Liang[20]提出的針對陸地衛星的參數化方程得到，具體公式如下：

式中：b1、b3、b4、b5、b7為對應波段的地表反照率。

3 結果驗證與分析

3.1 太陽輻射估算結果驗證

3.1.1 正午時刻太陽輻射瞬時值驗證

由于研究區所處位置地形復雜，受天氣、地面站點數量以及遙感影像獲取時間限制，實地觀測數據較為困難，文中以紅原站2016年日輻射觀測數據（含正午時刻輻射強度）為基準對估算結果進行驗證。為使結果符合實地情況，同時又能準確驗證其誤差，本文選取了5月12日、7月15日、10月3日及12月22日共4個日期來反映2016年春、夏、秋、冬四季變化情況。綜合上述太陽直射、天空散射輻射及周圍地形反射輻射模型，利用公式（1）得出研究區2016年不同季節正午瞬時太陽輻射，以2016年5月12日為例（圖2）。

圖2 2016年5月12日研究區正午瞬時太陽輻射估算結果

文中取氣象站點紅原站所在位置像元周圍3×3窗口的地表瞬時太陽輻射均值作為站點的估算結果，將該結果與紅原站正午時刻輻射強度觀測值進行對比（表2）。

表2 研究區正午時刻瞬時太陽輻射估算結果驗證

3.1.2 太陽輻射日總量驗證

研究表明，瞬時太陽輻射由于各參數取得難度較大且計算方法不同存在較大偏差，而更多的是在時間尺度上以日或更長的時間為單位來獲取總量[21]。太陽輻射在一天當中呈現周期變化規律，某一時刻太陽輻射的瞬時值和一天中輻射最大值具有線性函數關系。本文采用Jackson 的方法將瞬時輻射值向上推算為日總量數據[22-23]。具體公式如下：

式中：Iz為 太陽輻射日總量；Im為估算的正午時刻紅原站瞬時輻射最大值；n為日照時長，由氣象站點提供；t為時間變量。根據上述結果計算得到研究區太陽輻射的日總量輻射值（表3）。

表3 研究區太陽輻射日總量估算結果驗證

從表2、表3 結果上看，估算精度較好。正午時刻瞬時太陽輻射誤差均值為13.56 W/m2，日總量誤差均值為1.49 MJ/m2，整體結果良好。由于本文估算正午時刻的瞬時太陽輻射，進而求得日總輻射量，因此隨著不同季節太陽高度角的變化，誤差結果呈現群分現象。這是因為部分參數是由溫度、濕度等氣象數據計算得到，而研究區內可使用氣象站點較少，故采用均值代替逐時觀測數據代入公式（3）計算的大氣透過率，未考慮大氣自身氣溶膠含量和水汽含量等差異帶來的影響。

但總體而言，該方法能較好地估算山地地表太陽輻射，精度相對可靠，其結果能反映山地地表太陽輻射的空間變化，可以滿足后期的分析需求。

3.2 太陽輻射時空特征分析

經過輻射模型計算，2016年5月12日、7月15日、10月3日和12月22日正午時刻太陽高度角分別為68°44′、70°51′、51°17′和30°10′。山區太陽輻射太陽最大輻射值分別是1239.89 W/m2、1188.08 W/m2、1325.31 W/m2、1372.88 W/m2，最小值是7.33 W/m2、39.79 W/m2、9.83 W/m2、4.34 W/m2， 平 均 值 是947.89 W/m2、 939.69 W/m2、792.52 W/m2、575.58 W/m2。

研究區輻射分布受地貌影響強烈，具有明顯的地形規律。從結果上看，輻射分布與該區地貌特征規律相符合，呈現山脊比山谷大且陽坡大于陰坡的分布特點，最大值出現在東北方向坡面上。

3.2.1 坡度與太陽輻射的關系

山地坡面朝向不同，使得坡地上各時刻接收的太陽輻射通量各不相同[12]。文中利用ArcGIS 將研究區坡度按一定間隔劃分等級，并提取各等級像元數，統計相應級別的均值，使用ArcGIS 空間分析功能得到坡度變化與研究區不同時期太陽輻射的關系（圖3）。如圖所示，總體上隨著坡度的增大而減小，不同季節太陽輻射呈遞減趨勢。由于所選日期較為接近且川西春、夏兩季氣候相似，變化趨勢大體相當；秋季隨著太陽直射位置的移動，太陽輻射強度的減弱變化趨勢漸緩；冬季太陽輻射強度的減弱變化趨勢最慢。

圖3 研究區不同季節坡度與太陽輻射的關系

通過統計分析發現，春、夏兩季在坡度超過20°以后輻射量下降趨勢明顯加快，而秋、冬兩季太陽輻射受坡向影響，在坡度45°～55°之間略微上升。坡度超過55°以后，輻射量呈大斜率減小。各季節太陽輻射變化趨勢說明山區受太陽直射輻射、坡度、坡向及周圍地形遮蔽作用影響越大，輻射量變化也就越明顯[12]。

3.2.2 坡向與太陽輻射的關系

與坡度分析相似，將研究區坡向劃分等級。為便于統計，坡向按16 方位制劃分（圖4）。如圖所示，各季節正午瞬時太陽輻射最大值均以東北方向為中心，主要集中分布于22.5°～67.5°之間的坡面，5月和7月最為接近，10月和12月隨季節變化輻射量逐漸減小。這與正午太陽所處位置有關，直接輻射的高低與太陽所在方位密切相關，太陽直接輻射大，坡面的太陽輻射也就最大。坡面太陽輻射最小值出現在西南坡向，這是由于隨著太陽入射角的變小，同時受山體陰影遮蔽，接收的太陽直接輻射相對減小，太陽輻射也就最小。另外，坡向對太陽輻射的影響隨季節變化的趨勢也較為顯著。冬季在北回歸線以北太陽高度角最低，太陽直射強度最低，陰坡同時受地形遮蔽影響，輻射量也最小。

圖4 研究區坡向與太陽輻射的關系

3.2.3 植被與太陽輻射的關系

利用Landsat8 影像在ENVI 軟件中計算研究區植被指數（Normalized Difference Vegetation Index，NDVI），經過歸一化處理，按均值-標準差將其劃分為5個等級（圖5）。由圖上看出5月、7月變化趨勢大體一致，10月、12月在植被指數介于0.2～0.4 時，太陽輻射有略微上升趨勢，超過0.4 后下降趨勢較為明顯。上述特征表明隨著植被指數增大，地表植被覆蓋就越高，地面接收的太陽輻射總體上呈現遞減趨勢。這是因為植被對地表的遮蔽同時能夠吸收和反射一部分太陽輻射使得到達地面的太陽輻射減少。

圖5 研究區植被指數與太陽輻射關系

4 結論與討論

本文利用已有山區太陽輻射模型，結合DEM、氣象數據和遙感影像，綜合考慮山區各個因素，估算了川西高原研究區2016年不同時期正午瞬時山地太陽直接輻射、天空散射輻射和周圍地形輻射，進而得到研究區正午瞬時太陽輻射，并對估算結果進行驗證分析，主要結論如下：

（1）經過氣象站點觀測數據的驗證，山區瞬時太陽輻射估算值的平均誤差為13.56 W/m2，日總量誤差均值為1.49 MJ/m2，模型精度較好，符合研究區實際情況。

（2）山地太陽輻射分布與地形地貌密切相關，總體上陽坡大于陰坡，山脊大于山谷。

（3）整體上，山區太陽輻射隨坡度增大而遞減。同時受坡向、植被覆蓋度等因素影響，坡面輻射量在一定坡度有略微上升趨勢。當坡度大于太陽天頂角時，輻射量下降明顯。

（4）山區受植被覆蓋影響，隨著地表植被的增加，地表輻射量呈下降趨勢。

本文在前人研究的基礎上，結合氣象數據和遙感對山區地表太陽輻射進行估算，結果精度較高，參數可靠。但在估算過程中，還存在一些不足之處。如部分所用參數采用氣象站點觀測值均值代替，不可避免的帶來估算誤差，有待改進；受數據源限制，僅采用一個站點觀測數據對結果進行驗證，精度需進一步提高；研究區受天氣影響較為突出，在計算中未考慮云霧、氣溶膠對大氣透過率的影響等。上述問題將在后續的研究工作中逐步解決。