一种基于半球图像法测量叶面积指数的方法与流程

本发明属于生态学研究技术领域，更为具体地讲，涉及一种基于半球图像法测量叶面积指数的方法。

背景技术：

叶面积指数是一个无量纲量，可认为叶面积与地面积之间的比值。它反映了叶片的数量和分布情况，是研究冠层中光线穿透情况，冠层生产力，冠层下土壤水分蒸发、蒸腾损失总量，冠层截留，及土壤温度的基础因子。

目前叶面积指数的测量方法有直接测量法和间接测量法。直接测量法通常来说包括收集样本和叶面积测量两个步骤，虽然精确度高，但是具有很强的破坏性和工作量大的缺点，仅限于小冠层测量，对于树林、森林大冠层基本上是不可能的。间接测量法通过光线传播模型能相对简单地测量出光线的穿透情况，并基于lambert-beer定律对叶面积指数做出精确的估计。间接测量方法主要有斜点样方法、基于辐射测定法和孔隙度反演法，具有快速、简单等优势。其中斜点样方的问题是需要大量的专用探针且操作相当繁琐，因而不适应冠层较大的植物的参数获取；基于辐射测定法常常需要借助激光雷达(根据搭载平台分为地基、车载、机载、星载等类型)获取目标对象的空间信息，处理数据复杂，成本高昂。而孔隙度反演法的代表方法——半球图像法，通过鱼眼镜头获取接近180°冠层图像信息，经过图像处理获取冠层孔隙度，进而反演出叶面积指数。自even等将其用于植被的冠层的辐射监测以来，半球图像法以其超大视角范围和获取鱼眼图像设备价廉的优势得到了广泛的研究和应用。近年来，随着数码相机和数字图像处理技术的完善和进步，半球图像法的功能也日益完善，目前已成为各领域的重要研究方法。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种基于半球图像法测量叶面积指数的方法，在受时空约束条件下，快速、准确测量出叶面积指数。

为实现上述发明目的，本发明一种基于半球图像法测量叶面积指数的方法，其特征在于，包括以下步骤：

(1)、图像采集及预处理；

采用鱼眼摄像头拍摄植被冠层图像，然后对植被冠层图像先进行灰度化处理，再使用固定阈值法进行二值化处理，最后进行先膨胀后腐蚀的处理，从而得到轮廓清晰的二值化图像；

(2)、确定植被冠层图像的有效成像区域的几何中心o1(u0,v0)和圆形有效成像区域的半径r；

(2.1)、在二值化图像中，以左上角第一个像素点作为像素坐标系的原点o，建立像素坐标系xoy；

(2.2)、设置有效数据集，初始化为空；

(2.3)、从左往右逐行逐列扫描二值化图像，当扫描到某一像素点(x,y)的像素值为255时，进一步判断像素点(x+1,y)的像素值是否也为255，如果这两个像素点的像素值都为255，则将像素点(x,y)识别为轮廓边界点，并记录像素点(x,y)坐标保存至有效数据集，否则继续扫描，直到最后一个像素点扫描结束；

从右往左逐行逐列扫描二值化图像，当扫描到某一点(x,y)的像素值为255时，进一步判断像素点(x-1,y)的像素值是否也为255，如果这两个像素点的像素值都为255，则将像素点(x,y)识别为轮廓边界点，并记录像素点(x,y)坐标保存至有效数据集，否则继续扫描，直到最后一个像素点扫描结束；

(2.4)、在像素坐标系xoy中，构建圆方程x²+y²＝ax+by+c；

(2.5)、以x²+y²作为标签，x,y,1作为特征值，利用有效数据集训练基于机器学习的多元变量线性回归模型，得到参数a,b,c；

(2.6)、将参数a,b,c代入圆方程，得到圆心坐标半径为

(3)、对鱼眼镜头进行投影方程标定；

(3.1)、选取8×11的黑白棋盘板作为标准面，设置上、下四个相邻黑白棋格交线的交点为标定点，即黑白棋格图像的角点；在与标准面垂直的轴线上固定好镜头位置，使得透镜光轴垂直于标定板平面，拍摄标定板图像；

(3.2)、用激光笔标识好镜头中心在标准面上的投射点o2；

(3.3)、在世界坐标系中测定鱼眼镜头中心与标准面的距离l和黑白棋盘各标定点距离o2的距离h，计算出各标定点的入射角

(3.4)、使用harris角点检测法提取出各角点的像素坐标(u,v)，然后计算出各角点到鱼眼图像有效区域中心o1(u0,v0)的像素距离

(3.5)、调整标定板与鱼眼镜头中心的距离，在不同距离下重复(3.1)-(3.4)的操作，从而记录n组(r,θ)数据；

(3.6)、将提取的n组(r,θ)数据带入标定方程r＝fθ，再基于最小二乘法拟合，得到其中，θ＝(θ1,θ2...θn)^t,r＝(r1,r2...rn)^t，t表示转置；

(4)、采集待测点的植被冠层图像并处理；

(4.1)、采用鱼眼摄像头拍摄待测点的植被冠层图像，并进行灰度化处理；

(4.2)、对植被冠层的灰度图像进行二值化处理：

(4.2.1)、以图像圆心o1为起点，将灰度图像中的圆形有效成像区域按照拍摄的方位角和入射角分别进行分块，其中，灰度图像按拍摄的方位角一共分为m1个区域块，每个区域块又按入射角分为m2个区域块，共分为m1×m2个扇形块；

(4.2.2)、利用改进的otus算法获取每一个扇形块的二值化阈值t，并按照二值化阈值t进行二值化分割,得到待测点的植被冠层二值化图像；

(5)、提取多个角度的冠层孔隙度；

(5.1)、根据投影方程r＝fθ，在待测点的植被冠层二值化图像中按照(0～θ0)入射角范围取n个环，每个环对应宽度dθ＝θ°/n的入射角区间，由内到外，各环的半径为r1＝f·θ0/n,r2＝2f·θ0/n,r3＝3f·θ0/n......rn＝f·θ0，各环的入射角依次用θ1＝1/(2n)θ0,θ2＝3/(2n),θ0θ3＝5/(2n)θ0......θn＝(2n-1)/(2n)θ0表示；

(5.2)、计算各入射角下的冠层孔隙度t(θi)；

其中，t(θi)表示入射角θi下的冠层孔隙度，nleaf(θi)表示在入射角θi下背景像素点个数，nsky(θi)表示在入射角θi下叶片像素点个数；

(6)、利用反演公式计算叶面积指数；

叶面积指数的反演公式为：

其中，w′i为各个环的归一化权重因子，满足

本发明的发明目的是这样实现的：

本发明一种基于半球图像法测量叶面积指数的方法，先通过鱼眼摄像头拍摄植被冠层图像并进行预处理，再确定出植被冠层图像的有效成像区域的几何中心和圆形有效成像区域的半径，从而对鱼眼镜头进行投影方程标定；然后采集待测点的植被冠层图像并处理，再根据投影方程提取出多个角度的冠层孔隙度，进而反演出叶面积指数。

同时，本发明一种基于半球图像法测量叶面积指数的方法还具有以下有益效果：

(1)、一般的lai测量仪器拍摄条件最好满足万里无云，且避免太阳直射条件下进行；本发明的阈值分割算法具有自适应的优点，对光照条件的要求没有特别的限制，因此对测量环境要求较低，适合在各种天气条件下进行测量；

(2)、基于辐射测量的方法需要同时获取冠层上方和冠层底部的辐射量，当冠层较为高大时获取冠层上方的辐射不易实现，这种情况往往需要在附近空旷地测得的辐射量替代冠层上方的辐射量；且需要搭配激光雷达来获取冠层的三维空间信息，成本较高；与基于辐射测定的方法相比，本发明计算叶面积指数的方法只需要在冠层下方拍摄一次，对采集设备要求较低，只需要一个装配鱼眼镜头即可完成采集，具有价格低廉，且便于携带的优点；

(3)、本发明测量叶面积指数的方法使用冠层孔隙度来替代冠层间隙率，避免了基于辐射测定法测量冠层上方、下方辐射量的不便和对外部设备的依赖。只需通过数字图像处理技术，就可以很轻松地获取冠层孔隙度信息；

(4)、本发明采用的改进otus算法对图像分块进行二值化处理，减少因光照、鱼眼相机曝光度等因素对图像二值化品质的影响，提高了冠层孔隙度提取的精确度，进而提高叶面积指数计算的准确性。

附图说明

图1是本发明基于半球图像法测量叶面积指数的方法流程图；

图2是为了提取鱼眼图像轮廓信息的经过图像处理后的图像；

图3是使用圆方程拟合拟合鱼眼图像有效区域的效果图；

图4是建立投影方程图；

图5是某一次测量采集到的冠层图像；

图6是叶片对光线的反射率和透射比与波长的关系图；

图7是冠层图像分块示意图；

图8是各种阈值分割效果对比图；

图9是本发明算法测量结果与lai2200c测量结果的对比图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于半球图像法测量叶面积指数的方法流程图。

在本实施例中，如图1所示，本发明一种基于半球图像法测量叶面积指数的方法，包括以下步骤：

s1、图像采集及预处理；

采用鱼眼摄像头拍摄植被冠层图像，然后对植被冠层图像先进行灰度化处理，再使用固定阈值法(固定阈值设为30)进行二值化处理，最后进行先膨胀后腐蚀的处理，如图2所示，从而得到轮廓清晰的二值化图像；

s2、确定植被冠层图像的有效成像区域的几何中心o1(u0,v0)和圆形有效成像区域的半径r；

s2.1、在二值化图像中，以左上角第一个像素点作为像素坐标系的原点o，建立像素坐标系xoy；

s2.2、设置有效数据集，初始化为空；

s2.3、从左往右逐行逐列扫描二值化图像，当扫描到某一像素点(x,y)的像素值为255时，进一步判断像素点(x+1,y)的像素值是否也为255，如果这两个像素点的像素值都为255，则将像素点(x,y)识别为轮廓边界点，并记录像素点(x,y)坐标保存至有效数据集，否则继续扫描，直到最后一个像素点扫描结束；

从右往左逐行逐列扫描二值化图像，当扫描到某一点(x,y)的像素值为255时，进一步判断像素点(x-1,y)的像素值是否也为255，如果这两个像素点的像素值都为255，则将像素点(x,y)识别为轮廓边界点，并记录像素点(x,y)坐标保存至有效数据集，否则继续扫描，直到最后一个像素点扫描结束；

s2.4、在像素坐标系xoy中，构建圆方程x²+y²＝ax+by+c；在本实施例中，如图2所示，可以看出由于光线分布不均的影响，有的地方轮廓信息提取不完整，所以需要上述步骤进一步筛选有效的轮廓信息来进行圆方程拟合，最后通过圆方程拟合鱼眼图像有效区域的效果图，如图3所示；

s2.5、以x²+y²作为标签，x,y,1作为特征值，利用有效数据集训练基于机器学习的多元变量线性回归模型，得到参数a,b,c；

s2.6、将参数a,b,c代入圆方程，得到圆心坐标半径为

s3、对鱼眼镜头进行投影方程标定；

研究表明鱼眼镜头成像时遵循的模型可以近似为单位球面投影模型，可以将鱼眼相机的成像过程分解成两步：第一步，三维空间点线性地投影到一个球面上，它是一个虚拟的单位球面，它的球心与相机坐标系的原点重合；第二步，单位球面上的点投影到图像平面上。鱼眼镜头大多采用等距投影模型设计，这种投影方式最大的特点是来自入射角θ的一点m在成像平面上的像点(m')到图像中心o的距离r与入射角θ成线性关系,用方程表示为：r＝f·θ；那么，对鱼眼镜头进行投影方程标定的具体过程为：

s3.1、选取8×11的黑白棋盘板作为标准面，设置上下相邻四个黑白棋格交线的交点为标定点，即黑白棋格图像的角点；在与标准面垂直的轴线上固定好镜头位置，使得透镜光轴垂直于标定板平面，拍摄标定板图像；

s3.2、用激光笔标识好镜头中心在标准面上的投射点o2；

s3.3、在世界坐标系中测定鱼眼镜头中心与标准面的距离l和黑白棋盘各标定点距离o2的距离h，计算出各标定点的入射角

s3.4、使用harris角点检测法提取出各角点的像素坐标(u,v)，然后计算出各角点到鱼眼图像有效区域中心o1(u0,v0)的像素距离

s3.5、调整标定板与鱼眼镜头中心的距离，在不同距离下重复(3.1)-(3.4)的操作，从而记录n＝352组(r,θ)数据；

s3.6、如图4所示，将提取的352组(r,θ)数据带入标定方程r＝fθ，再基于最小二乘法拟合，得到其中，θ＝(θ1,θ2...θn)^t,r＝(r1,r2...rn)^t，t表示转置；

s4、采集待测点的植被冠层图像并处理；

s4.1、采用鱼眼摄像头拍摄待测点的植被冠层图像，如图5所示，共计4幅，并进行灰度化处理；在本实施例中，叶片对光线的反射率和透射比与波长的关系如图6所示，因为在小于490nm区域的光线受叶片反射和透射作用最小，这使得叶簇在天空背景下是黑色的，从而更容易区分天空和植物。蓝光波段在400-450nm，满足我们应用需求。为了更好地分割天空像元和树叶像元，因此选择rgb冠层图像的蓝色分量进行灰度化处理。

s4.2、对植被冠层的灰度图像进行二值化处理：

s4.2.1、以图像圆心o1为起点，将灰度图像中的圆形有效成像区域按照拍摄的方位角和入射角分别进行分块，其中，灰度图像按拍摄的方位角一共分为m1个区域块，每个区域块又按入射角分为m2个区域块，共分为m1×m2个扇形块；

在本实施例中，如图7所示，按拍摄的方位角一共分为8个区域块，每个区域块又按入射角分为10个区域块，共分为80个扇形块；这样分块的目的好处在于某一分块区域灰度分布直方图都存在明显的波峰波谷，从而可以更好地应用otus算法进行后续阈值分割；

s4.2.2、利用改进的otus算法获取每一个扇形块的二值化阈值t，然后将每一个扇形块中像素值大于等于阈值t的灰度值设置为255，并归为天空背景；将像素值小于阈值t的灰度值设置为0，归为叶片，得到待测点的植被冠层图像的二值化图像。

具体过程为：

1)、统计每个扇形块的灰度直方图，对每个扇形块分别计算阈值t＝0,2,…,255情况下的类间方差

其中，w1为天空像素占比，ni为某一扇形块内灰度值为i的像素点个数，sum为某一扇形块内总像素点个数；w2为叶面像素占比，w2＝1-w1；μ1为天空像素的平均灰度值，pi为某一扇形块内灰度值为i的像素占比，μ2为叶面像素的平均灰度值，

2)、在每个扇形块内，使得类间方差最大的阈值t作为该扇形块的二值化阈值t，然后使用该二值化阈值t进行二值化分割。

在本实施例中，图8给出了各种阈值分割效果对比图，其中，图(a)rgb原图；(b)固定阈值法分割，阈值为100；在圈出的两个地方，错误地将树叶像元识别为天空像元，阈值的选择对固定阈值分割影响较大，而阈值的选择往往只有通过人眼看到视觉效果人为调节，不具备实际应用的要求；(c)otus算法,因为无效区域大量黑色像素的影响，使得otus算法计算出的阈值往往偏低，往往容易将灰度较暗的天空像元识别为树叶像元；(d)基于高斯加权和的自适应阈值分割，邻域大小为81×81；该算法对于树干只提出了轮廓部分，很多细节丢失。这是因为高斯加权和的自适应阈值分割算法要求处理的图片，灰度不能太复杂，尽量符合混合高斯模型；(e)本文提出的阈值算法，在原有的otus算法上进行了改进，未将有效区域外的像素值参与计算，消除了无效区域大量黑色像素导致阈值偏低的影响，同时结合鱼眼镜头光线分布与入射角和方位角的关系，对图像进行分块处理，具有自适应性，阈值分割效果较好，满足本文对阈值分割的需求。

s5、提取多个角度的冠层孔隙度；

s5.1、根据投影方程r＝fθ，在待测点的植被冠层二值化图像中按照(0～θ0)入射角范围取5个环，每个环对应宽度dθ＝75°/5＝15°的入射角区间，由内到外各环的半径为r1＝f*15,r2＝f*30,r3＝f*45,r4＝f*60,r5＝f*75，各环的入射角依次用θ1＝7°,θ2＝23°,θ3＝38°,θ4＝53°,θ5＝68°表示。

s5.2、在假设叶片是黑色不透光的前提下，冠层孔隙度是光线穿透冠层不被捕获的概率，具体计算各入射角下的冠层孔隙度t(θi)为：

其中，t(θi)表示入射角θi下的冠层孔隙度，nleaf(θi)表示在入射角θi下背景像素点个数，nsky(θi)表示在入射角θi下叶片像素点个数；

s6、利用反演公式计算叶面积指数；

叶面积指数反演的理论公式为：

其中，p(θ)为入射角θ下的冠层间隙率，在利用半球图像法计算叶面积指数时通常使用冠层孔隙度t(θ)代替冠层间隙率p(θ)。

在实际计算中，通过利用5个环(θ1＝7°,θ2＝23°,θ3＝38°,θ4＝53°,θ5＝68°,dθ＝15°)的孔隙度信息对上面的公式离散化处理，得到实际的叶面积指数反演公式：

其中，w′i为各个环的归一化权重因子，满足未归一化前的权重值wi＝2sinθicosθiδθi，归一化的权重因子与归一化前权重因子两者之间的关系

如表1所示，各环的权重因子如下：

表1

在本实施例中，采集到学校四个样本地不同植被类型近50组冠层图像，并同时用lai-2200c测量每个实验点的叶面积指数。对采集到的冠层图像，采用本发明的算法进行叶面积指数的估算。然后对本发明的算法结果和lai-2200c的测量结果进行线性回归、相关性等分析，如图9所示。从全体数据分析结果(图9(a))来看，本文基于半球图像法lai的计算结果与lai-2200c的lai测量结果之间相关系数达到了0.90244。证明本论文的算法与lai-2200c的测量结果之间具有线性关系，说明本算法的精确度较高(约0.9)。

从更深层次上分析，两者的较高相关性从源头上是半球图像法采用的孔隙度变量与lai-2200c的间隙率变量之间较高的相关性。在将冠层看作泊松分布的假设下，两者从不同角度去量化树叶对光束的拦截概率。基于这点出发，可以用线性模型回归训练的参数对本文算法结果进行修正，即xr＝kx+b，其中xr为本文算法测量lai修正后的结果，k,b分别为图9(a)线性回归模型的斜率参数和截距参数。修正后的结果具有较高的精确度和准确度。

通过与lai2200c测量结果进行比较，二者相关系数达到0.9，肯定了本发明提供的基于半球图像法测量叶面积指数方法的可行性和可信性。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。