基于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）。空间分辨率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）山区受植被覆盖影响，随着地表植被的增加，地表辐射量呈下降趋势。

本文在前人研究的基础上，结合气象数据和遥感对山区地表太阳辐射进行估算，结果精度较高，参数可靠。但在估算过程中，还存在一些不足之处。如部分所用参数采用气象站点观测值均值代替，不可避免的带来估算误差，有待改进；受数据源限制，仅采用一个站点观测数据对结果进行验证，精度需进一步提高；研究区受天气影响较为突出，在计算中未考虑云雾、气溶胶对大气透过率的影响等。上述问题将在后续的研究工作中逐步解决。