一种基于3D温度特征的作物水分胁迫检测方法

一种基于3d温度特征的作物水分胁迫检测方法
技术领域
1.本发明属于作物水分胁迫检测领域，特别是涉及一种基于3d温度特征的作物水分胁迫检测方法。


背景技术：

2.气候变化和水资源短缺极大地威胁着全球粮食安全，并限制了作物生产，因此，需要对整个生长季节的水分状况进行有效监测，以便有效地最大限度地提高作物的产量和质量。植物对水分亏缺的直接反应是关闭叶片气孔，导致冠层温度升高。作物表面温度的检测对于及时评估水分胁迫具有重要意义。近年来，热红外相机广泛应用于作物的水分胁迫检测，尽管热图像包含作物和背景信息，但由于散射分布和微小的温差，识别作物和提取其表面温度很困难。因此，一些研究人员将热图像与其他图像(如颜色和光谱图像)相结合来识别作物冠层。然而，由于缺乏精确的图像配准方法，将热彩色成像系统应用于作物冠层温度提取的有效性受到了限制。并且，目前大多数研究都集中在作物二维温度特性的测量上，例如，中国发明专利申请“一种农田表面温度监测数据精量甄别方法及系统”(cn113748867a)只对作物冠层多点温度或整个冠层温度获取，没能从作物温度的空间分布差异性考虑，即对作物水分胁迫状态在何位置检测方面的研究很少；中国发明专利申请“基于冠层温度的作物水分胁迫诊断方法及系统”(cn110020793a)没能从作物在时序上的生理活动变化考虑，具体的对作物水分胁迫状态在何时检测方面的研究很少；此外，水分胁迫条件下作物的温度特性与冠层结构有关，获取作物的三维信息，这有助于将温度、叶片形状和作物冠层结构结合起来，作为综合分析的输入。


技术实现要素：

3.针对以上技术问题，本发明的目的是提供一种基于3d温度特征的作物水分胁迫检测方法，利用摄影测量的技术将可见光图像与热红外图像精配准，利用作物不同空间分布的叶片在时序上的温度变化差异确定作物水分胁迫指数在何位置检测。
4.为了实现上述目的，本发明提供了如下技术方案：
5.一种基于3d温度特征的作物水分胁迫检测方法，包括如下步骤：
6.s1、热红外相机、双目可见光相机立体标定；
7.热红外相机和双目可见光相机同时分别拍摄多张不同摆放位置的棋盘格图像，并将棋盘格图像分为双目可见光相机左、右摄像头标定组和双目可见光相机左摄像头与热红外相机标定组，标定获得热红外相机和双目可见光相机的标定参数；所述标定参数包括双目可见光相机左摄像头的内部参数矩阵、双目可见光相机右摄像头的内部参数矩阵、热红外相机的内部参数矩阵、双目可见光相机左摄像头与双目可见光相机右摄像头之间的外部参数、热红外相机与双目可见光相机左摄像头之间的外部参数、双目可见光相机左摄像头的畸变参数、双目可见光相机右摄像头的畸变参数以及热红外相机的畸变参数；其中，
8.s2、热红外图像和可见光图像的粗配准；
9.s2.1、通过热红外相机和双目可见光相机垂直于地面对作物进行图像采集，连续获取多组同一株作物一天中不同时间段的包含湿参考面的左、右可见光图像和热红外图像；
10.s2.2、对每组左、右可见光图像和热红外图像进行处理，根据步骤s1获得的标定参数分别对左、右可见光图像进行立体校正、消除畸变和行对准，获得校正后的左、右可见光图像；
11.s2.3、采用surf算法分别检测各组校正后的左、右可见光图像的特征点；运用三角测量原理计算出左、右可见光图像之间相对应的特征点对的世界坐标，并通过公式1将特征点对的世界坐标反投影到热红外图像上；
[0012][0013]
式中，s为非零尺度因子，即由特征点对计算的世界坐标的z分量；coordinates
thermal
为热红外图像上反投影点的齐次坐标；k为步骤s1获得的热红外相机的内部参数矩阵；r为步骤s1获得的热红外相机相对于双目可见光相机左摄像头的三维旋转矩阵；t为步骤s1获得的热红外相机相对于双目可见光相机左摄像头的平移向量；coordinates
world
是特征点对计算的世界坐标的齐次坐标；
[0014]
s2.4、使用ransac算法计算出左可见光图像的特征点与热红外图像上的反投影点之间的单应性矩阵，完成所有组的热红外图像与左可见光图像粗配准；
[0015]
s3、热红外图像和可见光图像的精配准；
[0016]
s3.1、使用ransac算法计算左可见光图像中的特征点与热红外图像中经过偏差补偿的反投影点之间单应性矩阵，对热红外图像与左可见光图像进行重新配准；
[0017]
其中，δ
x
，δy为热红外图像中的反投影点坐标与热红外图像中真实对应点的坐标(x
0 thermal
，y
0 thermal
)之间的位置偏差，初始值为(0，0)；
[0018]
s3.2、按一定步长增加位置偏差值，并重复步骤s3.1，直至结构相似性指数ssim的值达到最高，获得最优单应性矩阵，完成所有组的热红外图像与可见光图像精配准；所述结构相似性指数ssim通过公式2计算获得；
[0019][0020]
式中，x、y分别为热红外图像和左可见光图像；u
x
、uy分别表示图像x，y的均值；分别表示图像x，y的方差；σ
xy
表示图像x，y之间的协方差；c1、c2为常数；
[0021]
s4、双目三维重建生成点云并与温度数据融合；
[0022]
s4.1、根据步骤s2.2中同组的校正后的左、右可见光图像，利用立体匹配算法生成视差图；然后，利用相似三角形原理结合双目可见光相机左可见光图像的rgb颜色值，通过公式3获得作物冠层的原始三维点云数据；
[0023][0024]
式中，z，x，y为点云三维坐标；d为视差图的深度值；b为双目可见光相机的基线，单位为mm，f为双目可见光相机左摄像头的焦距，单位为像素；d是视差图的视差值；(u，v)是左可见光图像的像素坐标；(u0，v0)是双目相机左摄像头图像平面的主点像素坐标；
[0025]
s4.2、通过公式4将步骤s4.1获得作物冠层的原始三维点云数据与同组的热红外图像中的温度数据进行融合，使作物冠层的三维点云数据既包含rgb颜色信息又包含温度信息；
[0026]crgb
＝hc
ir
ꢀꢀꢀ
公式4
[0027]
式中，c
rgb
为双目可见光相机左摄像头的像平面坐标；h为步骤s3获得的双目可见光相机左可见光图像与热红外图像之间的最优单应性矩阵；c
ir
为热红外相机的像平面坐标；
[0028]
s5、确定作物水分胁迫指数的最优检测位置；
[0029]
s5.1、通过统计滤波滤除步骤s4获得的融合后的作物冠层的三维点云数据的离群点，然后采用基于hsv色彩空间的颜色模型分割出绿色冠层点云；通过meanshift算法在绿色冠层点云中寻找各冠层叶片质心，而后以质心为初始点进行区域生长，聚类分割得到每个冠层叶片点云并标记；
[0030]
s5.2、以拍摄时间为横坐标，每个叶片的平均温度为纵坐标，绘制同一植株的每个冠层叶片的时间——温度变化曲线；同时，绘制同组数据中所获得的湿参考面的时间——温度变化曲线；找出与湿参考面的时间——温度变化曲线最为相近的n个冠层叶片作为最能代表整株作物水分亏缺状况的叶片，通过公式6计算这n个冠层叶片的平均温度，并将其作为作物冠层平均温度；
[0031][0032]
其中，t
canopy_ave
为最能代表冠层蒸腾作用的n个冠层叶片的平均温度，单位为℃；n为叶片个数；ti为第i个叶片的平均温度，单位为℃；
[0033]
s6、计算作物水分胁迫指数cwsi值；
[0034]
通过公式7作物水分胁迫指数：
[0035][0036]
式中，t
crop
为步骤s5获得的作物冠层平均温度，单位为℃；t
dry
为干参考面温度，单位为℃；t
wet
为湿参考面温度，单位为℃。
[0037]
所述步骤s1中，
[0038]
所述双目可见光相机左摄像头的内部参数矩阵包括左摄像头焦距和左摄像头成像原点坐标；
[0039]
所述双目可见光相机右摄像头的内部参数矩阵包括右摄像头焦距和右摄像头成像原点坐标；
[0040]
所述热红外相机的内部参数矩阵包括热红外相机焦距和热红外相机成像原点坐标；
[0041]
所述双目可见光相机左摄像头与双目可见光相机右摄像头之间的外部参数包括双目可见光相机右摄像头相对于双目可见光相机左摄像头的三维旋转矩阵和平移向量，代表了双目可见光相机的右摄像头与左摄像头之间的转换关系；
[0042]
所述热红外相机与双目可见光相机左摄像头之间的外部参数包括热红外相机相对于双目可见光相机左摄像头的三维旋转矩阵和平移向量，代表了热红外相机与双目可见光相机左摄像头之间的转换关系；
[0043]
所述双目可见光相机左摄像头的畸变参数、双目可见光相机右摄像头的畸变参数以及热红外相机的畸变参数均包括径向畸变参数和切向畸变参数，径向畸变参数用于修正双目可见光相机摄像头由于光学透镜的特性使得成像存在的径向畸变，切向畸变参数用于修正由于装配方面的误差使得传感器与光学镜头之间非完全平行造成的切向畸变。
[0044]
所述步骤s2中，以30分钟为时间间隔，在11:00-15:00时间段内获取9组同一株作物的包含湿参考面的左、右可见光图像和热红外图像。
[0045]
所述步骤s5.2中，n大于等于总冠层叶片数的1/3。
[0046]
与现有技术相比，本发明的有益效果在于：
[0047]
本发明利用热红外相机获取作物的温度图像和温度矩阵，利用双目可见光相机获取作物的左、右可见光图像，进行双目立体视觉匹配并计算出可见光双目三维点云；通过摄影测量的方法将热红外图像与双目可见光相机左摄像头拍摄的可见光图像进行精配准；然后将经过单应变换的热红外相机的温度信息与双目可见光相机的三维点云进行信息融合，生成作物的3d温度特征；并通过基于区域生长的meanshift聚类分割算法分割出叶片；最后利用作物不同空间分布的叶片在时序上的温度变化差异确定作物水分胁迫指数检测的最优位置。
附图说明
[0048]
图1为基于3d温度特征的作物水分胁迫检测方法流程图；
[0049]
图2a为双目可见光相机拍摄的棋盘格图像；
[0050]
图2b为热红外相机拍摄的棋盘格图像；
[0051]
图3为双目可见光相机和热红外相机标定流程图；
[0052]
图4为图像粗配准工作流程框图；
[0053]
图5为图像精配准工作流程框图；
[0054]
图6为作物3d温度特征获取过程框图；
[0055]
图7a～7d为双目可见光相机左可见光图像与热红外图像粗配准过程示意图；
[0056]
图8a～8b为双目可见光相机左可见光图像与热红外图像精配准过程示意图；
[0057]
图9a～9b为基于hsv色彩空间的颜色模型分割出的冠层区域点云；
[0058]
图10为采用基于区域生长的meanshift聚类算法分割出的叶片点云示意图。
具体实施方式
[0059]
下面结合附图和实施例对本发明进行进一步说明。
[0060]
如图1所示，一种基于3d温度特征的作物水分胁迫检测方法，包括如下步骤：
[0061]
s1、热红外相机、双目可见光相机立体标定；
[0062]
热红外相机和双目可见光相机同时分别拍摄多张不同摆放位置的棋盘格图像，并将棋盘格图像分为双目可见光相机左、右摄像头标定组和双目可见光相机左摄像头与热红外相机标定组，标定获得热红外相机和双目可见光相机的标定参数；所述标定参数包括双目可见光相机左摄像头的内部参数矩阵、双目可见光相机右摄像头的内部参数矩阵、热红外相机的内部参数矩阵、双目可见光相机左摄像头与双目可见光相机右摄像头之间的外部参数、热红外相机与双目可见光相机左摄像头之间的外部参数、双目可见光相机左摄像头的畸变参数、双目可见光相机右摄像头的畸变参数以及热红外相机的畸变参数；其中，
[0063]
所述双目可见光相机左摄像头的内部参数矩阵包括左摄像头焦距和左摄像头成像原点坐标；
[0064]
所述双目可见光相机右摄像头的内部参数矩阵包括右摄像头焦距和右摄像头成像原点坐标；
[0065]
所述热红外相机的内部参数矩阵包括热红外相机焦距和热红外相机成像原点坐标；
[0066]
所述双目可见光相机左摄像头与双目可见光相机右摄像头之间的外部参数包括双目可见光相机右摄像头相对于双目可见光相机左摄像头的三维旋转矩阵和平移向量，代表了双目可见光相机的右摄像头与左摄像头之间的转换关系；
[0067]
所述热红外相机与双目可见光相机左摄像头之间的外部参数包括热红外相机相对于双目可见光相机左摄像头的三维旋转矩阵和平移向量，代表了热红外相机与双目可见光相机左摄像头之间的转换关系。
[0068]
所述双目可见光相机左摄像头的畸变参数、双目可见光相机右摄像头的畸变参数以及热红外相机的畸变参数均包括径向畸变参数和切向畸变参数，径向畸变参数用于修正双目可见光相机摄像头由于光学透镜的特性使得成像存在的径向畸变，切向畸变参数用于修正由于装配方面的误差使得传感器与光学镜头之间非完全平行造成的切向畸变。
[0069]
立体标定需要用到一个棋盘格和卤素灯。棋盘格的黑色方块用黑色的乙烯胶带覆盖着，由于它对红外辐射吸收更高，因此通过卤素灯照射，其温度上升更快，使得棋盘格在热红外图像中成像效果更好。使用热红外相机和双目可见光相机同时拍摄不同摆放位置的棋盘格照片，如图2a和图2b所示，共拍摄20-30组图像以提供冗余标定。使用matlab r2018a软件中的stereo camera calibrator工具箱对相机进行双目标定，如图3所示。
[0070]
s2、热红外图像和可见光图像的粗配准；
[0071]
如图4所示，基于摄影测量原理的图像配准算法首先需要计算出双目可见光图像中的特征点对所对应的世界坐标，并计算出世界坐标在热红外图像上的反投影点的坐标；然后计算出左可见光图像中的特征点和热红外图像中的反投影点之间的最优单应变换。
[0072]
s2.1、通过热红外相机和双目可见光相机垂直于地面对作物进行图像采集，连续获取多组同一株作物一天中不同时间段的包含湿参考面的左、右可见光图像和热红外图像；
[0073]
本发明实施例以30分钟为时间间隔，在11:00-15:00时间段内获取9组同一株作物的包含湿参考面的左、右可见光图像和热红外图像。
[0074]
s2.2、对每组左、右可见光图像和热红外图像进行处理，根据步骤s1获得的标定参数分别对左、右可见光图像进行立体校正、消除畸变和行对准，获得校正后的左、右可见光图像，如图7a所示；
[0075]
s2.3、采用surf算法分别检测各组校正后的左、右可见光图像的特征点，如图7b所示；运用三角测量原理计算出左、右可见光图像之间相对应的特征点对的世界坐标，并通过公式1将特征点对的世界坐标反投影到热红外图像上，如图7c所示；
[0076][0077]
式中，s为非零尺度因子，即由特征点对计算的世界坐标的z分量；coordinates
thermal
为热红外图像上反投影点的齐次坐标；k为步骤s1获得的热红外相机的内部参数矩阵；r为步骤s1获得的热红外相机相对于双目可见光相机左摄像头的三维旋转矩阵；t为步骤s1获得的热红外相机相对于双目可见光相机左摄像头的平移向量；coordinates
world
是特征点对计算的世界坐标的齐次坐标；
[0078]
s2.4、使用ransac算法计算出左可见光图像的特征点与热红外图像上的反投影点之间的单应性矩阵，完成所有组的热红外图像与左可见光图像粗配准，如图7d所示。
[0079]
s3、热红外图像和可见光图像的精配准；
[0080]
使用单应变换对两幅图像进行配准时，假设目标场景接近平面。然而，作物的表面结构非常复杂，深度变化很大。因此，当需要更高的配准精度时，如图5所示，需要在全局粗配准后采取进一步的措施。
[0081]
s3.1、使用ransac算法计算左可见光图像中的特征点与热红外图像中经过偏差补偿的反投影点之间单应性矩阵，对热红外图像与左可见光图像进行重新配准；
[0082]
其中，δ
x
，δy为热红外图像中的反投影点坐标与热红外图像中真实对应点的坐标(x
0 thermal
，y
0 thermal
)之间的位置偏差，初始值为(0，0)；如图8a所示，左可见光图像与热红外图像存在偏差。
[0083]
s3.2、按一定步长增加位置偏差值，并重复步骤s3.1，直至结构相似性指数ssim的值达到最高，获得最优单应性矩阵，完成所有组的热红外图像与可见光图像精配准，配准效果如图8b所示；所述结构相似性指数ssim通过公式2计算获得；
[0084][0085]
式中，x、y分别为热红外图像和左可见光图像；u
x
、uy分别表示图像x，y的均值；分别表示图像x，y的方差；σ
xy
表示图像x，y之间的协方差；c1、c2为常数。
[0086]
s4、双目三维重建生成点云并与温度数据融合。
[0087]
如图6所示，生成作物的3d温度特征，需要将热红外图像的温度信息与双目可见光图像生成的三维点云数据进行信息融合。
[0088]
s4.1、根据步骤s2.2中同组的校正后的左、右可见光图像，利用立体匹配算法生成视差图；然后，利用相似三角形原理结合双目可见光相机左可见光图像的rgb颜色值，通过公式3获得作物冠层的原始三维点云数据；
[0089][0090]
式中，z，x，y为点云三维坐标；d为视差图的深度值；b为双目可见光相机的基线(即左、右摄像头之间的距离)，单位为mm，f为双目可见光相机左摄像头的焦距，单位为像素；d是视差图的视差值；(u，v)是左可见光图像的像素坐标；(u0，v0)是双目相机左摄像头图像平面的主点像素坐标。
[0091]
s4.2、通过公式4将步骤s4.1获得作物冠层的原始三维点云数据与同组的热红外图像中的温度数据进行融合，使作物冠层的三维点云数据既包含rgb颜色信息又包含温度信息；
[0092]crgb
＝hc
ir
ꢀꢀꢀ
公式4
[0093]
式中，c
rgb
为双目可见光相机左摄像头的像平面坐标；h为步骤s3获得的双目可见光相机左可见光图像与热红外图像之间的最优单应性矩阵；c
ir
为热红外相机的像平面坐标；
[0094]
s5、确定作物水分胁迫指数的最优检测位置；
[0095]
s5.1、通过统计滤波滤除步骤s4获得的融合后的作物冠层的三维点云数据的离群点，然后采用基于hsv色彩空间的颜色模型分割出绿色冠层点云(如图9a所示为冠层的颜色点云，如图9b所示为冠层的温度点云)；通过meanshift算法在绿色冠层点云中寻找各冠层叶片质心，而后以质心为初始点进行区域生长，聚类分割得到每个冠层叶片点云并标记，如图10所示；
[0096]
叶片是植物进行蒸腾作用等生理过程的主要载体，是和外界环境进行物质交换的重要器官，能够准确定量的描述叶片对作物生长检测和精准管理具有重要意义。区域生长算法思想为通过判断点云法线之间的夹角和曲率是否满足阈值决定近邻点与当前点是否为一类，若法线夹角和曲率均小于阈值则作为种子点继续循环迭代，若仅满足法线夹角阈值则只归类，邻域内不再有种子点时，分类完成，并从原始点云中去除，重新选择初始种子点开始归类。meanshift聚类算法即均值漂移算法，是通过点云密度梯度上升追踪样本点附近密度中心，从而根据密度中心进行归类的算法。
[0097]
s5.2、以拍摄时间为横坐标，每个叶片的平均温度为纵坐标，绘制同一植株的每个冠层叶片的时间——温度变化曲线；同时，绘制同组数据中所获得的湿参考面的时间——温度变化曲线；找出与湿参考面的时间——温度变化曲线最为相近的n个冠层叶片作为最能代表整株作物水分亏缺状况的叶片，通过公式6计算这n个冠层叶片的平均温度，并将其作为作物冠层平均温度；
[0098][0099]
其中，t
canopy_ave
为最能代表冠层蒸腾作用的n个冠层叶片的平均温度，单位为℃；n为叶片个数；ti为第i个叶片的平均温度，单位为℃。
[0100]
优选地，n大于等于总冠层叶片数的1/3。
[0101]
已知对植株温度起直接作用的是植株蒸腾，在相同外界环境条件下，作物不同叶位的叶片由于叶绿素含量和叶片功能期长短差异造成了蒸腾作用强度也不同。蒸腾作用强的叶片带走较多热量，叶片温度较低且变化平缓和湿参考面的温度变化趋势相近；相反蒸腾作用弱的叶片温度较高且变化较快。
[0102]
s6、计算作物水分胁迫指数cwsi值；
[0103]
通过公式7作物水分胁迫指数：
[0104][0105]
式中，t
crop
为步骤s5获得的作物冠层平均温度，单位为℃；t
dry
为干参考面温度，单位为℃；t
wet
为湿参考面温度，单位为℃。
[0106]
所述湿参考面的制作方法为：将一个2cm厚的聚苯乙烯泡沫板漂浮在尺寸为15cm
×
10cm
×
15cm的装有水的塑料容器中，使其覆盖大部分水面；泡沫板表面包裹一层0.5mm厚的吸水性无纺布和粘胶混合布，外层覆盖有2mm厚的聚酯无纺吸水布，它的作用是吸收水分保持表面湿润，模拟作物完全蒸腾作用下的叶片表面温度即t
wet
。t
dry
模拟作物气孔完全关闭不进行蒸腾作用时的温度，其采用空气温度加5℃的方法计算。