树木年轮资料由于空间分布广泛、分辨率高, 能准确记录自然环境变化过程, 用于重建历史气候变化情况, 在气候学、林学、生态学等多个学科交叉研究中发挥着重要作用^[1]。21世纪之前, 树木年轮气候学的主要研究对象是乔木, 直到三齿蒿(Artemisia tridentata Nutt.)成为第一个被研究的灌木, 揭示了灌木年轮对气候变化的高敏感性, 有关灌木年轮的研究才日益受到关注^[2—3]。灌木对研究无乔木或少乔木分布区的气候、环境变化具有重要意义, 如今已成为中国寒旱区树木年轮学的主要研究对象^[4]。

目前, 利用灌木年轮信息进行气候、环境响应和重建的研究还相对较少, 这些研究通常以TRW作为指标, 个别研究使用了由TRW估算的BAI指标^[5], 而年轮密度^[6]和稳定同位素指标^[7]则少有涉及。研究中对气候因子涉及较多, 对环境因子如NDVI尚无涉及。TRW指标应用的前提假设是树木年轮各向生长速率基本相同, 如西伯利亚落叶松(Larix sibirica Ledeb.)的年轮就近似同心圆环, 当树木因环境影响偏向一侧生长时, 研究通过测量伸张木或硬压木垂直方向上的TRW来解决该问题^[8]。灌木不同于乔木, 严酷的自然环境常常导致灌木年轮扭曲变形或单侧停止生长, 这样TRW指标就难以满足灌木年轮学的研究需求^[9]。为此, 学者利用BAI指标代替TRW指标用于灌木年轮学研究, 但目前尚缺乏有效获取BAI的测量手段^[10]。常用的专业年轮分析系统如WinDENDRO和LINTAB等只能测量TRW和年轮密度, 这严重限制了灌木年轮学的发展^[11]。

随着电子扫描、图像处理等先进技术的不断发展, 目前关于年轮图像提取、年轮参数准确测量的研究报道越来越多。数字图像处理技术中的边缘检测、图像分割等方法已被成功应用于树木年轮信息提取。研究发现, 边缘检测中的Canny算子对年轮提取效果好于其它算子的提取效果^[12—13]。U-Net卷积神经网络能够实现年轮图像的准确自动识别^[14]。在年轮参数测量方面, AutoCAD和GIS已被成功用于年轮参数测量^[15—16]。以上研究极大的丰富了年轮测量的方法和手段。研究中如果能将年轮图像分类结果和专业测绘软件如GIS有机结合, 相信定能打破现有年轮测量技术瓶颈, 推动灌木年轮学发展。

鉴于此, 本研究以多枝柽柳为研究对象, 拟采用U-net深度学习的方法, 实现多枝柽柳年轮早材区域的自动分割, 利用GIS为图像赋予坐标的方法, 完成多枝柽柳年轮各项参数的测量。突破专业年轮测量软件的功能限制。利用BAI与早材起始生长线长度的比值计算TRGR指标, 解决TRW指标对灌木径向生长代表不充分的问题。研究结果有助于解决年轮指标去趋势难的问题。结果可推动灌木年轮学的发展。

1 材料与方法 1.1 多枝柽柳TRGR的测量方法 1.1.1 多枝柽柳取样与年轮盘预处理

样品采集时间为2021年7月, 采样地在新疆石河子市147团, 玛纳斯河灌区莫索湾总干渠西侧, 古尔班通古特沙漠南缘, 地理位置(44°42′18″N, 86°14′3″E)。采集时随机选择12棵多枝柽柳植株, 并从其灌丛中选择一个分支, 沿地表向上截取约100cm长的茎杆。在实验室用切割机将多枝柽柳茎杆切成约2cm厚的轮盘, 选择一个切面, 依次使用400目、600目、1200目的砂纸进行打磨直至切面光滑平整。使用爱普生(EPSON)V39扫描仪, 根据实际轮盘大小设置正方形扫描区域的边长, 选择2400点每英寸的扫描分辨率, 在黑暗条件下对处理好的轮盘进行扫描, 最后共得到430张轮盘扫描图像。需注意切割后的多枝柽柳轮盘样品容易开裂, 切割后应尽快处理并扫描, 若不能及时处理时, 可将切割样品湿润后放入密封袋内进行短期保存。

1.1.2 U-net早材语义分割

选择宁霄等^[14]设计构建的基于U-Net卷积神经网络的年轮盘图像分割网络, 对多枝柽柳轮盘扫描图进行图像分割。使用LabelMe对260幅多枝柽柳轮盘扫描图的早材区域进行标注用作为训练样本集。研究为减少训练的时间和成本, 训练过程中使用512×512像素的轮盘图片进行训练。利用训练好的模型对多枝柽柳轮盘扫描图进行图像分割, 提取年轮盘早材区域。

1.1.3 GIS测量TRGR

基于孙甜甜等^[17]对梭梭年轮的测量方法, 研究对其方法进行了简化(图 1)。具体是取消了在GIS中制作规定大小的正方形图框, 改为直接在扫描时确定扫描范围大小, 例如扫描一张分辨率为2400点每英寸, 边长为3.5 cm的正方形区域, 研究根据该图大小在GIS中创建相应倍数的渔网, 渔网参数设置为：左下为0, 右上为350000(为扫描区域的整数倍)。渔网投影设置为WGS_1984_World_Mercator。将扫描后的多枝柽柳轮盘图像通过GIS图像配准工具展布在建立的渔网上, 输出配准后的轮盘图像每个像素大小是105.836105 m, 而3.5 cm扫描区域每边上各有3307个像素点, 每个像素点大小是105.836105×10^-7 m, 则轮盘实际大小和配准后图像比例为1 ∶ 10000000, 上述步骤说明, 配准只是按比例放大了多枝柽柳轮盘图像; 在ENVI中对U-net分割结果图进行重采样, 使其恢复到轮盘扫描图的原始分辨率。对重采样的分割结果在GIS中同样完成图像配准并赋予坐标。在ENVI中对具有坐标的早材分割结果图进行图像提取, 波段阈值设置为大于等于1, 使用栅格转矢量工具, 将提取出来的早材区域保存为矢量格式。在GIS中将年轮早材矢量图转换为线状图。使用GIS编辑工具, 选择各年早材生长起始线, 删除早材停止生长线, 手动完善个别不连续的早材起始线。使用平滑工具对提取的早材起始生长线进行平滑, 经过实验发现, 平滑值用原图分辨率/512(分割结果图分辨率)×100000(图像配准的放大倍数)×0.01较为适合。计算各年早材起始生长线长度。将早材起始线状图转换成面状图, 计算各年年轮BAI。利用公式1计算多枝柽柳各年TRGR。

(1)

式中, BAI_n+1表示第n+1年年轮的轮环面积, 单位为mm²; L_n表示第n+1轮圆环内周长, 单位为mm; TRGR_n+1表示第n+1年年轮生长速率, 单位为mm/a。

1.2 NDVI时间序列分析

研究从美国国家航空航天局NASA的EOS/MODIS数据网(https://ladsweb.modaps.eosdis.nasa.gov)下载2000年至2020年空间分辨率250 m, 时间分辨率8天, 轨道号为h23v04、h24v04的MOD09Q1数据共1914幅。该数据2000和2001年有部分缺失, 分别为2000年第1—49期和2001年第169、177期数据。利用MRT(MODIS Reprojection Tools)软件进行图像拼接和投影转换(椭球体WGS-84, 投影Mercator)。利用MOD09A1数据地表反射率250 m的状态信息, 获取每幅影像晴空像元的NDVI^[18]。

通过分析2000年缺失的数据属于冬季, 采样点尚被冰雪覆盖, 此时, NDVI几乎为零, 因此, 研究采用零值填补缺失的各期数据。2001年夏季缺失的2期数据则采用前后两期图像值的平均值来替代。为使2020年的数据能够被拟合, 研究利用2020年的数据充当了2021年的数据。处理后的数据经Timesat Savitzky-Golay滤波时间序列拟合, 获取2000—2020年点样尺度NDVI季节曲线^[19]。该曲线共有13个参数, 研究选取季节峰值、季节振幅、季节初期的增长速率、大季节积分、小季节积分、季节开始的值和季节结束的值等7个关键参数(图 2)分别与多枝柽柳TRGR进行相关性分析, 以检验TRGR测量结果。

2 结果与分析 2.1 多枝柽柳TRGR测量结果

基于U-net深度学习后训练的模型, 研究选择一组未参与训练数据集的多枝柽柳轮盘扫描图像进行年轮早材区域分割(图 3)。将分割结果按图 1流程处理, 测量出各年年轮相关参数, 通过公式1计算TRGR, 结果如图 4所示, 同一茎杆上相邻位置样品测量的TRGR基本一致, 表明测量结果稳定, 符合同一茎杆上一定距离内相同年份TRGR相一致的理论预期。

U-net图像语义分割基本能够实现多枝柽柳年轮早材区域的准确提取, 但对一些年轮特征极不明显、侵蚀或有侧枝发生的部位, 早材提取结果仍有部分不完全的现象, 此时需根据多枝柽柳轮盘扫描配准图, 对早材起始生长线进行一定程度的修补。有关使用GIS测量年轮参数准确性的问题, 孙甜甜等^[17]已利用WinDENDRO专业年轮测量软件进行过详细验证, 研究对此未作展开。

2.2 NDVI时序分析结果对多枝柽柳TRGR测量结果的检验

NDVI作为叶绿素含量的表征, 能直接反映植物的长势, 理论上是较温度、降水、帕尔默干旱严重指数等^[20]众多环境因子与树木径向生长关系最为直接和密切的因子^[21]。研究发现NDVI与树木径向生长存在明显正相关关系, 研究据此利用NDVI和年轮指标的相关性, 来衡量年轮指标的优劣^[22—23]。

研究利用S-G滤波对2000—2020年点样尺度NDVI时间序列进行拟合, 获取了拟合曲线各季节关键参数(表 1)。研究将测量的TRGR与NDVI季节曲线关键参数间进行相关性分析, 发现10个受测样本中有8个轮盘的TRGR与各参数间存在正相关关系(表 2), 最高相关系数达到0.98(P < 0.01)。针对20年生多枝柽柳TRGR与NDVI季节曲线关键参数间表现出无相关性的问题, 分析发现这与2008年植株茎杆一侧形成层死亡, 其后死亡一侧周皮逐年风化脱落, 在2008和2009年部分木质部缺失关系密切(图 5)。

为此, 研究除去2008—2009年的干扰, 将该植株TRGR与NDVI时间序列季节参数分段进行相关性分析, 结果两者间均具有较高相关性(表 2)。2000—2007年TRGR与NDVI时间序列季节参数间的整体相关性要高于2010—2020年(图 6), 分析原因与2015年石河子市修建景观河, 玛纳斯河向石河子垦区蘑菇湖水库的输水减少, 向玛纳斯河下游输水增多有关^[24]。莫索湾干渠径流增加为采样点的多枝柽柳提供了额外的地下水分补给, 抵消了2015年高温对采样点多枝柽柳生长的抑制^[25], 在一定程度上导致多枝柽柳生长好于采样点其他植物的生长。

综上所述, 测量的多枝柽柳TRGR变化与采样点历年环境变化和植被长势间存在高度一致性, 测量结果与事实相一致, 结果符合理论预期。

2.3 TRGR指标优势分析 2.3.1 TRGR与TRW测量结果的比较分析

树木年轮学研究中, 研究者通常利用乔木TRW与NDVI的相关性重建历史NDVI, 并用于气候学研究。但是在无乔木或少乔木分布区, 建群种植物往往是小乔木或灌木, 这些植物年轮的生长往往偏心程度大、不规则, 生长锥打孔抽样对于各环重心偏移较大的树干所测年轮宽度准确性不高^[3]。研究以多枝柽柳轮盘为例, 模拟生长锥打孔, 测量了36个方向各向TRW与季节峰值的相关性, 结果如图 7所示, 各方向测量的TRW序列与季节峰值的相关系数存在很大差异。初步证明传统测量乔木年轮宽度的方法不适用于灌木年轮学研究。

针对多枝柽柳偏心生长现象, 为进一步分析TRGR对于TRW指标的优势, 研究构建了四种树木径向生长模型, 同心圆、重心线性、重心非线性和重心紊乱变化模型(图 8)。模拟生长锥取样, 分析各方向TRW与标准轮宽(同心圆模型轮宽)的相关性, 结果如图 9所示, 随着模型中各轮重心变化从有序到无序, 四种模型各方向的轮宽与标准轮宽相关系数平均值依次降低(图 10), 表明对于各轮重心位置变化较大的灌木, 通过生长锥取样准确测量年轮宽度的难度较大。但本研究测量的TRGR不会受到重心变化的影响, 因此, TRGR指标更适合于灌木年轮学的研究。

2.3.2 TRGR与BAI测量结果的比较分析

研究将测量的多枝柽柳BAI指标与采样点NDVI时间序列季节参数间进行相关性分析, 结果如表 3所示, 其中, 年轮数少的轮盘BAI序列与NDVI季节参数间呈负相关, 年轮数多的轮盘呈正相关。2016年以来采样点NDVI峰值呈逐年降低趋势(表 1), 1—4a多枝柽柳BAI因其内在递增规律, 未能反映出该期环境的变化特征。7—10a多枝柽柳年轮BAI序列呈递增变化趋势恰能与2000以来采样点NDVI总体呈现增加的变化趋势相一致。理论上, 第n-1年的BAI会影响第n年维管形成层的长度, 进而影响第n年的积材, 故BAI指标存在累积效应。研究为进一步阐明BAI指标的不足, 构建了一个已知年轮宽度的模型(图 11)来说明该问题, TRGR与已知同心圆TRW测量结果基本一致, 两者均能用于树木年际间生长状况的比较, 但BAI指标因存在递增规律性, 不利于年际间树木生长状况比较。因此, TRGR较BAI指标优势明显。

3 讨论

基于前人研究成果, 本研究将U-net深度学习网络与专业测绘软件GIS相结合, 实现了对多枝柽柳年轮扫描图像的自动化识别和测量, 研究构建的方法可以替代专业年轮分析系统的大部分功能。应用中除了可以使用U-net神经网络进行图像分类外, 还可以使用几乎所有图像监督及非监督分类方法, 图像分类结果均可基于GIS软件, 完成各种年轮参数测量。但该方法也存在不足：研究基于完整的树盘展开, 取样会对植株茎杆造成不可逆的损伤; 树盘的打磨、扫描、标注等工作繁重, 所需时间长; 涉及GIS、ENVI、Python等多个软件的联合应用。X射线计算机断层扫描(CT扫描)作为一种无损无侵入的技术, 近年来也被广泛用于年轮宽度测量, 在树木年代学中发挥着重要作用^[26]。树木CT二维扫描图像也可以使用数字图像处理技术完成图像分割^[27]。今后如果能将CT扫描与本研究测量方法相结合, 通过系统集成, 研发野外便携式年轮无损测量设备, 将会有效克服研究方法的不足, 推动树木年轮测量技术革新。

TRW、BAI等现有年轮指标中既包含环境变化信息, 也包含树木生长发育及物理结构等非环境变化信息。在树木年轮学研究中如何减少非环境趋势的干扰, 突显指标中的环境趋势, 一直是该学科致力解决的关键问题^[28—29]。在前人总结研究的基础上^[30], 本文将其归纳为三种趋势, 即胸径趋势、年龄趋势和环境趋势(图 12)。其中, TRW指标有随胸径增加而减小的趋势, BAI指标有随胸径增加而增大的趋势, 两者都具有树木结构性趋势。TRGR指标计算的是单位长度形成层的年活动强弱, 该指标不受胸径趋势的影响, 树木年轮指标的年龄趋势属于低频信号趋势, 树木寿命越长, 其对短时间指标序列的影响就越小。研究对比发现, TRGR指标在应用中更具优势, 原因是该指标既能准确反映年轮整体径向生长信息又能消除胸径趋势, 从而更加清晰的表征环境变化趋势。TRGR指标在应用中只需考虑年龄趋势的去除, 为解决年轮指标去趋势难的问题明确了工作重点。研究中由于生长速率指标本身已去除了胸径趋势, 低频年龄趋势对短时间年轮序列测量结果影响不大的原因, 单株年轮生长速率测量结果与采样点NDVI关键拟合参数间均具有高度相关性(表 1)。研究分析中TRGR指标与NDVI关系密切, 该结果为多枝柽柳年轮信息重构NDVI或反之提供了有利证据。

年龄趋势是一种生长趋势, 树木通常在幼龄和老龄时长的慢, 在青年时长的快^[31]。因此, 今后研究中可基于龄级生长速率分析法, 通过Gompertz或Logistic生长方程拟合各年龄滑动平均生长速率^[32], 获取各年龄生长速率基准值, 利用实测生长速率与龄级生长速率基准值之比消除年龄趋势, 将有望解决TRGR指标去除年龄趋势的问题, 据此所得的TRGR年表将会更具环境解释力。

4 结论

研究采用的TRGR指标能够反映多枝柽柳年轮整体径向生长情况和用于计算积材, 可以取代TRW指标作为寒旱区灌木年轮学研究的新指标。TRGR指标反应单位维管形成层的年活动强弱, 相比于BAI指标更适于树木年际间生长状况的比较。构建的TRGR指标理论上可有效提高交叉定年准确性, 解决年轮指标中环境信息不易有效分离的问题。研究测量的多枝柽柳TRGR与点样尺度NDVI时间序列拟合参数间均具有较高正相关性, 结果为利用灌木年轮信息重构NDVI或反之提供了有利证据。研究构建的年轮测量方法对树木年轮信息无损测量具有重要参考价值。研究结果可推动寒旱区灌木年轮学发展。