范永祥 冯仲科 闫 飞 申朝永 关天蒴 苏珏颖
(1.季华实验室, 佛山 528000; 2.北京林业大学精准林业北京市重点实验室, 北京 100083;3.贵州省第三测绘院, 贵阳 550004)
可靠的森林资源信息是评估森林发展状况及制定森林发展计划的重要依据[1]。样地调查作为一种重要的森林资源信息清查方式,利用抽样技术在森林中选取适当数量、尺寸且具有代表性的样地,通过对样地每木检尺、归纳、总结便可反演森林发展现状及消长变化[2],是全国森林资源清查(简称一类调查)、森林资源规划调查(简称二类调查)及作业设计调查(简称三类调查)的重要手段[3]。近年来,随着遥感技术及计算机数据处理能力的发展,星载、机载等遥感手段被用于完成不同尺度森林调查及监测,但海量森林遥感数据的处理仍依赖于可靠、详细、完整的森林样地调查数据,以便对星载数据、机载数据及反演模型进行校准[4]。高效、精确地完成森林样地调查对掌握森林资源现状与发展动态至关重要。
森林样地调查通过每木检尺对样地中所有立木属性如树种、胸径、树高、第一枝下高、冠幅及立木位置等测定,经统计便可获取平均及总体的样地林分属性;而其他的属性(如叶面积指数、茎曲线、生物量及碳储量等)可利用树种、胸径、枝下高、冠幅及树高等建立反演模型进行估计[5-6]。每木检尺的精度直接决定了样地调查的精度及各反演模型的精度。传统样地调查中,使用一些简易的工具进行每木检尺等工作。其中,胸径采用胸径尺、轮尺等测量,树高采用布鲁莱斯测高器、超声波测高器等测定[6]。近年来,一些简单集成测距及测角传感器的新型设备被研发并用于森林样地清查中[7-11],这些设备虽然解决了每木检尺需要多人、多种设备配合测量的问题,但仍然是一类以人为主要观测者的设备,故仍然是一类费时、费力、主观性强的测量设备。但由于人工对测量结果的监督作用,胸径等测量值仍被认为是可靠的,经常作为评估其他调查方法的参考数据[5]。
随着遥感技术及计算机运算能力的发展,以面阵相机、ToF(Time of flight)相机及激光雷达(Light detection and ranging, LiDAR)等为主要观测传感器的近景遥感手段在森林样地调查中得到广泛应用[4,12-13]。该类型设备按采集方式可分为静态采集设备及动态采集设备。典型的静态采集设备如地面激光雷达(Terrestrial laser scanning, TLS)等,该类型设备需安置于样地中进行静态作业,根据安置站数可分为单次扫描及多次扫描两种作业模式,其中单次扫描模式容易受遮挡影响而无法获取样地完整的点云信息[14-15],多扫描模式则面临多站数据拼接难、作业效率低、不适合在较大样地中作业等问题[16-17]。动态采集设备将视觉传感器数据与IMU(Inertial measurement unit)及GNSS(Global navigation satellite system)等传感器数据融合,使设备可以在移动中进行数据采集并获取样地三维数据[18-20];特别是,SLAM(Simultaneous localization and mapping)技术的出现,使移动测量系统无需依赖于GNSS信号实现了高郁闭度林下定位及建图[21-23]。相比于传统设备,该类设备通过在样地进行数据采集后利用计算机对数据进行处理,实现每木检尺及样地属性提取,减少了调查人员的工作量且避免了主观观测误差[4]。然而,专业遥感设备仍面临价格昂贵、质量较大及无人工监督导致鲁棒性差等问题。
智能手机作为一种廉价、便携及高普及性的设备,是一种理想的森林样地调查设备。近年来,手机芯片、传感器等技术得到长足进步,使消费级面阵相机、ToF相机、激光雷达及IMU等传感器均可被嵌入智能手机设备上,Tango、AR Core、AR Kite等SLAM算法可在手机上运行并支持增强现实的实现,这些均为智能手机在森林样地调查中提供了新的机会[24-26]。文献[27]利用内嵌有ToF相机的智能手机扫描并获取了森林样地三维点云,获取立木胸径及位置均方根误差分别为1.15~1.91 cm及0.200~1.135 m。文献[28]利用内嵌有LiDAR的Apple iPad扫描并获取了样地三维点云,并提取到立木胸径估计值均方根误差为3.10~3.40 cm及0.109~0.218 m。文献[29]利用内嵌有LiDAR的iPhone及iPad开发了ForestScanner,用于实时估计立木胸径及位置,胸径及立木间相对距离估计值均方根误差为2.3 cm及3.1~4.4 m。课题组利用内嵌有ToF相机的智能手机研发了增强现实森林样地调查系统,该系统利用RGB-D SLAM系统实时数据实现胸径、立木位置及树高测量,利用基于树干回环检测提高了立木位置精度,利用增强现实场景实时人工监督测量结果并交互修正误差,胸径及立木位置的估计值均方根误差为1.00 cm及0.078~0.086 m[30-32]。然而,目前主流智能手机并非内嵌有ToF相机及LiDAR等传感器,这将影响智能手机在森林样地调查中的推广,故亟需研发基于单目SLAM系统(仅使用面阵相机及IMU作为数据源)的森林样地调查系统;此外,该系统在胸径估计中采用单帧深度图进行圆拟合估计,不仅鲁棒性差且不适合对倾斜立木进行精确胸径与立木位置估计。本文利用带有单目SLAM系统的手机构建森林样地调查系统;构建基于单目SLAM数据过滤胸高点云数据的方法;基于点到圆柱体表面距离及圆柱体切线到圆柱体表面距离构建胸径与立木位置精确估计算法;利用增强现实技术进行结果实时显示与操作交互。然后,该系统在5块32 m×32 m的方形样地中进行测试,并对立木位置及胸径估计结果进行评估。
本文单目SLAM系统以面阵相机为观测传感器、IMU为运动传感器,其中:面阵相机的优势在于长时间、长距离位姿估计误差较小,IMU的优势在于短时间、短距离内的位姿估计较准(即使在运动快速变化时),两种传感器数据融合可有效弥补各自的不足,是目前较为主流的SLAM系统之一。若设0~K时刻面阵相机对路标点{m}的观测序列为
Z0:K={Z0,Z1,…,ZK}
IMU运动估计序列为
u1:K={u1,u2,…,uK}
则SLAM过程将相机观测数据与IMU数据进行融合,并对各时刻位姿
x1:K={x1,x2,…,xK}
及路标点{m}进行实时估计(图1);从概率论角度而言,SLAM过程表达为后验概率推断
图1 SLAM过程贝叶斯网Fig.1 Bayesian network of SLAM process
(1)
SLAM解决方案正是通过对该后验概率推断采用不同的简化假设(如马尔可夫性假设等)并利用各种滤波及滑窗优化技术实现对位姿及路标点的后验推断[33-34]。
深度图的估计以SLAM过程获取的位姿及稀疏地图为输入,在朗伯面假设、光度一致性假设及可视性约束下,通过构建光度一致性非线性优化模型实现[35-36]。该非线性模型可简单表示为
(2)
式中IR(i,j)、Ik(i,j)——当前Patch及第k帧相应Patch中(i,j)偏移像素点光度
ck——第k帧颜色尺度
E——当前帧附近k帧相应Patch光度残差之和
通过对该光度残差平方和极小化,完成对当前Patch中心像素的深度等进行优化估计并通过区域生长获取深度图及其置信图。通过该方法获取立木树干深度图及由深度图转换的点云,如图2所示。由于森林中灌木等干扰、深度图估计中对树干圆柱体相切位置点估计不精确等原因,深度图中存在大量噪点,从而在利用该深度图进行立木位置及胸径等估计时造成影响,本文将构建相应的过滤算法及胸径等估计算法解决该问题。
图2 利用多视图影像重建立木树干深度图及单帧点云Fig.2 Depth images and single-frame point cloud of tree trunk reconstructed by multi-view stereo technology
本文测树系统基本工作框架及界面参考课题组已构建的RGB-D SLAM测树系统工作架构进行,该系统主要由单目SLAM系统、测树核心系统及测树用户交互系统3个模块组成,其中:SLAM系统用于位姿估计与稠密深度图估计,测树核心系统用于接收用户指令并利用SLAM输出数据完成每木检尺等工作,测树用户交互模块用于为用户提供增强现实窗口并通过手势等完成与测树核心系统交互。本测树系统增强现实交互界面如图3所示,每木检尺操作工作流程如图4所示。相比于RGB-D SLAM测树系统,新型测树系统算法利用多次、多视角测量数据对立木位置及胸径进行更精确估计,故在交互界面中增加“确认”按钮以协助该功能的完成。
图3 测树系统增强现实交互界面Fig.3 Augmented reality interface of tree measurement system
图4 森林样地调查中测树系统手势操作及计算流程图Fig.4 Gesture operation and calculation flowchart of tree measuring system in plot survey
在进行立木位置与胸径估计前,本文以被点击像素点坐标以及该时刻深度图为输入,先对深度图中噪点进行过滤,然后对深度图中立木圆柱体切点进行分类与精确估计,具体流程为:根据测量时设备与立木距离经验值,移除深度图中深度小于0.2 m或大于1.3 m的估计值;根据像素周围点深度的平滑度过滤树干周围灌木等遮挡噪点及树干圆柱体相切位置不合理点,平滑度计算公式为
(3)
其中,d为当前像素的深度,dij为x、y方向偏移量分别为i、j像素的深度,本文中S>0.008时深度估计被认为无效;查询距离被点击像素最近的、深度有效的像素点标记为种子点,查询该种子点周围8个像素中深度有效的像素并标记为新的种子点,然后继续索引新种子点周围有效深度像素及种子标记,直至无新种子被添加为止,最后将所有种子作为树干上的点;将深度图左右两侧边缘点作为切点候选点,利用相机位姿及内参将候选点变换到世界坐标系后,将各点按照极坐标角度排序,将中位数点作为提取到的树干圆柱体切线点(图5)。将深度图中各点及切点转换到局部坐标系后,该扫描预处理过程获取数据最终可表示为{T,{x},xL,xR,θ},其中:T表征该时刻相机的位置t及姿态R,{x}为位于树干圆柱体表面的估计点集合,xL、xR为树干左右两侧切点,θ为观测点对树干圆柱体的观测角,公式为
图5 树干圆柱体深度图噪声过滤结果Fig.5 Noise filtering result of tree trunk cylinder depth map
(4)
(5)
其中
式中θk——第k次扫描观测角
nk——第k次扫描中树干圆柱体总点数
x=tk+τnk,i
(6)
其中
式中τ——自变量
x——切线上任意一点
tk——第k次扫描时相机所在位置
nk,i——切线方向向量
由该切线构建的损失函数为
(7)
其中
式中jDk——立木胸径
利用所构建损失函数进行优化后,便可获取到树干圆柱体参数,然后将估计结果构建到增强现实场景中进行监督即可。
选择位于北京市近郊小西山的西山试验林场(39°58′N,116°11′E)为研究区域,该林场主要以人工林为主体。本文选择其中4块32 m×32 m的方形样地为研究对象。所选样地地面具有较少灌木且容易到达,不同径阶组比较均衡。表1总结了不同样地的基本属性。
表1 样地属性描述统计Tab.1 Summary statistics of plot attributes
为充分利用回环检测修正立木位置,本文数据采集中使用如图6所示测量轨迹(图中棕色圆圈为立木,数字表示观测顺序,绿色箭头表示行径方向,红色箭头表示闭环观测),即从样地中心出发并测量距离样地中心最近的立木;行经到西北角并对行径路线上邻近立木进行每木检尺;按仿航线法逐行完成每木检尺;回到样地中心,并对行径路线上立木进行二次观测。为充分研究新型设备胸径观测次数对测量精度的影响,本文试验中选择单次观测、正交观测、对称观测及环绕观测4种观测方式分别完成同一样地的每木检尺,4种观测方法对胸高圆柱体观测视角如图7所示(图中棕色圆圈为胸高圆柱体,红色方框为设备观测位置,黑色虚线为设备观测方向)。
图6 每木检尺中行径轨迹Fig.6 Path track in forest plot survey
图7 不同观测方式中观测视点位置与方向Fig.7 Viewpoint position and direction by different observation methods
为验证本文中新型设备在森林样地调查中的精度,使用胸径尺测量值作为胸径参考值、使用全站仪(Leica flexline TS06plus)对立木位置进行测量并作为参考值。使用偏差BIAS、均方根误差(Root mean squared error,RMSE)、相对偏差Brel及相对均方根Rrel对各测量值的偏差及变异性进行评估。
新型设备在5块样地中立木位置估计误差统计结果如表2所示,结果显示:X轴及Y轴方向的估计值偏差范围为-0.014~0.020 m,显然无偏;两轴方向估计值均方根误差小于0.09 m,显然变异性较小;立木位置最大观测误差小于0.25 m,且鲁棒性较好(图8、9);相比于单次观测,正交观测、对称观测及环绕观测方法中位置估计均方根误差、最大观测误差等略有减小,其中环绕观测均方根精度最高,但精度提升不显著。显然,从效率角度讲,单次观测方法最佳,且精度能满足日常需求;从精度及效率角度而言,正交观测及对称观测均可取得较高精度;此外,随着观测次数的增加,可一定程度上提高精度,但观测次数超过2次(如环绕观测)时,精度提升效果不显著。
表2 立木位置估计值精度
图8 立木位置误差分布Fig.8 Position errors of all trunks in plots
图9 立木位置估计误差箱状图Fig.9 Errors of tree position observations
使用不同观测方法获取5块样地中胸径估计误差统计结果,如表3所示。结果显示:不同测量结果与参考值一致且具有较小偏差(-0.85~-0.03 cm,图10),且均方根误差较小(1.32~2.51 cm);在不同径阶组下,胸径估计值的偏差、均方根表现较为一致(图11);随着观测次数的增加,测量值的偏差逐渐减小,显然胸高圆柱体数据的完善有助于对胸径进行更高精度的估计(图10、11);随着观测次数的增加,测量结果鲁棒性提高(图11),这是由于当树干无法近似为圆柱体时,单次测量值很难获取树干胸高圆柱体完整信息,不同侧面观测数据加入优化算法后有效抑制数据出现噪声。
表3 立木胸径估计值精度Tab.3 Accuracies of DBH estimations
图10 胸径估计值散点图Fig.10 Scatter plot of DBH estimation
图11 胸径估计值误差箱状图Fig.11 Errors of DBH observations under different DBH
(1)以内嵌有普通相机及IMU的智能手机为硬件系统、单目SLAM系统为数据源,研究立木胸径及位置的测量算法,以利用内嵌这些普通传感器的手机完成森林样地调查。首先基于平滑度设计了从多视图几何重构深度图中过滤胸高圆柱体数据的方法,并获取高鲁棒性胸高圆柱体表面点云数据及切线;然后,利用多组胸高圆柱体表面点云及切线数据设计并构建了胸高圆柱体拟合损失函数,实现对立木胸径及位置精确拟合的算法;最后,以该算法为基础在手机平台上构建了交互式增强现实测树系统,并用于森林样地调查。
(2)新型测树系统在5块32 m×32 m方形样地中进行了测试,每块样地设计并使用了4种不同的观测方法对立木胸高圆柱体观测。经与参考值进行对比表明,新型测树系统可高鲁棒性测量立木位置及胸径,且通过两次以上对胸高圆柱体观测可提高胸径估计精度(特别是不可近似为圆柱体的立木树干)。新型测树系统可完成森林样地调查。