基于多角度卫星观测数据的森林生物量反演方法

1.本发明属于空间信息技术领域，涉及基于多角度卫星观测数据的森林生物量反演方法。


背景技术：

2.森林地上生物量(agb)是量化地表碳循环的关键植被特征，因此被广泛应用于气候变化、生物多样性等研究领域。大尺度的森林agb图被认为是监测植被动态变化的有效手段，因此如何获取精度较高的agb图是遥感科学和生态研究领域的热点研究课题。随着遥感技术的发展，种类繁多的遥感卫星传感器获取的对地观测数据为绘制大尺度森林agb图提供了机会，如光学遥感数据、雷达数据(radar)，和激光雷达数据(lidar)。
3.radar和lidar作为主动遥感技术，由于具备较强的穿透能力，已被广泛用于大尺度森林agb估算。虽然radar和lidar是森林agb监测的有效手段，但这类主动仪器存在一些局限性，如较差的时间和空间覆盖范围，这些问题在短期内无法解决。多角度观测遥感在植被结构参数反演中的能力已经被广泛认可，如已经成功用于植被聚集指数，森林高度等植被结构参数的反演；此外，已有的多角度遥感观测具有长时间序列和全球尺度的观测能力，如中分辨率成像光谱仪(modis)。因此，多角度遥感有望成为森林agb制图的有效手段。


技术实现要素：

4.本发明要解决的问题是在于提供基于多角度卫星观测数据的森林生物量反演方法，解决以上背景技术中存在的问题，本技术基于modis多角度观测中涵盖有植被的结构信息，并且具有长时间序列和全球尺度的观测能力，modis多角度观测数据结合本发明提出的方法能够实现地区或者全球尺度森林地上生物量的制图，绘制更为完整和全面，精度更高。
5.为解决上述技术问题，本发明采用的技术方案是：基于多角度卫星观测数据的森林生物量反演方法，包括以下步骤；
6.步骤s10：基于半变异函数选择出森林覆盖且均质的modis像元；
7.步骤s11：根据步骤s10选择出的森林覆盖且均质的modis像元，采用核驱动模型计算相应的modis像元在近红外波段热点观测方向的反射率；
8.步骤s12：根据步骤s10选择出的森林覆盖且均质的modis像元，对modis像元内的森林进行采样测量，测量获取树木的胸径(dbh)；根据采样的dbh，基于生长异速方程计算样地内树木的agb；最后，将采样得到的所有林木的agb求和作为对应modis像元的地面观测agb。
9.步骤s13：根据步骤s11和步骤s12提供的modis像元在近红外波段的热点反射率和地面实测agb，基于最小二乘回归法建立agb与modis在近红外波段的热点观测方向反射率的数学关系模型，利用该关系模型即可开展相关森林类型的agb预测。
10.在本发明一实施例中，所述步骤s10的具体实现步骤如下：
11.步骤s101：modis像元均质性的评价是基于地统计学中的半变异函数实现。选择晴
空条件下的美国陆地卫星(landsat)地表反照率产品作为半变异函数的输入参数，当基于半变异函数计算的基台值小于5.0e-04时，相关modis像元为均质像元，即像元被森林均匀覆盖。半变异函数的数学表达如下
[0012][0013]
其中，γe(h)代表landsat地表反照率半变异函数估计值，z
xi
代表在位置xi处的landsat地表反照率像元值，z
xi+h
代表距离像元xi为h时的landsat地表反照率像元值，n(h)为所有间距为h的有效像元对数目。本发明中，采样步长h设置为50m。
[0014]
步骤s102：以步骤s101计算得到的半变异函数估计值γ(h)为输入，利用球状模型拟合半方差函数，计算获取基台值。球状模型的具体计算公式如下所示：
[0015][0016]
其中，c0+c为基台值。通过设置变程a，并结合输入的半变异函数估计值γ(h)即可计算获得基台值；本发明中变程a设置为500m。
[0017]
在本发明一实施例中，所述步骤s11的具体实现步骤如下：
[0018]
基于核驱动模型计算modis在近红外波段的热点观测方向反射率的具体过程如下：
[0019]
核驱动模型的通用表达为
[0020][0021]
其中，是λ波段的方向反射率，本发明涉及波段为近红外波段。f
iso
、f
vol
和f
geo
是核驱动模型中决定brdf形状的三个参数，具体表示了同性散射、体散射和几何光学散射所占方向反射率的权重；本发明中，这三个参数由modis官方产品mcd43a1参数产品提供。k
geo
和k
vol
分别是体散射核和几何光学核，这两个核是关于太阳和观测方向的角度函数，描述了太阳天顶角为θ，观测天顶角为相对方位角为时，植被冠层内部的散射现象，包括了多次散射和单次散射。k
geo
和k
vol
的数学表达公式如下所示：
[0022][0023][0024]
其中，是热点函数，c1和c2分别指热点高度和宽度，对于近红外波段，c1＝0.6，c2＝3.3
°
；ξ为太阳光入射方向和观测方向的夹角，其计算公式如下所示：
[0025][0026]
本发明涉及的热点观测方向反射率的几何配置参数为θ＝45
°
，
[0027]
在本发明一实施例中，所述步骤s12的具体实现步骤如下：
[0028]
以树木胸径dbh为输入，基于生长异速方程计算树木agb的具体步骤如下：
[0029]
生长异速方程的数学表达为
[0030]
agb＝exp(b0+b1(ln(dbh)))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0031]
输入测量的dbh就可以获得相应树木的agb；其中，b0和b1为生长异速方程的计算系数，由树种类型决定。本发明中，对于阔叶林：b0＝-2.2094，b1＝2.3867；对于针叶林：b0＝-2.5356，b1＝2.4349。
[0032]
在本发明一实施例中，所述步骤s13的具体实现步骤如下：
[0033]
以步骤s11计算获取的modis热点观测方向反射率和步骤s12提供的modis像元内的地面观测agb数据作为输入，采用最小二乘回归法，得到一个关于agb和modis热点观测之间的一元线性回归关系公式：
[0034]
y＝ax+b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0035]
其中，y为输出的agb，x为输入的modis在近红外波段的热点方向反射率，a和b分别为基于最小二乘法拟合获得的系数和常数项。
[0036]
与现有技术相比，本发明具有的优点和积极效果如下。
[0037]
1、本发明基于modis在近红外波段的热点方向反射率，利用最小二乘法构建森林agb与modis热点反射率的数学关系模型，进而实现基于modis多角度观测数据的森林agb的预测和制图，本发明逻辑清楚，方法简单明了，适应性强，基于modis热点观测反射率获取区域或者全球尺度连续森林agb图的方法，结果可以提高遥感领域相关应用的精度，如需要森林结构信息的生态应用，对研究森林植被定量遥感领域具有重要的研究意义和应用价值；
[0038]
2、在本发明中，利用半变异函数评价相关modis像元的均质性，筛选出被森林覆盖且均质的modis像元，然后，基于核驱动模型(rtclsr)计算森林覆盖且均质modis像元在近红外波段的热点观测方向反射率；最后，基于最小二乘回归法建立森林覆盖且均质modis像元在近红外波段的热点反射率与地面实测agb的一元线性回归关系，得到关于森林agb和modis在近红外波段热点观测之间的数学关系，基于该数学关系即可开展相关森林类型的agb预测，实现基于modis多角度卫星观测数据的森林agb反演。
附图说明
[0039]
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0040]
图1是本发明基于多角度卫星观测数据的森林生物量反演方法的流程示意图；
[0041]
图2是本发明modis像元均质性评价示意图；
[0042]
图3是本发明基于多角度卫星观测数据的森林生物量反演方法的预测生物量精度图。
具体实施方式
[0043]
需要说明的是，在不充突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0044]
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。
[0045]
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明中的具体含义。
[0046]
下面结合附图对本发明的具体实施例做详细说明。
[0047]
中等分辨率成像光谱仪modis是搭载在terra和aqua卫星上的一个重要的传感器，是卫星上唯一将实时观测数据通过x波段向全世界直接广播，并可以免费接收数据并无偿使用的星载仪器，全球许多国家和地区都在接收和使用modis数据。它的主要目的是：实现从单系列极轨空间平台上对太阳辐射、大气、海洋和陆地进行综合观测，获取有关海洋、陆地、冰雪圈和太阳动力系统等信息；进行土地利用和土地覆盖研究、气候的季节和年际变化研究、自然灾害监测和分析研究、长期气候变率和变化以及大气臭氧变化研究等；进而实现对大气和地球环境变化的长期观测和研究的总体(战略)目标。
[0048]
modis最大空间分辨率可达250米，扫描宽度2330公里，modis是当前世界上新一代"图谱合一"的光学遥感仪器，有36个离散光谱波段，光谱范围宽，从0.4微米(可见光)到14.4微米(热红外)全光谱覆盖，modis的多波段数据可以同时提供反映陆地表面状况、云边界、云特性、海洋水色、浮游植物、生物地理、化学、大气中水汽、气溶胶、地表温度、云顶温度、大气温度、臭氧和云顶高度等特征的信息，因此被广泛应用于地表、生物圈、固态地球、大气和海洋的动态监测。
[0049]
如图1-图3所示，本发明涉及一种基于多角度卫星观测数据的森林生物量反演方法，方法的具体流程如下：
[0050]
步骤s10：基于半变异函数选择出森林覆盖且均质的modis像元，如图2；
[0051]
步骤s11：根据步骤s10选择出的森林覆盖且均质的modis像元，采用核驱动模型计算相应的modis像元在近红外波段热点观测方向的反射率；
[0052]
步骤s12：根据步骤s10选择出的森林覆盖且均质的modis像元，对modis像元内的森林进行采样测量，测量获取树木的胸径(dbh)；根据采样的dbh，基于生长异速方程计算样地内树木的agb；最后，将采样得到的所有林木的生物量求和作为对应modis像元的地面观测agb。
[0053]
步骤s13：根据步骤s11和步骤s12提供的modis像元在近红外波段的热点反射率和地面实测agb，基于最小二乘回归法建立agb与modis在近红外波段的热点观测方向反射率的数学关系模型，利用该关系模型即可开展相关森林类型的agb预测。
[0054]
在本发明一实施例中，所述步骤s10的具体实现步骤如下：
[0055]
步骤s101：modis像元均质性的评价是基于地统计学中的半变异函数实现。选择晴空条件下的美国陆地卫星(landsat)地表反照率产品作为半变异函数的输入参数，当基于半变异函数计算的基台值小于5.0e-04时，相关modis像元为均质像元，即像元被森林均匀覆盖。半变异函数的数学表达如下
[0056][0057]
其中，γe(h)代表landsat地表反照率半变异函数估计值，z
xi
代表在位置xi处的landsat地表反照率像元值，z
xi+h
代表距离像元xi为h时的landsat地表反照率像元值，n(h)为所有间距为h的有效像元对数目。本发明中，采样步长h设置为50m。
[0058]
步骤s102：以步骤s101计算得到的半变异函数估计值γ(h)为输入，利用球状模型拟合半方差函数，计算获取基台值。球状模型的具体计算公式如下所示：
[0059][0060]
其中，c0+c为基台值。通过设置变程a结合输入的半变异函数估计值即可计算获得基台值，本发明中变程a设置为500m。
[0061]
在本发明一实施例中，所述步骤s11的具体实现步骤如下：
[0062]
基于核驱动模型计算modis近红外波段的热点观测方向反射率的具体过程如下：
[0063]
核驱动模型的通用表达为
[0064][0065]
其中，是λ波段的方向反射率，本发明涉及波段为近红外波段。f
iso
、f
vol
和f
geo
是核驱动模型中决定brdf形状的三个参数，具体表示了同性散射、体散射和几何光学散射所占方向反射率的权重；本发明中，这三个参数由modis官方产品mcd43a1参数产品提供。k
geo
和k
vol
分别是体散射核和几何光学核，这两个核是关于太阳和观测方向的角度函数，描述了太阳天顶角为θ，观测天顶角为相对方位角为时，植被冠层内部的散射现象，包括了多次散射和单次散射。k
geo
和k
vol
的数学表达公式如下所示：
[0066][0067][0068]
其中，是热点函数，c1和c2分别指热点高度和宽度，对于近红外波段，c1＝0.6，c2＝3.3
°
；ξ为太阳光入射方向和观测方向的夹角，其计算公式如下所示：
[0069][0070]
本发明涉及的热点观测方向反射率的几何配置参数为θ＝45
°
，
[0071]
在本发明一实施例中，所述步骤s12的具体实现步骤如下：
[0072]
以树木胸径dbh为输入，基于生长异速方程计算树木agb的具体步骤如下：
[0073]
生长异速方程的数学表达为
[0074]
agb＝exp(b0+b1(ln(dbh)))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0075]
输入测量的dbh就可以获得相应树木的agb；其中，b0和b1为生长移速方程的计算系数，由树种类型决定。本发明中，对于阔叶林：b0＝-2.2094，b1＝2.3867；对于针叶林：b0＝-2.5356，b1＝2.4349。
[0076]
在本发明一实施例中，所述步骤s13的具体实现步骤如下：
[0077]
以步骤s11计算获取的modis在近红外波段的热点观测方向反射率和步骤s12提供的modis像元内的地面观测agb数据作为输入，采用最小二乘回归法，得到一个关于agb和modis热点观测之间的一元线性回归关系公式：
[0078]
y＝ax+b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0079]
其中，y为输出的agb，x为输入的modis在近红外波段的热点方向反射率，a和b分别位基于最小二乘法拟合获得的系数和常数项。
[0080]
实施例1：
[0081]
在一台配置有英特尔e3-1220@3.1ghz/8m双核处理器的联想thinkstation台式计算机上，以matlab编程语言作为实现工具进行了案例实施，也可以其他控制程序实现，不做具体的限定；在该实施例中，以美国东海岸的巴特利特实验森林、哈佛森林、豪兰研究森林、哈伯德布鲁克实验森林和佩诺布斯科特实验森林为研究区，根据图1所示的方法实施流程，对本发明方法进行了评价验证。
[0082]
图3是评价结果，从结果可以看出modis热点观测方向获得的反射率与森林生物量具有很强的相关性，能够解释高达63％的生物量变化。因此，本专利提出的基于多角度卫星观测数据的森林生物量反演方法有效可行。
[0083]
以上对本发明的一个实施例进行了详细说明，但所述内容仅为本发明的较佳实施例，不能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖范围之内。