一种基于遥感的林下光照强度估测方法

摘要

本发明公开了一种基于遥感的林下光照强度估测方法，包括以下步骤：步骤1：确定待测区域与光照相关的参数；所述的光照相关的参数包括一天当中可能受太阳辐射的总时间h

说明书

技术领域

本发明涉及一种基于遥感的林下光照强度估测方法。

背景技术

林下光照强度对林下环境异质性及动植物的生理生态过程有着重要作用，受到的关注日益增加，但目前仍处于起步阶段。由于受到太阳高度角和方位角的日变化和季节性变化、复杂的大气状况和局部地面环境的影响，林下光照不仅存在日变化和季节变化，同时也表现出复杂的三维空间异质性，所以测量林下光照强度十分困难。传统的直接测量方法(如照度计)只能监测到林下有限数量测量点的光照强度在监测时间内的变化。间接测量方法(如半球面影像方法)虽具有简单、快速等优点，但一次也只能测量某一点的光照强度。上述两类方法在不同程度上存在工作量大、过程冗繁、破坏性强、耗时耗力等问题。遥感技术凭借长期、连续、稳定、高效的观测优势，在生态环境监测研究中得到了广泛应用，但至今，并没有一种通用、可重复测量森林区域光照强度的方法；同时，对林下光照的测量研究极少。因此，找到一种快速、客观、可重复、以及可信度较高的林下光照监测方法具有重要意义，对认识开发我国森林事业具有重要的作用。

因此，有必要设计一种基于遥感的林下光照强度估测方法。

发明内容

本发明所要解决的技术问题是提供一种基于遥感的林下光照强度估测方法，该基于遥感的林下光照强度估测方法易于实施。

发明的技术解决方案如下：

一种基于遥感的林下光照强度估测方法，其特征在于，包括以下步骤：

步骤1：确定待测区域与光照相关的参数；

所述的光照相关的参数包括一天当中可能受太阳辐射的总时间h_t(又称为太阳辐射时间h_t)以及太阳位置参数；

太阳位置参数包括太阳高度角α和方位角α_s；

α用于水平表面太阳直射辐射计算，参见公式14，以及用于局部入射角方向散射辐射(W/m²)计算，参见公式19等等；

α_s用于山体阴影影响系数计算，参见公式18；

h_t用于计算太阳时角h_s，参见公式6。

步骤2：计算冠层上方太阳总辐射量；

待测区域(即具体地形环境)中冠层上方太阳总辐射量由局部入射角方向的太阳直射辐射、散射辐射以及地表反射辐射组成，即：

G_i＝B_i+D_i+R_i；

其中G_i、B_i、D_i和R_i分别为冠层上方太阳总辐射量、局部入射角方向的太阳直射辐射量、散射辐射量以及地表反射辐射量；

利用卫星遥感影像估测大气气溶胶、水汽、云层对太阳辐射的吸收与散射特性，根据地表太阳辐射估算模型逐时计算大气状况下具体地形环境中冠层上方局部入射角方向的太阳直射辐射量B_i和散射辐射量D_i；

利用地表反照率计算地表反射辐射量R_i；

步骤3：利用卫星遥感影像反演叶面积指数，基于比耳—朗伯定律，根据冠层上方太阳总辐射计算得出林下光照强度。

可以进一步生成区域尺度林下光照强度时空分布图。

步骤2中，局部入射角方向太阳直射辐射量的计算公式如下：

B_i＝I_n·δ_i·hillshade；式中δ_i为光线局部入射角，hillshade为山体阴影影响系数；具体地，有：

光线局部入射角：

山体阴影影响系数：

hillshade＝(cosα·cos Slop)+[sinα·sin Slop·cos(α_s-Asp)]；

式中Slop为地形的坡度，Asp为地形的坡向(坡度是指坡面的倾斜程度，坡向是指地形坡面的朝向)。

δ_s为太阳赤纬角，有δ_s＝23.45sin(360°·(284+N)/365)，式中N为积日(一年中的第几天)；

L为地理纬度；

h_s为太阳时角。

步骤2中，局部入射角方向散射辐射量D_i的计算方法如下：

(1)在没有山体遮挡情况：

式中

A_LN＝α_s-Asp>s-ASP<＝π)

A_LN＝α_s-Asp-2π>s-ASP＞π)；

A_LN＝α_s-Asp+2π>s-ASP＜-π)

式中，(a)为太阳直射透过率，其中，I_o为大气外层辐射量：

I_o＝S₀·(1+0.0344cos(360°·N/365))；

式中S_o为太阳辐射强度常数1367(w/m²)；

I_n为地表法向太阳直射辐射：I_n＝0.9751·I_o·τ_R·τ_Gas·τ_Wv·τ_Ae·τ_Cloud；

式中τ_R为瑞利散射系数；τ_Gas为综合气体影响系数；τ_Wv为水汽影响系数；τ_Ae为气溶胶影响系数；τ_Cloud为云影响系数。

(b)D_h为水平表面散射辐射量(W/m²)：

式中T_LK为Linke大气浑浊指数(Linke>

E_sc为地球与太阳之间的平均距离和瞬时距离的比值的平方；

(2)在有山体遮挡(山体阴影)情况(δ_i＜0andα＞＝0°)：

D_i＝B_hF(Slop)

式中：

B_h为水平表面太阳直射辐射计算：

B_h＝I_n·sinα；

F(Slop)为顾及山体坡度的散射系数，

F(Slop)＝view(Slop)+0.25227·(sin Slop-Slop·cos Slop-πsin²(Slop/2))

view(Slop)为可视域范围内天穹比例：

view(Slop)＝(1+cos Slop)/2。

步骤2中，地表反射辐射量R_i的计算方法如下：

式中albedo为地表反照率。

步骤3中，林下光照强度采用下式计算：

I_under＝e^-k·LAIG_i；

其中，(1)k为消光系数；

消光系数k由地面调查点获取，根据实地样点测量得到的冠层下方和冠层上方的光照强度比值(I_under/G_i)与叶面积指数呈反向指数关系，利用非线性拟合关系计算得到。详见后文中公式25、26以及图1、2。

(2)LAI为叶面积指数，有LAI＝0.22RSR+0.59，其中，

ρ_B3、ρ_B4和ρ_B5分别表示Landsat陆地遥感卫星影像的红波段、近红外波段和中红外波段的反射率；

max()和min()分别表示最大和最小函数。上述参数从Landsat陆地遥感卫星影像中获取；叶面积指数中涉及的RSR，以及各个波段反射率，最大最小值都可利用相关软件直接从Landsat陆地遥感卫星影像中读取。

有益效果：

本发明的基于遥感的林下光照强度估测方法，填补了现有技术中林下光照强度估测方法的技术空白，是一种低成本、大尺度且可信度较高的测量方法。

针对林下光照强度的测定，现有技术多集中在实地测量技术范畴内，如照度计或基于鱼眼镜头的半球面影像方法都只能进行单个样方的测量，大范围的样方的数据获取费时费力，且很难短时间内重复操作，从而导致结果包含空间异质性信息。本发明利用免费的卫星遥感数据，结合数字高程模型和大气传输模型，估测复杂地形下森林冠层下方的光照强度，突破实地测量技术的瓶颈，能实现区域范围内林下光照信息估测，具有简单、低成本、可信度较高等特点。

本发明技术方案主要分两部分，即复杂地形下的太阳辐射估测，森林冠层叶面积指数估测，单从两部分内容而言，都有相关文档或技术参考，但将这两部分综合在一起估测林下光照的研究鲜有报道，本发明关键就是打通这个连接点，实现既定目标，同时对上述两部分技术都有所改进，提高估测精度。

(1)这个连接点由公式(26)实现，根据比耳—朗伯定律，建立林下光照强度与叶面积指数LAI、冠层上方太阳辐射之间的非线性关系，消光系数k由地面调查点获取。该方法将卫星遥感数据和地面观测数据相结合，以期获得更客观更高精度的估测结果。

(2)复杂地形下太阳辐射估测的改进：首先，在水平表面散射强度的计算上，引入了由国际太阳能权威专家Jaro Hofierka提出的大气浑浊指数(Linke turbidity factor，公式(15))，该指数综合考虑了大气中水汽和细小颗粒对太阳光散射的影响，结合地形因素，能更准确估测复杂地形下太阳散射强度。其次，在计算法向太阳直射强度时，综合考虑了瑞利散射、稳定气体、水汽、气溶胶以及云量等方面的影响，结果更加符合实际天气情况。最后，在计算地表反射强度时，利用Landsat卫星遥感数据的多波段组合方法计算地表反照率，求算精确到每个像元。

(3)森林冠层叶面积指数估测的改进：本发明利用Landsat卫星遥感影像中红波段B3、近红外波段B4和中红外波段B5数据组合反演LAI，比传统的基于归一化植被指数NDVI的反演方法更加精确，由于林分范围内植被类型基本一致，本发明方法的稳定性更好，即使在林冠边缘LAI变化复杂的区域，其反演效果也比传统方法的误差小。现阶段获取大面积森林叶面积指数最有发展潜力的技术之一是激光雷达，但因费用较高，推广有一定的难度，本发明的遥感数据是免费的Landsat卫星影像，虽然精度不及激光雷达，但具有成本低、可重复操作等优点，易于推广实施。

附图说明

图1为研究地1的LAI与(I_under/G_i)数学拟合关系示意图；

图2为研究地2的LAI与(I_under/G_i)数学拟合关系示意图；

图3为2016年5月4日研究地1林下光照总辐射强度(w/m²)分布图；

图4为2016年11月1日研究地2林下光照总辐射强度(w/m²)分布图；

图5研究样地林下光照总辐射强度野外调查与遥感估测线性关系；

图6为复杂地形下太阳辐射示意图；

图7为太阳高度角、方位角示意图；

图8为本发明流程图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：如图1～8，本发明利用卫星遥感影像数据估测复杂山地环境中森林冠层下方太阳辐射强度，同时顾及大气和地形效应，提高估测精度，首先，根据天体自传理论确定实验区一天当中可能受太阳辐射的总时间以及太阳位置等参数；然后，利用卫星遥感影像估测大气气溶胶、水汽、云层，以及其他气体成分对太阳辐射的吸收与散射特性，根据地表太阳辐射估算模型逐时计算“真实”大气状况下水平地表接受到的太阳直射和散射辐射；再次，结合数字高程模型DEM计算复杂地形下的光线入射角方向的太阳辐射参数；最后，利用卫星遥感影像反演叶面积指数，基于比耳—朗伯定律，根据冠层上方太阳总辐射计算得出林下光照强度，可以进一步生成区域尺度林下光照强度时空分布图。

太阳辐射时间及位置参数提取

太阳辐射时间h_t由日出太阳角度h_sr、地理纬度L、太阳赤纬角δ_s确定：

h_sr＝arccos(-tan>s)>

δ_s＝23.45sin(360°·(284+N)/365)>

(3)

式中N为积日(一年中的第几天)，round()为取整函数。

太阳高度角α和太阳方位角α_s由太阳时角h_s、地理纬度L和太阳赤纬角δ_s共同确定：

α＝arcsin(sin L·sinδ_s+cos>s>s)>

h_s＝h_sr-(π·t/h_t)>

式中t为时间间隔，例如0.5小时或1小时。

基于遥感的地表太阳辐射估算

主要包括地表法向太阳直射辐射和水平表面散射2方面计算。

(I)地表法向太阳直射辐射(W/m²)计算：

大气外层辐射(W/m²)：

I_o＝S₀·(1+0.0344cos(360°·N/365))>

式中S_o为太阳辐射强度常数1367(w/m²)

地表法向太阳直射辐射：

I_n＝0.9751·I_o·τ_R·τ_Gas·τ_Wv·τ_Ae·τ_Cloud>

式中τ_R为瑞利散射系数；τ_Gas为综合气体影响系数；τ_Wv为水汽影响系数；τ_Ae

为气溶胶影响系数；τ_Cloud为云影响系数。

(1)瑞利散射系数：

式中M为大气质量，M_p为经过高程校正的大气质量，

p＝exp(-0.0001184z)·p_o

Z为地面高程，p_o为标准气压。

(2)气体影响系数：

(3)水汽影响系数：

τ_Wv＝1-{2.4959γ[(1.0+79.034γ)^0.6828+6.385γ]^-1}>

式中γ＝w·M，w为大气可降水量(cm)。

(4)气溶胶影响系数：

式中κ_aλ为波长为380nm和500nm时的气溶胶光学厚度。但是，由于采用MODIS或FY-3A卫星提供的气溶胶光学厚度数据波长分别为470nm、550nm和650nm，因此需要进行一定的转换：

式中AOT₅₅₀和AOT₃₈₀分别为波长为550nm和380nm的气溶胶光学厚度；a为参数，Iqbal等研究者给出的建议值为1.3。

(5)云影响系数：

τ_Cloud＝1-CI>

式中CI为云总量数据，无量纲。

水平表面太阳直射辐射计算：

B_h＝I_n·sinα>

法向太阳直射辐射计算中大气可降水量w、0.55nm气溶胶光学厚度参数AOT₅₅₀以及云量CI等由MODIS系列卫星或FY-3A系列卫星影像估算得出。(II)水平表面散射辐射(W/m²)计算：

式中T_LK为Linke大气浑浊指数(Linke>

E_sc为地球与太阳之间的平均距离和瞬时距离的比

值的平方。Linke指数范围1～7，晴朗无云的天气，Linke为1～2，雨雪或雾

霾等天气，Linke取值较大。

复杂地形下太阳辐射估算

复杂地形环境中森林冠层上方太阳总辐射量由局部入射角方向的太阳直射辐射、散射辐射以及地表反射辐射3方面组成。

(A)局部入射角方向太阳直射辐射[W/m²]计算：

B_i＝I_n·δ_i·hillshade>

式中δ_i为光线局部入射角，hillshade为山体阴影影响系数。

光线局部入射角：

式中Slop为地形的坡度，Asp为地形的坡向。

山体阴影影响系数：

hillshade＝(cosα·cos Slop)+[sinα·sin Slop·cos(α_s-Asp)]>

3.3.2局部入射角方向散射辐射(W/m²)计算：

(1)在没有山体遮挡情况：

式中

A_LN＝α_s-ASP>s-ASP<＝π)

A_LN＝α_s-ASP-2π>s-ASP＞π)

A_LN＝α_s-ASP+2π>s-ASP＜-π)

(2)在有山体遮挡(山体阴影)情况(δ_i＜0andα＞＝0°)：

D_i＝B_hF(Slop)>

式中F(Slop)为顾及山体坡度的散射系数，view(Slop)为可视域范围内天穹比例：

F(Slop)＝view(Slop)+0.25227·(sin Slop-Slop·cos Slop-πsin²(Slop/2))

view(Slop)＝(1+cos Slop)/2

(B)复杂地形下地表反射辐射(W/m²)计算：

假设入射角表面呈漫反射状态：

式中albedo为地表反照率。

复杂地形环境中冠层上方太阳总辐射量[W/m2]由局部入射角方向的太阳直射辐射、散射辐射以及地表反射辐射组成，即：

G_i＝B_i+D_i+R_i>

林下光照强度估算

叶面积指数(LAI)是冠层结构中重要的参数之一，LAI是无量纲参数，定义为单位地表面积上所有叶片表面积的一半。已有研究者在实地应用光学仪器(如LAI—2000)测量LAI，虽然得到的精度较高，但仍然没有突破大范围、快速测量的瓶颈；也有研究者利用卫星遥感影像多波段数据反演区域尺度的LAI。本发明利用Landsat卫星遥感影像中红波段B3、近红外波段B4和中红外波段B5数据组合反演LAI，并取得了较好结果(R²＝0.75)：

LAI＝0.22RSR+0.59 (24)

根据比耳—朗伯定律，叶面积指数LAI和冠层下方与冠层上方的光照强度比值(I_under/G_i)呈反向指数关系，即：

式中G_i为公式(22)得出的冠层上方太阳总辐射量，k为消光系数，取决于叶倾角和光束方向，实验中利用公式(24)得到的LAI与实地样点测量的冠层下方与冠层上方的光照强度比值之间的数学关系，林下总光照强度即：

I_under＝e^-k·LAIG_i>

大气参数与地表反照率遥感估算

利用MODIS卫星遥感数据估算公式(11)中的大气可降水量w：

根据MODIS 17、18和19通道为大气吸收通道，2和5通道为大气窗口通道，按照比耳—朗伯定律，推导出大气水汽透过率为通道19和通道2的反射率比值τ_w＝ρ_B19/ρ_B2，对于复合性地表，大气可降水量w即为：

w＝(0.02-lnτ_w/0.651)²>

从而计算出MODIS遥感影像中每个像元的大气可降水量。

利用MODIS卫星遥感数据估算0.55nm气溶胶光学厚度AOT₅₅₀：

利用MODIS1、3和7通道数据按照经典的暗像元法反演地表真实反射率，再结合表观反射率、MODIS自带的经纬度坐标数据以及气溶胶模式等参数，带入6S辐射模型得到理论表观反射率值，当理论表观反射率与MODIS数据文件中的表观反射率相对应时即得到550nm的光学厚度AOT₅₅₀，用以计算公式(12)的气溶胶影响系数。

利用卫星遥感数据检测公式(13)中的云量CI：

由于云与雪在可见光波段具有相似的光谱特性，从而影响云监测的精度，因此，可以先把雪分离出来，然后进行云检测以提高检测精度。利用MODIS多波段组合法进行云检测，首先，利用MODIS2、5、7、和26通道数据检测冰雪，即：(ρ_B2-ρ_B5)/(ρ_B2+ρ_B5)＜0.15ANDρ_B7＜0.05ANDρ_B26＜0.05>

式中ρ为MODIS影像中波段B_i反射率，

然后再利用MODIS 3、8和26通道数据检测云：

ρ_B3＞0.2OR[ρ_B26＜0.02ANDρ_B8＞0.17]>

云量CI即为研究区范围内有云面积与总面积之比。

利用Landsat系列卫星影像估算公式(21)中的地表反照率：

albedo＝0.356ρ_B1+0.130ρ_B3+0.373ρ_B4+0.085ρ_B5+0.072ρ_B7-0.0018>

式中ρ为Landsat影像中波段B_i反射率。

针对2个具体的研究区实施本发明：

研究区概况：

研究地1位于湖南省东北部的福寿山林场(28°3′00″—28°32′30″N，113°41′15″—113°45′00″E)，地处罗霄山脉连云山支脉，地势南高北低，平均海拔1200余米，平均坡度为22—27度，呈现群山重叠的地貌。年平均气温12.1℃极端最高气温33.4℃最低气温为—15℃，年日照1500小时，无霜期217天。主要植被类型是典型的中亚热带常绿阔叶林，上层乔木树种主要有杉木、黄山松、青冈栎、苦槠、檫木、江南桤木、山核桃及壳斗科植物。

研究地2位于辽宁省抚顺市岗山国家森林公园(41°39′00″—41°32′30″N，125°23′00″—125°31′00″E)，山峦起伏、沟壑纵横，地势北高南低，主峰海拔1373米，森林植物群落垂直分布特征明显，样地主要植被类型是典型的温带针叶林，上层乔木树种主要有人工落叶松、红松以及沙松。

野外调查方法

在固定样地内设置60个30m×30m的样方，用差分GPS确定每个样方的位置，在每个样方中心和对角线四分位处用数码相机外接鱼眼镜头(广角为183°，正投影)获取半球面林冠影像，影像方向与磁北方向重合。用Gap Light Analyzer(GLA，V2.0，图像处理软件)对冠层照片进行分析并计算得出林下光照强度，计量单位是w/m²。研究地1的调查时间为2016年5月4日，天气情况为多云或有薄雾；研究地2的调查时间为2016年11月1日，天气情况较为晴朗。60个样方数据中40个作为训练样本数据，其余20个作为验证样本数据。

遥感数据获取

MODIS卫星遥感数据从美国国家航空航天局(NASA)官方网站https://ladsweb.modaps.eosdis.nasa.gov/免费下载，下载与研究地范围相符的MODIS一级数据产品，包括Terra和Aqua系列的MOD03经纬度坐标数据，MOD02定标辐射数据，影像生成时间与野外调查时间一致。

数字高程模型数据(DEM)获取分两种途径，一是从地形图(方里网)上矢量化生成(经费来源湖南省自然科学基金项目，编号:2015JJ2201)，二是从中国科学院地理空间数据云网站http://www.gscloud.cn/免费下载，下载与研究地相符的30m分辨率栅格形式DEM。

Landsat卫星遥感数据从中国科学院地理空间数据云网站http://www.gscloud.cn/免费下载，下载与研究地范围相符的多波段遥感影像，尽量选择与野外调查时间相近、少云或晴朗天气下生成的影像。

基于本发明方法对2个研究区域的估测过程：

水平地表太阳辐射估算

借助专业遥感图像处理软件The Environment for Visualizing Images(ENVI,v5.3)对MODIS、Landsat卫星影像进行处理。(1)对MOD02数据进行辐射定标，根据各个波段的中心波长信息得到反射率数据，单位值为0～1无量纲。同时根据MOD03经纬度坐标数据进行图像几何校正，用以矫正非系统因素产生的误差，同时将影像投影到平面上，使其符合地图投影系统。(2)在ENVI软件中利用FLAASH模块对MODIS、Landsat卫星数据进行大气校正，消除大气等因素对地物反射率的影响，获得真实地物反射率等信息。(3)按照公式(27)～公式(30)计算大气可降水量w，0.55nm气溶胶光学厚度AOT₅₅₀，云量CI，地表反照率albedo等参数，再根据公式(8)、(15)计算法向太阳直射I_n和水平表面散射D_h。

复杂地形下太阳辐射估算

借助ENVI软件对将数字高程模型(DEM)的投影转换成与Landsat影像一致的投影系统，提取坡度、坡向与海拔等地形特征参数。按照公式(16)～公式(22)计算复杂地形下森林冠层上方太阳总辐射量。

林下光照强度估算

借助ENVI软件对Landsat影像进行波段计算，按照公式(23)(24)估算出叶面积指数LAI，根据野外调查的20个训练样本数据，建立林下光照强度I_under与消光系数k、叶面积指数LAI、冠层上方光照强度G_i之间的数学关系，用SPSS软件按照公式(25)进行数据拟合，求得研究地1的消光系数k＝0.577(R²＝0.702)，研究地2的消光系数k＝0.449(R²＝0.725)，然后按照公式(26)计算得出整个区域的林下光照强度空间分布。

估算结果：参见图3和图4。

结果验证

借助专业数据统计分析软件Statistical Product and Service Solutions(SPSS，V19)进行林下光照强度遥感估测与野外调查结果(验证样本)的线性拟合精度分析(图5)，利用配对T检验(paired t-test)法比较3种传统估测方法(表1)与野外调查结果的差异，以证明本估测方法的优越性。

表1：3种太阳辐射估测方法

注：以上3种方法都是估测山地环境中森林冠层上方太阳辐射，再结合叶面积指数估测林下光照强度。

表2研究地1验证样本数据的估测方法与野外调查结果比较

表3研究地2验证样本数据的估测方法与野外调查结果比较

林下光照强度差值正态分布显著性Sig.大于0.05，说明利用配对T检验法具有统计学意义。配对T检验的原假设是光照强度差值的分布符合平均值为0的t分布，研究地1数据中，IMBM、DSCT和IARS方法的双侧显著性P都小于0.05，说明3种估测方法与野外调查结果存在显著差异，而估测本方法与野外调查结果显著性差异不明显，研究地2数据分析也能得出相似结果，结果说明，本估测方法估测值比其他方法更加接近野外调查值。从相对均方根误差来看，IMBM、DSCT和IARS相对均方根误差在10～20％之间，而本估测方法的相对均方根误差在8.5％左右，明显高于其他方法的结果。因此，结果表明，本估测方法可以较好地得到地面太阳辐射参数，为区域尺度中林下光照强度的估算提供较为精准的数据支持。

林下光照强度估测的误差来源于三个方面：一是不同种类的林冠层结构和叶面积指数计算误差；二是气溶胶、水汽等大气影响因素在时空尺度上反演精度的问题，三是遥感数据和野外调查结果的空间尺度效应问题。

本估测方法的最大优势在于林下太阳辐射强度的估算仅仅利用了遥感数据(免费)和数字高程模型DEM作为模型输入参数，而野外调查值只用于算法精度验证。