一种基于光合概念模型的作物生物量估测方法

本发明属于无人机遥感监测及图像处理，更具体地，涉及一种基于光合概念模型的作物生物量估测方法。

背景技术：

1、地上生物量（above ground biomass，agb）作为作物生态系统监测中的重要指标之一，不仅与作物的生长发育及产量形成等密切相关，同时也是开展全球气候变化、碳循环、物质流动和能量交换等研究的重要内容，其时序变化能够直接表征作物进行净光合作用积累有机物的能力，反映作物的生长态势和营养状况。基于agb的长势评估是对作物生长状况及其对环境条件响应的定量分析，快速且精准地获取作物agb信息对营养诊断、田间施肥等精细化管理和农业生态环境保护具有科学意义。

2、传统人工agb测量方法需要进行破坏性采样，然后将样本进行清洗、分解、烘干、称重和记录等操作，不仅费时、费力、效率低下，且人的主观性差异对数据结果有较大的影响，无法应用于田间尺度的作物长势监测。遥感技术以其快捷、无损、覆盖范围广的优势，成为了大范围agb估测的重要手段。而在所有的遥感平台中，无人机遥感平台以其简单易操作、可搭载多种传感器、可自由设定数据获取间隔等特点，尤其受到精准农业研究和从业者的青睐。

3、现有研究表明，通过若干波段反射率数学组合计算的植被指数是小麦、玉米、水稻等不同作物agb估测的有效方法。作物叶面积指数（leaf area index，lai）、agb等长势参数与光谱植被指数间存在较强的相关关系，但可能受到作物种类、生长阶段、环境因子等多种因素的影响导致关系不稳定。因此，有必要采用遥感技术构建时-空-谱特征深度融合的作物通用agb估测模型，对直接影响粮食安全、经济收益和生态环境的大田作物进行时空精细化长势监测。

技术实现思路

1、针对相关技术中的上述技术问题，本公开提出一种基于光合概念模型的作物生物量估测方法，能够克服现有技术存在的上述不足。

2、为实现上述技术目的，本公开的技术方案是这样实现的：

3、本公开目的在于提供一种基于光合概念模型的作物生物量估测方法，包括以下步骤：

4、s1数据采集与处理：根据试验区需要分析的作物情况，规划无人机图像采集方案，执行采集方案来获取作物的冠层遥感影像，所述冠层遥感影像包括rgb影像和多光谱影像；对rgb影像进行图像拼接及三维重建，对多光谱影像进行几何处理以及辐射定标处理，然后利用rgb影像和多光谱影像的不同波段计算灰度共生矩阵纹理；

5、s2提取光谱与几何信息：

6、s2.1提取作物冠层高度：根据rgb影像预处理得到的dsm，获取无植被覆盖数字地形模型dtm，利用数字表面模型dsm和无植被覆盖数字地形模型dtm叠加相减以获取作物冠层高度模型，从而得出作物的冠层高度；根据rgb影像中田块的位置勾画rgb感兴趣区域，通过rgb感兴趣区域来提取计算出田块作物的冠层高度均值，将冠层高度均值作为田块的冠层高度；

7、s2.2提取冠层反射率信息：对于预处理完成的多光谱影像，根据多光谱影像中田块的位置勾画多光谱感兴趣区域，通过多光谱感兴趣区域来提取计算田块的反射率均值，将反射率均值作为田块的反射率；

8、s2.3计算植被指数：根据多光谱影像提取田块的光谱反射率，计算植被指数；

9、s3构建与简化作物agb估测通用模型：

10、s3.1构建agb估测理论模型：根据净同化率（net assimilation rate，nar）、生物量、lai等参数的物理含义来推导构建作物agb的计算模型；

11、s3.2挖掘作物nar的遥感替代指标：

12、s3.2.1:首先，利用遥感技术获取表征作物某种生长状况的植被指数、覆盖度、纹理、冠层高度、冠层体积等，其中：植被指数用于表征作物色素含量，覆盖度用于表征作物生长密度，纹理用于表征作物群体之间或作物与土壤背景之间的相互关系，冠层高度和冠层体积用于表征作物在不同方向上的体量；

13、s3.2.2：然后，通过获取不同时期以上参数的时序变化趋势，进一步提取各参数的变化率特征，形成nar的可替代遥感指标变量库，并以变量库中作物全生育期的指标和实测lai为数据基础，与实测agb进行对比，通过比较agb的计算精度，挖掘精度较高时的替代变量，从而反推出可用于替代作物nar的遥感指标；

14、s3.3估测几何光谱信息一体化的叶面积指数：根据叶面积密度的物理含义推导构建作物lai的计算模型，从而得到作物的lai；

15、s3.4估测光合概念模型的作物agb值：利用s3.2得到的nar和s3.3得到的lai，并根据agb计算公式得到作物不同时期的agb值，从而参量化生物量的理论计算模型；

16、s3.5简化模型：在无人机数据可获取的前提下，选取处于最大生物量时的时间，采用图像分割方法将s3.4所取得的agb值随时间变化曲线分割为三角形与梯形，选取若干个重要时间的数据代入即可完成多生育期agb的估测，重要时间为作物出苗、生物量最大以及需要监测时的播种后天数，具体的重要时间在控制实验中获得。

17、优选地，所述s1具体包括以下过程：

18、s1.1：首先根据作物种植面积、地形地貌、天气情况、数据获取要求规划无人机图像采集方案，所述规划无人机图像采集方案包括无人机的航拍方式、航拍路线和设置的航带重叠度；

19、s1.2：然后根据所述无人机图像采集方案，利用带有rgb和多光谱传感器的无人机系统获取作物的冠层遥感影像，所述冠层遥感影像包括rgb影像和多光谱影像；

20、s1.3：对rgb影像进行图像拼接及三维重建，从而生成数字表面模型（digitalsurface model，dsm）和数字正射影像（digital orthophoto map，dom）；对多光谱影像进行几何处理以及辐射定标处理，通过几何处理对每个镜头进行物镜畸变矫正，再对不同波段进行光学配准，对不同波段的影像进行辐射定标；利用rgb和多光谱影像的不同波段计算灰度共生矩阵的纹理，所述纹理包括均值、方差、同质性、对比度、相异性、熵、二阶矩、相关性。

21、优选地，所述s1.3中的图像拼接及三维重建具体方式为：

22、首先，对齐rgb影像照片（align photos）构建稀疏点云：基于人工标注的地面控制点，采用sift（scale-invariant feature transform）算子提取影像特征点完成初步影像间的匹配；基于特征点对计算影像在三维空间的相对位置，并解算相机的姿态参数，使用光束法区域网平差进行优化并构建出稀疏点云；

23、其次，构建密集点云：根据稀疏点云以及地面控制点坐标，采用多视角立体算法，将稀疏点云中的点作为种子点，进一步进行加密过滤，从而生成密集点云；

24、最后，生成模型：基于密集点云，通过加密和网格化插值生成数字表面模型（digital surface model，dsm）和数字正射影像（digital orthophoto map，dom）。

25、优选地，所述s1.3中，优选采用 agisoft photoscan professional v1.4.5（agisoft llc, st. petersburg, russia）软件对rgb影像进行图像拼接及三维重建。

26、优选地，所述s1.3中，辐射定标处理即将mca 影像中数据dn（digital number）值转换成具有实际物理意义的地物表面反射率值，辐射定标采用统一分段经验线性法完成；

27、采用8块已知反射率的定标毯（反射率分别为3%、6%、12%、24%、36%、48%、56%和80%）来进行辐射定标过程，首先根据定标毯勾画感兴趣区域，统计感兴趣区域内的平均nd值，再结合已知反射率完成辐射定标过程，具体计算公式如下：

28、；

29、上述公式一中，y是辐射定标后的反射率，x是不同波段的dn值，i是从1到12的波段数，常数a和b是所用定标毯dn值和反射率拟合方程的相应参数。

30、优选地，所述s1.3中，纹理的计算过程中，所述感兴趣区域内的所有像素点的纹理值均值可以作为对应田块的纹理值。取均值方法可以和反射率、冠层高度等相一致，画的区域内的所有点的均值就作为感兴趣区域所对应田块的该指标值。此处可采用公知常识或现有技术。

31、优选地，所述s2.3具体包括以下过程：

32、s2.3.1提取光谱反射率：根据多光谱影像提取田块区域性的光谱反射率；

33、s2.3.2计算植被指数：按照下列过程具体计算植被指数：

34、归一化差异植被指数ndvi = (r800－r670)/(r800+r670)，

35、归一化红边差异植被指数ndre = (r800－r720)/(r800+r720)，

36、红边叶绿素指数cired edge = r800/r720－1，

37、绿光叶绿素指数cigreen = r800/r550－1，

38、可见光大气阻抗植被指数vari = (r550－r670)/(r550+r670－r490)；

39、其中，r490、r550、r670、r720、r800分别代表多光谱影像的蓝色波段、绿色波段、红色波段、红边波段、近红外波段的反射率。

40、优选地，所述s3.1具体包括以下过程：

41、s3.1.1定义参数：nar用于反映绿色植物的光合作用效率，nar的定义为作物单位叶面积在单位时间内的干物质积累量，单位为g/cm2/day；生物量主要指地上部生物量，生物量的定义为植被在一定时间内单位地面面积所产有机质的累积量（即w/a），单位为g/cm2；叶面积指数又叫叶面积系数，是指单位土地面积上植物叶片总面积占土地面积的倍数（即la/a），该参数无量纲；

42、s3.1.2：根据定义，nar可以表示为：

43、；

44、其中，w为干物质积累量，la表示叶面积，lai表示叶面积指数，agb表示地上生物量，a表示占地面积，t表示时间；

45、s3.1.3：上述公式二经过化简可以表示为：

46、；

47、s3.1.4：由此可以得到以nar和lai为桥梁的agb参量化计算模型为：

48、；

49、s3.1.5：则根据微积分定义可以将agb表达为：

50、；

51、s3.1.6：在时间t不连续的情况下，可以将上式近似表达为：

52、；

53、基于上述推导，即可构建不连续时间条件下作物agb计算的机理模型，在获取数据时间次数较多的情况下，可以精确计算作物agb。但因为模型中的nar无法直接获取，lai的田间测量十分费力，因此下一步的重点为利用遥感技术精确获取nar和lai。

54、优选地，所述s3.3具体包括以下过程：

55、s3.3.1：叶面积密度（leaf area density，lad）用于表征植被垂直方向上不同分层内叶面积差异的参数，其定义为单位地面面积、单位高度范围内的单面叶面积总和，其在高度上的积分为叶面积指数；

56、根据定义，lai的计算公式可以写为：

57、；

58、公式七中，h为作物冠层高度变量，h为作物冠层的最大高度；

59、s3.3.2：为简化计算过程，在不同时期，均可以采用不同的二次函数对叶面积进行拟合，拟合方程如下：

60、；

61、其中，h为作物冠层高度变量，a1, a2, a3为二次函数中的系数和常数；

62、s3.3.3：将lad模拟的二次函数带入lai计算公式中可以得到：

63、；

64、s3.3.4：根据不定积分求解可以得到：

65、；

66、s3.3.5：为进一步简化以上计算模型，引入最大叶面积密度 ladmax，假设 ladmax在h/n处获得，则对于二次函数，取最值时有：

67、；

68、s3.3.6：根据上述公式十一即可得到a1，a2和h之间的关系：

69、；

70、s3.3.7：将上述公式十二代入二次函数即可得到函数的最大值为：

71、；

72、s3.3.8：将a1, a2和h之间的关系带入不定积分求解得到的lai式中，可以得到lai的计算公式为：

73、；

74、s3.3.9：结合以上推导可以得到lai与 ladmax和h之间的简化关系为：

75、；

76、其中，α和β均为常数；

77、因此，最后得到的lai计算模型中，公式的变量为 ladmax和h；

78、其中，h可以通过无人机获取的rgb数据在拼接时生成的dsm和dtm相减得到； ladmax与植被覆盖度密切相关，而植被指数可以用以估测植被覆盖率，因此，从这一方面分析，植被指数可以用来表征 ladmax；

79、s3.3.10：最终lai的计算可以等价与植被指数和最大冠层高度的乘积，即：

80、；

81、优选地，所述s3.5具有包括以下过程：

82、s3.5.1：因一般作物的agb在整个生育期的变化趋势为先增加后降低，因此，在最大生物量时的时间以及对应的无人机数据可获取的情况下，采用图像分割方法将agb随时间变化曲线分割为三角形与梯形，因此生物量的计算模型可以简化为：

83、；

84、其中，si和sii分别为分割图像中三角形和梯形的面积，ts为生物量最大时对应的时间；

85、s3.5.2：代入相应参数后可以得到：

86、；

87、其中，pnar为nar的有效遥感替代指标，t0为生物量大于0的播种后天数，ts为生物量最大时对应的播种后天数，pnars为ts时期对应的nar有效遥感替代指标，vis为ts时期对应的植被指数，hs为ts时期对应的冠层高度；

88、s3.5.3：因此，只需要获取几个重要时间的数据即可完成多生育期agb的估测，重要的时间为作物出苗、生物量最大以及需要监测时的播种后天数，这些具体的时间可以在控制实验中获得。

89、优选地，所述作物包括水稻、小麦、玉米等。

90、本公开有益效果：本公开通过挖掘作物叶面积密度的有效光谱替代指标，结合作物叶面积在空间上的几何分布规律，使得基于物理定义的叶面积指数计算模型转化为遥感特征参量化模型，构建融合遥感光谱和几何特征的作物叶面积指数估测模型，解决光谱指数在中高植被覆盖度时饱和以及叶面积指数被低估的问题，充分融合作物垂直分布特征，提升叶面积指数的估测精度；

91、本公开根据作物净同化率的有效遥感替代指标和几何光谱信息融合估测的叶面积指数结果，通过物理过程定义的机理模型和遥感估测的经验模型相结合，构建多时期半经验半机理的作物生物量估测光合概念模型，并进行模型简化，可便捷获取田间作物生长态势的时空精细化分布状况，解决作物生物量估测模型精度不稳定和迁移性差的问题。