1 引言

太阳辐射是植物进行光合作用所需外部能量的唯一来源, 植物在进行光合作用时对太阳辐射能的吸收和利用具有一定的选择性, 通常把波长在400~700 nm范围内能被植物光合作用利用的那部分太阳辐射称为光合有效辐射(Photosynthetically Active Radiation, PAR)(姚济敏等, 2005; 王文玲等, 2013)。PAR是一种重要的气候资源, 到达植被表面的PAR是植被生物量形成的基本能源(周允华等, 1996; 李英年和周华坤, 2002)。PAR是影响植物光合作用的几个主要环境因子(气温、相对湿度、光合有效辐射及CO₂浓度等)之一, 它与植物的光合速率有着显著的相关关系(Chen et al, 2012)。因此, PAR也是陆地生态系统碳循环研究和定量估算中的重要参数。植被净初级生产力(Net Primary Productivity, NPP)是表征陆地生态系统碳循环的重要参数, 也是陆地生态系统碳汇的主要决定因子(朱文泉等, 2007; Field et al, 1998)。NPP是生态系统中绿色植物通过光合作用所产生的有机物总量扣除自养呼吸后的剩余部分, 光合作用是植物与环境间物质与能量交换的核心过程, 大多数生态模型都以这一过程为出发点(刘荣高等, 2004)。作为NPP估算模型的一个重要输入参数, PAR估算的不确定性是目前遥感估算NPP精度不高的重要原因之一(孙睿和朱启疆, 2001)。

光合有效辐射可通过地面光谱仪直接测得。然而, 目前大部分地面观测站仅展开总辐射量的观测, 进行光合有效辐射观测的为数不多。正因为此, 目前空间连续的光合有效辐射多由观测数据插值或经验公式估算得到(林忠辉等, 2002; 张运林和秦伯强, 2002), 如季国良等(1993)根据实际观测资料拟合得到了张掖地区PAR的经验公式。随着遥感技术的不断发展, 遥感成为获得较高分辨率PAR的一种可行性方法(董泰锋等, 2011)。遥感数据获取周期短、成本低, 因而通过遥感反演估算已成为当前获得大范围光合有效辐射的有效途径。目前国内外已先后提出了多种基于辐射传输模型的光合有效辐射遥感反演算法, 但不同算法获得的数据差别较大。Eck and Dye(1991)建立了370 nm通道的紫外反射率与PAR的线性估算模型, Pinker and Laszlo(1992)发布了第一个以月为单位的全球范围内250 km的PAR数据产品ISCCP-PL。美国国家航空航天局(NASA)发射的Terra与Aqua两颗卫星为遥感反演高精度高分辨率的PAR数据提供了新的途径。其星上搭载的MODIS传感器包含经过精确定标的36个波段, 可同时获得影响陆表光合有效辐射的气溶胶、云、臭氧、大气水汽含量及地表反射率等参数信息, 并具有百米至千米级的空间分辨率和至少一天两次的较高时间分辨率。刘荣高等(2004)以查找表的方式, 以MOIDS数据为基础反演估算了晴空下PAR, 其估算结果与气溶胶的空间分布一致, 说明华北平原晴空下的PAR主要受到气溶胶的影响。Liang et al(2006)利用大气辐射传输模型模拟出不同大气状况下大气层顶辐亮度与PAR的关系。谢小萍等(2009)利用MODIS的大气和陆地产品考虑了瑞丽散射、臭氧和水汽吸收、气溶胶散射以及地表和大气间的多次反射, 进行山东省PAR的估算, 其估算结果与禹城站和长白山站实测值有较好的相关性。

复杂地形条件下, 海拔、坡度、坡向的分异以及周围地形遮蔽的影响使得山区地表接受太阳辐射的入射角度不同, 同时受到周围环境反射辐射的贡献也不同(李新等, 1999)。由于地形条件下的辐射计算公式过于复杂, 地形参数获取技术缺乏, 长期以来, 人们往往对其忽略或简化(何洪林等, 2003)。数字高程模型(DEM)描述地形高程的变化, 为精确计算地形对太阳辐射的影响提供了前提条件(田辉等, 2007)。李新等(1999)在计算任意地形条件下太阳辐射模型时采用了遮蔽因子、全天各向同性可见因子等地形参数, 改进了太阳辐射空间分布模拟过程中的几个重要参数的算法; 李净和李新(2007)引进李新模型中的地形参数, 并利用DEM获取了计算太阳辐射的地形参数, 建立了山区太阳辐射模拟的模型。

岷江上游地处四川盆地丘陵向川西北高原的过渡地带, 属青藏高原东缘高山峡谷区。该区海拔差异悬殊, 平均坡度达28. 8°, 地形非常复杂。本研究以岷江上游地区的复杂地形条件为切入点, 建立适合复杂地形条件下的光合有效辐射估算模型, 旨在提高复杂地形下光合有效辐射遥感估算精度, 可为进一步提高植被生产力估算模型的精度提供数据支持。

2 研究区与数据基础 2.1 研究区概况

岷江是长江上游水流量最大的支流。岷江上游是指岷江位于都江堰以上的流域(图 1)。该区位于102°35′E103°56′E, 30°45′N33°09′N之间。北部地区主要以高原地貌为主, 南部为典型的高山峡谷地形。区域内海拔落差较大, 海拔最高为6250 m(四姑娘山), 最低为764 m(东南漩口地区)。坡度为0~88°, 平均坡度达到28. 8°。岷江上游受地形因素和天气条件的影响, 全年日照时数在1500~1800 h之间, 日照百分率仅38%。一年内日照时数以4月、8月最高, 分别为145 h和155 h左右, 2月、9月最低, 日照时数仅为115 h和105 h左右。

2.2 数据基础 2.2.1 遥感数据

研究采用的遥感数据主要包括MODIS大气标准产品数据(https: //ladsweb.nascom.nasa.gov/)和ASTERG DEM数据(http: //www.gscloud.cn/)。MODIS大气标准产品数据来自MODIS产品数据集, 具体产品类型在表 1中列出。其中气溶胶光学厚度、大气臭氧含量、大气水汽含量用于估算大气中太阳辐射的消减。由于MODIS各产品空间分辨率不同, 且大气臭氧含量和气溶胶光学厚度在大气中是连续变化的, 研究中采用反距离权重插值方法将大气臭氧含量产品和气溶胶光学厚度产品重采样到与大气水汽含量产品相同的1 km分辨率, 并捕捉栅格使其与相应栅格对应。受到天气的影响, 遥感数据经常会出现一些异常值, 研究中修正异常值的方法为将所有正常值的栅格数据转换为点数据, 再通过ArcGIS空间分析功能中的克里金插值方法进行插值。考虑到非晴空条件下, MODIS大气水汽产品和MODIS气溶胶产品受云的影响, 异常值较多, 不利于本研究的进行, 因此选择了2013年22个晴空天气进行光合有效辐射的模拟与估算。

ASTER GDEM数据由日本METI和美国NASA联合研制。其基于“先进星载热发射和反辐射计(ASTER)”数据计算生成, 是目前唯一覆盖全球陆地表面的高分辨率高程影像数据。研究中采用30 m分辨率的ASTER GDEM数据产品, 用于提取每一个栅格的海拔、坡度和坡向。

2.2.2 辐射监测数据

文中采用地面实测光合有效辐射对估算结果进行验证。中国生态系统观测研究网络(http: //rs.cern.ac.cn/)提供了实时监测数据和每日总光合有效辐射, 研究选用了位于岷江上游的茂县山地生态系统定位研究站的逐日光合有效辐射总量和每日逐时光合有效辐射量数据。该站点位于四川省茂县静州村附近(图 1), 站址中心位置为103°53′E, 31°40′N, 海拔1820 m, 该站隶属于中国科学院成都生物研究所。站点采用芬兰Vaisala的Milos520系统进行辐射数据和常规气象要素的连续观测。光合有效辐射是用美国Lincoln生产的Li-190SA传感器来观测的, 数据采集频率为1 min一次。为了保证观测数据的可靠性, 采用美国Keithley公司2700型数字万用表每两年进行一次辐射传感器的标定, 控制光合有效辐射表的标定不确定度在±5%以内, 满足世界气象组织的精度要求(胡波等, 2012)。

3 岷江上游光合有效辐射估算模型的构建 3.1 瞬时光合有效辐射估算的基本原理

地球表面接收到的太阳辐射是太阳照射到地球大气层顶的辐射经过地球大气消光作用后的部分(黄鹤等, 2015)。晴空条件下, 大气对太阳辐射的影响主要为大气组分的吸收作用, 包括氧气和二氧化碳等混合气体、水汽、臭氧、气溶胶等, 以及大气分子的瑞丽散射和气溶胶散射作用。地表接受到的太阳辐射又会因为地表的坡度和坡向而反射一部分回到大气中。一般认为光合有效辐射(PAR)的波段范围为0. 4~0. 7 μm。计算公式为:

$ {\rm{PAR}} = \int_{0.4}^{0.7} {{I_\lambda }{\rm{d}}\lambda } , $

(1)

式中: I_λ为太阳下行单色光辐射值(单位: W·m^-2·μm^-1), λ为真空中的波长(单位: μm), 它可以表示为直接辐射(I_b)和散射辐射(I_d)的总和, 即:

$ {I_\lambda } = {I_b} + {I_d}, $

(2)

式中: I_b为直接辐射(单位: W·m^-2·μm^-1), 表示辐射穿过大气直接照射到地面的太阳辐射; I_d为散射辐射(单位: W·m^-2·μm^-1), 表示辐射在大气中经过多次散射后又到达地面的太阳辐射。

3.1.1 单波段大气层顶瞬时辐射的计算

单波段大气顶层辐射I_0λ可表示为:

$ {I_{0\lambda }} = {E_\lambda } \cdot \cos \left( \theta \right) \cdot {d^2}, $

(3)

式中: E_λ为标准日地距离下, 大气顶层单色光辐射强度, 数据来源于Bird and Riordan (1986); θ为太阳天顶角; d为以天文单位表示的日地距离, 其计算公式如下(Lillesand et al, 2000):

$ d = 1 + 0.01674\sin \left[ {\frac{{2{\rm{\pi }}\left( {J - 93.5} \right)}}{{365}}} \right], $

(4)

式中: J为儒略天数(Julian Day), 表示影像获取日期在太阳历中的排列序号, 如1月1日的排列序号为1, 12月31日的排列序号为365。

3.1.2 直接辐射的估算

直射单色光通量I_bλ的估算中主要考虑大气中的瑞利散射和大气水汽、大气中氧气、臭氧、气溶胶的吸收和散射等作用, 计算公式(Gregg and Carder, 1990)为:

$ {I_{b\lambda }} = {I_{0\lambda }}\cos \theta {T_{\rm{r}}}{T_{\rm{a}}}{T_{\rm{w}}}{T_{\rm{o}}}{T_{{\rm{oz}}}}, $

(5)

式中: I_0λ为单波段大气层顶瞬时辐射, θ为太阳天顶角, T_r、T_a、T_w、T_o、T_oz分别为臭氧、氧气等混合气体、水汽、气溶胶和瑞丽散射的透过率。大气中臭氧、氧气等混合气体、水汽、气溶胶和瑞丽散射的透过率计算公式分别如下:

$ {T_{\rm{r}}}\left( \lambda \right) = \exp \left[ { - M'\left( \theta \right)/\left( {115.6406{\lambda ^4} - 0.335{\lambda ^2}} \right)} \right], $

(6)

$ {T_{\rm{a}}}\left( \lambda \right) = \exp \left[ { - {\tau _{\rm{a}}}\left( \lambda \right)M\left( \theta \right)} \right], $

(7)

$ {T_{\rm{w}}}\left( \lambda \right) = \exp \left\{ { - \frac{{0.2385{a_{\rm{w}}}\left( \lambda \right){w_{\rm{v}}}M\left( \theta \right)}}{{{{\left[ {1 + 20.07{a_{\rm{w}}}\left( \lambda \right){w_{\rm{v}}}M\left( \theta \right)} \right]}^{0.45}}}}} \right\}, $

(8)

$ {T_{\rm{o}}}\left( \lambda \right) = \exp \left\{ { - \frac{{1.41{a_{\rm{o}}}\left( \lambda \right)M'\left( \theta \right)}}{{{{\left[ {1 + 118.3{a_{\rm{o}}}\left( \lambda \right)M'\left( \theta \right)} \right]}^{0.45}}}}} \right\}, $

(9)

$ {T_{{\rm{oz}}}}\left( \lambda \right) = \exp \left[ { - {a_{{\rm{oz}}}}\left( \lambda \right){H_{{\rm{oz}}}}{M_{{\rm{oz}}}}\left( \theta \right)} \right], $

(10)

式中: M(θ)为大气路径长度, M_oz(θ)为臭氧传播路径长度, M′(θ)为气压归一化后的大气路径长度; a_oz(λ)为臭氧的吸收系数, H_oz为臭氧垂直方向上的臭氧质量; a_o(λ)为氧气的吸收系数; a_w(λ)是水汽吸收系数, w_v为大气含水量(单位: cm); τ_a(λ)为气溶胶光学厚度; λ为波长(单位: μm)。其中, a_oz(λ)、a_o(λ)、a_w(λ)根据Neckel and Labs的研究成果取值如下(表 2):

3.1.3 散射辐射的估算

散射单色光通量I_dλ的估算主要考虑, 辐射在传播过程中的瑞丽散射分量I_rλ和气溶胶散射分量I_aλ, 计算公式为:

$ {I_{d\lambda }} = {I_{r\lambda }} + {I_{a\lambda }}, $

(11)

瑞丽散射和气溶胶散射的计算公式分别如下:

$ {I_{r\lambda }} = {I_{0\lambda }}\cos \theta {T_{{\rm{oz}}}}{T_{\rm{w}}}{T_{\rm{a}}}\left( {1 - {T_{\rm{r}}}} \right) \times 0.5, $

(12)

$ {I_{a\lambda }} = {I_{0\lambda }}\cos \theta {T_{{\rm{oz}}}}{T_{\rm{w}}}F{w_{o\lambda }}{T_{\rm{r}}}\left( {1 - {T_{\rm{a}}}} \right), $

(13)

式中: I_0λ为单波段大气层顶瞬时辐射, θ为太阳天顶角, T_oz、T_w、T_r、T_a分别为臭氧、水汽、瑞丽散射和气溶胶的透过率; F为前向散射, w_0λ为单色光单次散射反射率。其中, F和w_0λ分别按以下公式计算:

$ F = 0.9302 \times \cos {\theta ^{0.2556}}, $

(14)

$ {w_{0\lambda }} = 0.945 \times \exp \left[ { - 0.095 \times {{\left( {{\rm{I}}{{\rm{g}}_{10}}\left( \lambda \right)/0.4} \right)}^2}} \right], $

(15)

式中: θ为太阳天顶角, λ为波长(单位: μm)。

3.2 瞬时光合有效辐射的地形校正 3.2.1 坡面上光线入射角的校正

在地表起伏较大的山地地区, 地表为存在一定坡度与坡向的倾斜面, 因此, 地表接受到的太阳光线的入射角与地面的坡度与坡向有关。在这种情况下, 坡面上太阳入射角r可以由坡度、坡向、太阳高度角和方位角校正, 计算公式为:

$ \begin{array}{l} r = \arccos \left( {S \cdot N} \right)\\ \;\; = \arccos \left( {\left[ {\begin{array}{*{20}{c}} {\sin {Z_{\rm{s}}}\mathit{cos}{{\rm{A}}_s}}\\ {\sin {Z_{\rm{s}}}\mathit{cos}{{\rm{A}}_s}}\\ {\cos s} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {\sin s\cos A}\\ {\sin s\sin A}\\ {\cos s} \end{array}} \right]} \right), \end{array} $

(16)

式中: S为太阳光线向量; N为地表法线向量; Z_s和A_s为太阳天顶角和方位角; A和s分别为地表坡度和坡向。

3.2.2 地形遮蔽的校正

在获得太阳高度角和方位角后, 山区太阳辐射的计算还必须考虑研究区域是否处于周围复杂地形的阴影之下, 是计算光合有效辐射必须考虑的一个重要因素。岷江上游地区地形复杂, 多沟壑, 地形遮蔽对太阳辐射的影响非常显著。本文采用地形遮蔽系数(V)来表示, 地形遮蔽系数可根据太阳高度角和太阳方位角方向上各点的高程确定。以某一像元为起点, 在一定的太阳方位角方向上进行追踪, 当该像元与目标像元之间的高度角大于当时太阳高度角, 则认为该像元未被遮蔽, 该像元的遮蔽系数为0;否则, 认为该像元被遮蔽, 遮蔽系数为1。当受到地形遮蔽时, 该像元只能接受到散射辐射。遮蔽系数V对太阳辐射的影响可表示如下:

$ \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{I_\lambda } = {I_{b\lambda }} + {I_{d\lambda }},}\\ {{I_\lambda } = {I_{d\lambda }},} \end{array}}&{\begin{array}{*{20}{c}} {V = 0}\\ {V = 1} \end{array}} \end{array}} \right. $

(17)

地形遮蔽系数V可利用ArcGIS空间分析功能的山体阴影计算结果二值化得到:

$ \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {V = 1,}\\ {V = 0,} \end{array}}&{\begin{array}{*{20}{c}} {0 < {\rm{Hillshade}} < 1}\\ {1 \le {\rm{Hillshade}} \le 255} \end{array}} \end{array}} \right. $

(18)

式中: Hillshade为山体阴影。在ArcGIS中, 山体阴影Hillshade通过为每个栅格单元指定太阳高角度等信息, 来计算该单元陆地表面的假定亮度值而得到, 其计算公式如下(Burrough and Mcdonnell, 1998):

$ \begin{array}{l} {\rm{Hillshade}} = 255.0 \times \left( {\left( {\cos \left( {{\rm{Zenith}}} \right) \cdot \cos \left( {{\rm{slope}}} \right)} \right)} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + \left( {\sin \left( {{\rm{Zenith}}} \right) \cdot \sin \left( {{\rm{slope}}} \right)} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. { \cdot \cos \left( {{\rm{Azimuth}} - {\rm{Aspect}}} \right)} \right)} \right), \end{array} $

(19)

式中: Hillshade为山体阴影, Zenith表示太阳天顶角, slope表示地表坡度, Azimuth表示太阳方位角, Aspect表示地表坡向。

3.2.3 邻近地形的反射辐射增强的校正

各种下垫面性质不同, 其吸收和反射太阳辐射的能力也是有很大差异的(张一平等, 2004)。复杂地形条件下, 周围地形反射过来的光合有效辐射也是重要的一部分(王晶晶, 2006)。下垫面对太阳辐射的性质可以通过地表反射率反映, 地表反射率是指地表对入射的太阳辐射的反射通量与入射的太阳辐射通量的比值, 表示下垫面对辐射能的吸收能力。本文中地表反射率采用刘三超等(2003)在黑河流域的不同地表覆被类型的实测地表反射率(表 3)。周围地形的反射辐射采用Hay(1980年)的计算方法(安顺等, 1988), 计算公式为:

$ {E_r} = {I_\lambda }\rho {\sin ^2}\left( {\frac{{slope}}{2}} \right) = 0.5\left( {{I_{d\lambda }} + {I_{b\lambda }}} \right)\rho \left( {1 - \cos \left( {slope} \right)} \right), $

(20)

式中: E_r为周围地形的反射辐射, ρ为坡面反射率, slope为地表坡度。

3.3 PAR日数据的时间尺度上推方法

在一些相关研究中, 瞬时太阳辐射量并不能直接作为输入参数, 而需要PAR的时间尺度一般是日单位或是更长的时间单位(Wang et al, 2010), 如:植被生产力的估算等往往需要PAR的日总量信息。本文采用Jackson的方法将瞬时PAR值向上推算为PAR日总量数据(Jackson et al, 1983), 具体方法如下:

(1) 由于太阳辐射在一天当中具有周期变化的规律, 因此, 某一时刻太阳辐射的瞬时值和一天中的最大值具有如下函数关系:

$ {S_i} = {S_{\rm{m}}}\sin \left( {\frac{{\pi t}}{N}} \right), $

(21)

式中: S_i为某一时刻的太阳辐射瞬时值; S_m为正午时太阳辐射的最大值; N为太阳升起到落下的时间间隔, 即日照时间; t为时间变量, 取值为0到N。当t=N/2时, 是S_i= S_m。

(2) 通过众多的观测数据发现, 在某一特定地区, PAR与太阳总辐射能的比率基本上是长期不变的(Weiss A et al, 1985)。因此, 一天中PAR随时间t的变化函数应该与太阳总辐射随时间t的变化函数一样, 即:

$ PA{R_i} = PA{R_{\rm{m}}}\sin \left( {\frac{{\pi {t_p}}}{N}} \right), $

(22)

式中: PAR_i为某一时刻的光合有效辐射瞬时值; PAR_m为正午时光合有效辐射的最大值; N为太阳升起到落下的时间间隔, t_p为时间变量, 取值为0到N。

那么, 对上式进行时间积分可得到PAR的日总量PAR_d的计算公式:

$ {\rm{PA}}{{\rm{R}}_d} = \int_0^N {{\rm{PA}}{{\rm{R}}_{\rm{m}}}\sin \left( {\frac{{\pi {t_p}}}{N}} \right){\rm{d}}{t_p} = \frac{{2N}}{\pi }{\rm{PA}}{{\rm{R}}_{\rm{m}}}} , $

(23)

其中, 每日的日照时间即太阳升起到落下的时间间隔N的计算方法如下:

$ N = a + b{\sin ^2}\left[ {\frac{{\pi \left( {J + 10} \right)}}{{365}}} \right], $

(24)

式中: J是一年中的儒略日期; a、b为常数, 取决于地理纬度, 可利用纬度(L)由下列公式计算得到:

$ \begin{array}{l} a = 12 - 5.69 \times {10^{ - 2}}L - 2.02 \times {10^{ - 4}}{L^2} + 8.25\\ \;\;\;\;\;\; \times {10^{ - 6}}{L^3} - 3.15 \times {10^{ - 7}}{L^4}, \end{array} $

(25)

$ \begin{array}{l} b = 0.123L - 3.10 \times {10^{ - 4}}{L^2} + 8.00 \times {10^{ - 7}}{L^3}\\ \;\;\;\;\;4.99 \times {10^{ - 7}}{L^4}. \end{array} $

(26)

4 结果与分析

研究基于处理后的MODIS标准大气产品, 采用上述模型分别估算了岷江上游2013年22个晴空条件下MODIS卫星白天过境时刻的瞬时光合有效辐射以及当日光合有效辐射总量, 并利用中国生

态系统研究网络茂县站的同时刻光合有效辐射和同期日光合有效辐射总量的监测数据对估算结果进行了验证分析。为了降低由噪声引起的反演结果的不确定性, 采用气象站点所在像元为中心的3×3像元窗口的光合有效辐射平均值作为该站点的反演结果。

4.1 反演结果的精度评价

通过模型估算的瞬时光合有效辐射与地面实时实测数据对比(图 2a)及经时间拓展的当日光合有效辐射总量与实测日总光合有效辐射的对比(图 2b)可以发现, 模型估算结果与实测数据的相关系数R分别达到了0. 8738和0. 8298, 表明模型反演结果与地面实测值有可信的解译度。因此该模型能够较好的模拟估算岷江上游复杂山区地形的光合有效辐射。

由模型估算结果和地面实测值基本统计信息(表 4)可以看出, 由本模型模拟估算瞬时光合有效辐射与地面实测值的最大误差为15. 1%;本模型推算的日光合有效辐射总量与地面实测日总量的最大误差为18. 7%;其他各统计特征值也很接近, 表明模型反演精度较高。

4.2 地形差异对地表接受到的光合有效辐射影响

为了检验岷江上游地区地形对光合有效辐射估算的影响, 在考虑和不考虑地形因素影响的情况下分别估算了岷江上游地区的光合有效辐射(图 3)。图 3a可以看出模型反演结果与地面实测值变化趋势基本一致, 表明模型有较高的准确度。在进行了地形因素校正之后, 模型的估算精度有明显的提高。从验证结果发现, 考虑了地形因子后, 模型反演误差较大的8~10天有明显提高。从总体看, 不考虑地形校正的模型估算结果明显高于地面实测值, 考虑了地形校正之后, 模型的估算结果有明显降低。这是因为模型中考虑地形的坡度、坡向和遮蔽作用减少了到达地表的光合有效辐射。从图 3b可以看出, 由瞬时光合有效辐射时间上推到光合有效辐射日总量后, 模型估算结果与地面实测值有较好的相关性。

由表 5可以看出, 当不考虑地形校正时, 瞬时光合有效辐射估算结果与地面观测值的平均绝对误差为37. 6 W·m^-2·s^-1, 平均相对误差为9. 14%;光合有效辐射日总量估算结果与地面观测值的平均绝对误差为0. 73 MJ·m^-2, 平均相对误差为7. 21%;考虑地形校正情况下, 瞬时光合有效辐射的平均误差都有减小, 平均绝对误差减小到26. 3 W·m^-2·s^-1, 平均相对误差减小到7. 26%;而光合有效辐射日总量的平均误差变化并不明显, 还有小幅的升高, 平均绝对误差为0. 76 MJ·m^-2, 平均相对误差为7. 26%。由此可见, 本模型能够较好的模拟估算岷江上游地区的光合有效辐射分布状况, 并在考虑地形校正后, 对模型估算瞬时光合有效辐射精度有一定的提高, 而对光合有效辐射日总量并没有提高, 这是因为对总光合有效辐射日总量, 地形的影响是随着时间变化的, 而本文中时间尺度上推方法中并没有考虑地形对光合有效辐射影响的日变化情况, 另外, 在瞬时光合有效辐射上推到光合有效辐射日总量过程中本身存在一定的误差。

5 结论与讨论

光合有效辐射是植被生态过程模型中的重要能量来源, 是大部分植被生产力估算模型的主要输入参数。本研究通过简化的辐射传输模型, 利用MODIS标准大气产品并结合地形校正建立了适合于复杂地形区域的瞬时光合有效辐射估算模型。通过估算2013年22个晴空条件下的岷江上游地区瞬时光合有效辐射, 并结合茂县生态研究站的地面监测数据进行验证分析得出如下结论:

(1) 本研究提出的基于MODIS大气产品和地形校正的晴空下光合有效辐射的反演估算方法具有较高精度。地形校正后, 模型估算的岷江上游地区瞬时光合有效辐射与地面实测值两者的相关系数达到了0. 8738, 平均误差为7. 26%。因此该模型估算结果是可信的。

(2) 为了满足生理生态模型对以日或更长时间尺度的光合有效辐射总量的数据需求, 通过将估算的瞬时光合有效辐射通过时间尺度上拓展得到日总光合有效辐射, 与实测数据相关系数达到了0. 8298, 估算的平均误差为7. 26%, 说明使用本方法将瞬时光合有效辐射上推到日总光合有效辐射适用于岷江上游地区, 适用于生理生态模型中。

(3) 在地形复杂区域, 地形校正能够提高有效光合有效辐射的估算精度。地形校正后, 瞬时光合有效辐射平均相对误差由9. 14%减小到7. 26%。

海拔、坡度、坡向的分异以及周围地形遮蔽的影响使得来自周围环境反射辐射的贡献不同。地形影响是光合有效辐射地表分布差异的主要原因。本模型估算的瞬时光合有效辐射结果与刘荣高等(2004)和李净和李新(2007)的研究结果一致, 说明该模型能够较准确的估算岷江上游地区瞬时光合有效辐射。但是, 研究中为了准确的进行估算地形起伏不明显地区的地形影响, 选用了30 m空间分辨率的DEM数据进行地形校正, 因此在进行MODIS数据和DEM数据叠加运算上, 存在一定的尺度转换误差。

需要指出的是, 研究中提出的复杂地形条件下的光合有效辐射遥感估算模型只适合于晴空大气模式, 模型中考虑了气溶胶、混合气体、水汽等大气影响因子, 没有考虑云的影响, 而云作为辐射传输中的重要参数且不确定性较大(王丽娟等, 2015), 云的空间分布差异复杂, 且日-云-地之间的几何关系十分难以模拟, 这在本研究中没有涉及, 也是随后工作重点需要考虑的问题。

研究模型估算的PAR为卫星过境时刻的瞬时PAR, 而地面站点监测数据为每小时的PAR平均值, 这使得利用站点监测数据对模型估算结果进行验证时存在一定的误差。另外, 研究中将瞬时PAR利用经验公式拟合得到日最大值和PAR日总量, 但实际中PAR日总量是随着一天内天气变化而变化的。多源遥感数据的综合应用, 特别是静止卫星与极轨卫星的结合, 应该有助于更好地拟合得到PAR日总量。但由于受研究区所处特殊地理环境及其气候特点的影响, 满足质量要求的卫星遥感数据的收集存在一定的困难, 而其解决方法受到限制。

研究中采用地面实测站点数据验证模型反演精度时, 是采用气象站点所在像元为中心的3×3像元窗口反演结果的平均值与站点实测数据进行对比验证, 在站点的最优代表空间和影像像元之间存在一定的尺度效应, 影响验证精度; 又在研究区域内只有一个可以作为验证的地面监测站点, 难免说服力不够, 随后工作中将扩大研究区域增加地面验证站点数据, 验证模型的适用性。

致谢 国家生态系统研究网络提供的实测数据为本研究的顺利进展提供了数据支持。