梁晓军,庞 勇,陈博伟
(中国林业科学研究院资源信息研究所,国家林业和草原局林业遥感与信息技术重点实验室 北京 100091)
林业遥感真实性检验与地面调查中[1-3],需要高精度的单木位置信息。样地单木位置的精确观测已有大量研究及方法探索。传统林业样地单木调查中,一般首先确立样地中心或角点从而通过拉皮尺计算每株树与中心点或角点相对距离计算出单木相对位置,此种方法通常较为费时,且在计算坐标时由于皮尺的震荡及人为观测的随机性,单木位置的坐标值误差分布具有随机性。使用全球导航卫星系统(Global Navigation Satellite System,GNSS)定位时,有使用手持GNSS接收机直接定位单木位置,也有探索实时动态载波相位差分技术(Real-time Kinematic,RTK)定位的方法,但由于每个GNSS接收机单独定位误差会在2~5 m,导致单木位置不能得到精准确立[4-7]。使用高精度全站仪(Total Station,TS)观测的坐标值系统误差在3 mm之内,但通常观测为树皮某处,也不能精确描述单木位置。
地基激光雷达(Terrestrial Laser Scanning,TLS)高精度的参数反演能力[8-13],与机载激光雷达(Airborne Laser Scanning,ALS)及人工测量数据结合可以提供精确的地面验证数据[14-16]。基于TLS数据使用Hough变换圆检测的方法在拟合单木胸径(Diameter at Breast Height,DBH)时展示出较高效率和更高精度[17],同时基于提取的直径可以获得单木DBH处圆心,此圆心可以更好地描述冠幅覆盖下的单木位置。而基于仿射变换空间校正的思想在全局校正空间点时具有很高的指导意义,通常在点校正过程中,只需要选择好校正点,使用一般仿射变换空间校正方法就可以将全局点的坐标校正。在匹配相对坐标与地理坐标时,基于影像配准的方法也有大量研究,如王书涵等[18]提出一种基于TS和快鸟影像协同的山地定位匹配方法,Sarah等[19]提出一种基于高分辨率影像的高精度冠幅制图方法。
本研究基于TLS提取的DBH圆心点,经过外业调查(In-situ Measurement,IM)的DBH数据验证筛选,将单木IM坐标位置以空间校正方法匹配至DBH中心点处。同时,基于影像林隙特征,以ALS提取的数字冠层高度模型(Canopy Height Model,CHM)产品作为基准,将匹配的IM相对坐标单木位置再匹配至空间地理坐标系下,为大面积遥感精确验证提供可靠基础单木位置数据。
研究区位于黑龙江佳木斯孟家岗林场(130°32′~130°52′ E,46°20′~46°30′ N),平均海拔250 m,以低山丘陵为主。地势西南低,东北高,坡度在10°~20°之间。林场主要森林类型为人工林,主要树种长白落叶松(Larix olgensisHenry),占人工林面积的60.3%。
实验中激光雷达数据分别选用2017年6月采集的TLS数据及同期飞行的ALS数据。其中,TLS采用的是Trimble TX8三维激光扫描仪,此设备测距噪声小于2 mm,最大测量距离120 m,激光波长为1 500 nm,样地扫描为一个中心站加4个边缘站的设置,并均匀布设10个圆形靶球及7个黑白靶标;ALS数据采集于中国林业科学研究院资源信息研究所研制的CAF-LiCHy机载遥感系统[20],此系统搭载了Riegl公司的全波形激光雷达扫描仪LMS-Q680i,能够提供高精度高密度离散点云数据,孟家岗实验飞行绝对高度1 000 m左右,航速50 m·s-1,脉冲重复频率300 KHz,航带旁向重叠率约62%,研究样地范围点云数据密度10 pts·m-2以上。同时为研究需要,分别对TLS及ALS数据生产了0.25 m空间分辨率的CHM产品。
针对研究内容,选择了高、中、低3个不同株数密度人工落叶松林样地。不同株数密度样地TLS扫描场景图如图1所示。
TLS及ALS数据采集的同时对4 cm以上径阶单木做了每木检尺,包括同期实测的每木DBH、树高,并记录包括枯死及弯曲等的单木生长状态。原始单木位置沿用2012年皮尺及罗盘等测定IM结果。使用索佳免棱镜CX全站仪在样地中心均匀测定约1/3单木DBH处树皮位置并记录树号,同时测定TLS的7个黑白靶标点以将TLS坐标系平面坐标匹配至TS坐标系下,此TS的测距精度在200 m测距范围内可达到2.2 mm (2 mm固定误差加0.2 mm比例误差)。如图2给出了其中一个样地TS测点及靶标点分布图。
图1 不同株数密度样地地基激光雷达扫描场景图Fig. 1 The preview scenes of sample plots with different stem density
图2 全站仪测点及靶标点分布Fig. 2 The distribution of mesured points and targets points from total station
从图2可以看出,样地中靶标设置基本均匀分布在样地中心围绕TLS中心站。TS测点基本均匀分布在样地中的设置有助于降低后续空间位置校正中的全局单木匹配误差。
表1给出了3个不同株数密度样地的概况,表中3行自上而下分别对应高中低3种样地密度。
表1 样地概况Table 1 The profile of sample plots
研究的技术流程如图3。
图3 技术流程图Fig. 3 The technical flowchart
首先将TLS数据结合TS靶标测点坐标将TLS及TS数据匹配至同一坐标系下,并输出高密度点云数据;其次对TLS数据做预处理,然后对处理的点云分层栅格化、Hough变换圆检测得到提取的单木DBH及中心点;其后通过比较叠加位置TLS提取DBH与TS位置对应实测DBH,得出准确的TLS单木相对坐标位置。最后,基于TS位置匹配误差小的外围点将TS全部测点做空间校正匹配至TLS平面点上,然后将样地内全部IM单木位置基于TS精确位置匹配的外围点匹配至TLS坐标系下。另一部分,分别对TLS和ALS数据生产统一空间分辨率的CHM产品,目视分析同一区域的影像林隙特征,将匹配的全样地相对坐标单木位置匹配至空间地理坐标系下,以支持大范围林业遥感验证。
TLS原始数据使用Trimble Realworks软件基于10 cm圆形靶球将单测站数据拼站,后将拼站结果中黑白靶标扫描位置与TS测定的黑白靶标位置按顺序整理后做位置匹配,得到基于TS原始坐标系的点云数据。其后使用lastools[21]工具集对输出点云做去噪、地面点及非地面点分类、归一化、裁剪处理得到预处理结果。
单木DBH的提取主要依赖于高精度高密度的预处理点云数据,预处理点云拼站误差2 mm之内,点云密度可达3 000 pts·m-2。使用SEAL[22]软件胸径的提取是对预处理点云在DBH高度处上下均等取共14层5 cm厚点云分层栅格化,后基于影像Hough变换圆检测,提取直径与圆中心,从而得到扫描场景内每木DBH值与圆心点。
如图4所示,A、B分别为同一坐标系下的TLS提取DBH及TS测点位置对应实测DBH示意图,在R≤DBH且TLS圆心点与TS点在DBH差小于2 cm的条件下,即|D1-D2|≤2 cm时认为是同株树的点。从而可以通过空间校正的方法将此两点匹配起来,所以记录了树号的TS精确位置为将B的圆心(红点位置)匹配至A圆心后的结果,准确的TLS点位以TS点号对应。这样做的目的是确保作为匹配基准的TLS中心点位有对应准确单木树号。基于前面所述TLS点结果,可以对整个样地的IM单木坐标做基于仿射变换的TS与IM的空间点校正,从而有了基于TLS提取DBH中心点的整样地精确相对坐标单木位置。
式(1)为空间校正中仿射变换的变换公式,此变换基于平移、缩放、旋转、倾斜对原始点位做空间位置变换。
式中,Y为校正后坐标矩阵,X为校正前坐标矩阵,R为旋转矩阵,L为剪切矩阵,S为缩放矩阵,M为平移矩阵。
矩阵R、L、S、M由多组TS及TLS点(X、Y)拟合获得,所以当选取的校准点位置误差大时,最终匹配的Y累计误差会增大;为了控制整个样地位置匹配精度,需要尽量均匀选取样地外围对应点,从而控制整体点的变换。
TLS与TS及IM位置匹配完毕后,TLS的CHM产品与匹配的IM单木位置直接叠加即为空间校正好的相对坐标地面单木位置结果。
图4 全站仪坐标点匹配至地基激光雷达提取中心点示意图Fig. 4 The schematic diagram of points mesured from total station matching to the center points extracted from terrestrial laser scanning
使用lastools工具集分别对预处理的ALS及TLS数据生产0.25 m分辨率CHM产品,以支撑影像匹配。
在TLS坐标系下,匹配好的精确IM单木位置为相对坐标系坐标,而ALS获取的CHM产品具有精确的地理坐标信息。本研究使用CHM产品的林隙特征,将地面相对坐标IM单木位置和ALS的CHM地理位置匹配起来。
如图5所示,通过对图中所示对应区域的TLS(左)及ALS(右)的CHM林隙特征(图中黑色区域中CHM为0的地面)对比分析,以特征中尖峰、拐角等位置,结合匹配后的相对坐标点对可以分别在两个CHM产品上构造矢量线及矢量矩形,将矩形角点做为控制点,基于此特征及控制点,将地面相对坐标位置通过空间校正的方法匹配至ALS数据的地理坐标系下,完成地面坐标与地理坐标系下的ALS数据精确匹配。基于ALS的CHM局部最大值对匹配后的IM单木地理坐标做再平移检查,即假定匹配的IM坐标有待再移动至更佳位置,在样地外扩10 m后开始对匹配后的IM坐标作1 m、旋转1°窗口移动[23],如果某处整体树高相关性最高就认为此时结果为位置最优。
图5 CHM林隙特征Fig. 5 The characteristics of forest gap in CHM
精度评价包括两方面。一是对相对坐标地面数据的精度评价,首先使用实测DBH数据评价TLS提取的单木DBH值精度,即TLS提取单木位置对应DBH值与实测值的差;其次是DBH提取正确率,基于IM单木位置对比TLS提取位置是否有过检或漏检;然后是最终匹配的IM单木位置上DBH与匹配的TLS位置对应DBH差在2 cm之内即认为IM单木位置匹配正确;二是对匹配的相对坐标地面数据与ALS数据的地理位置的精度评价,即对匹配后的TLS与ALS的CHM产品目视检查,在CHM产品上检查像元差,位置差在两个像元内即认为匹配完成。
基于Hough变换圆检测方法提取的单木DBH值,与实测DBH值比较结果显示,在3个不同株数密度样地中,单木DBH处无其他目标严重遮挡、无倒木情况下,提取正确率为100%。图6给出了一种密度样地TLS提取的单木DBH值与实测DBH比较误差分布图,在误差定量分析中可以看出,提取DBH值与实测值差的绝对值在1 cm之内。
图6 地基激光雷达提取胸径值与实测胸径值差的绝对值分布Fig. 6 Distribution of DBH absolute value difference between in-situ measured and extracted from terrestrial laser scanning
从相关性分析中可以看出(图7),TLS提取的DBH与实测DBH拟合相关性达到0.99以上。两种分析结果可为IM单木位置匹配至TLS提取DBH中心处提供可靠基础。
图7 地基激光雷达提取胸径值与实测胸径值相关性结果Fig. 7 Relationship between DBH value extracted from terrestrial laser scanning and measured DBH value
TS单木点与TLS提取DBH中心点未完全匹配上时,如图8所示,给出了加实测DBH的TLS点(蓝色圆)及Hough变换前分层栅格化的TLS提取DBH处图像,叠加坐标未做空间校正的TS点加实测DBH(绿色圆)有位置偏差结果。TS及TLS点位置匹配完好时,如图8中为匹配后的TS(红色圆)加实测DBH点、TLS点匹配完好叠加结果,可以看出,匹配位置的DBH差别不明显,说明3.1中TLS提取的DBH值可靠。
图8 全站仪测点与地基激光雷达提取点匹配前后结果Fig. 8 The result of unmatched and matched with points mesured from terrestrial laser scanning cases of total station points position
基于上述匹配TS点结果,选择在样地外围的TS与IM对应点对,对IM数据做全样地单木位置空间校正,最终将样地内全部单木位置匹配至TLS坐标系下。如图9,显示了TLS与TS、IM以DBH显示匹配结果,图9.a给出了TS与TLS结果,图9.b给出了IM与TLS结果。TS与TLS匹配后TS位置与TLS圆心位置差均在1 cm之内;基于匹配TS位置校正的IM位置与TLS圆心差在1 m之内,小于外业测量冠幅南北或东西向值。本研究中比较了3种株数密度样地的匹配结果,显示在低密度林位置匹配精度最高,达到98.6%,中密度林达到92.5%,高密度林达到91%,说明样地密度对匹配精度有不同程度影响。
图9 (a) 地基激光雷达提取胸径点与全站仪测点匹配结果;(b) 地基激光雷达提取胸径点与外业实测坐标匹配结果Fig. 9 (a) Matching results of DBH points extracted from terrestrial laser scanning compared with points measured from total station; (b) Matching results of DBH points extracted from terrestrial laser scanning compared with points measured from in-situ
图10分别给出了低中高(a、b、c)3个株数密度样地的最终地面相对坐标位置数据匹配至ALS地理坐标系下的结果,图中展示了匹配后的TLS设站S1-S5、IM单木点、样地GPS角点及匹配ALS的CHM叠加结果,可以看出地理位置匹配数据中单木点与ALS的CHM很好地匹配起来。比较结果显示,地面数据位置与ALS数据位置基于影像林隙特征匹配精度0.5 m之内,即匹配后TLS的CHM与ALS的CHM结果显示误差在两个像元,且通过CHM树冠最大值检测匹配IM单木位置结果显示,基于TLS提取DBH中心点匹配的方法结合影像林隙特征匹配的空间数据结果最优。
图10 地理空间单木位置匹配结果Fig. 10 The matched results of geospatial position of individual trees
TLS在单木DBH提取上已经有大量应用研究,使用Hough变换等方法精确地提取匹配DBH中心点的方法,使单木位置可以用TLS提取DBH中心点描述,最近有学者使用即时定位与地图构建(Simultaneous Localization And Mapping,SLAM)算法及移动激光雷达估测DBH及单木位置,但估测误差大于1 cm[24-27],本研究基于此重点讨论一种基于TLS提取DBH中心点的单木空间位置匹配方法。方法基于TLS的Hough变换圆检测提取单木DBH及中心点,提取的DBH与实测值比较精度高,但在其他目标遮挡严重、有分叉木的情况下,提取DBH精度仍然会受限;在匹配TS点时,控制了TLS提取DBH与实测值差小于1 cm,即小于平均DBH的10%误差,从而邻近TLS提取DBH中心与TS点做为同一单木点做匹配;在匹配IM与TS单木点时,IM单木位置会有记录错误的情况,此时需要依据TLS提取DBH中心点与邻近IM单木点判断匹配结果,所以控制TLS提取DBH与实测DBH差放大至2 cm;以DBH差别及匹配位置差别小的样地单木外围点匹配TS及IM点也是为了控制最终匹配误差。
本研究讨论的点校正方法基于空间校正进行,匹配过程先对TS坐标做匹配,后基于此匹配结果对IM坐标再做匹配,且实验了TLS单独与这两种数据匹配的高精度结果。
本研究最终结果是将IM单木位置、TLS提取DBH中心通过匹配的TS坐标匹配起来,然后匹配至地理坐标系下,此过程的匹配结果会带有IM位置的随机误差;冯仲科等[28]研究表明使用TS野外林业定位信息精度均在mm级,所以如果使用TS将每木坐标按一定规则定位,再使用本方法匹配全样地单木位置则结果更优;此匹配方法可以发展为更为自动化的方法,以减少人工干预,提高效率。
基于地基激光雷达胸径提取的单木空间位置匹配方法,以地基激光雷达提取胸径圆心点为基准,将外业调查单木位置数据依据全站仪测量点位置使用空间校正方法匹配至地基激光雷达胸径圆心处,使外业调查单木位置误差降低至1 m之内,且在3种株数密度样地的实验结果显示,高、中、低密度样地中位置正确匹配精度分别达到98.6%、92.5%、91%,说明株数密度对位置匹配精度影响程度递减。基于匹配的外业调查单木位置数据,依据影像林隙特征将此数据与机载激光雷达的数字冠层高度模型产品做匹配,匹配位置误差在0.5 m之内,实现了高精度单木位置匹配结果。
此方法的发展为单木位置匹配校准提供了一种可靠方法;同时发展了基于影像特征的地理匹配方法,并为地面调查与地理空间数据匹配提供了思路;最后也为林业遥感大范围数据验证提供了可靠的数据基础,从而为遥感单木分割、参数反演提供坚实基础,进而指导森林高效经营。