城市绿地地上生物量估测方法与流程

本发明涉及地上生物量测算领域，具体涉及一种城市绿地地上生物量估测方法。

背景技术：

城市绿地生物量提取，对于城市绿地资源监测和环境因子调查具有重要意义。同时，这些信息也可以用于掌握城市绿地内植物与环境的关系，城市绿地内植被生长发育和更新、演替的规律，对于城市绿地经营的可持续经营管理、生态系统碳循环研究及全球气候变化理解具有重要意义。常规的城市绿地地上生物量提取主要依赖于实测或以实测法为基础的统计分析方法，其精度往往不高，且难于在“面”上实用化推广。

目前，遥感影像估测地上生物量的方法主要是：1)从光学遥感影像提取与城市绿地相关的植被指数、波段之间的比值和光谱曲线特征作为自变量，将国家森林资源二类调查数据或样地实测数据作为因变量，建立生物量回归模型，从而估测城市绿地地上生物量，2)利用雷达、激光雷达、高光谱和高分影像，基于实测数据建立模型，进行反演估测城市绿地的地上生物量。以上两种方法都能大致反演城市绿地的地上生物量，但都存在问题。第一种方法，将多种植物混为一谈，仅仅建立一种模型对地上生物量进行反演，精度较低。第二种方法，虽然利用雷达数据增加了高度信息，但是雷达数据量庞大，操作复杂不易计算，增加了估测的数据量计算，且雷达的后向散射的强度随着生物量的增加而线性增加，达到一定生物量水平后，后向散射趋于饱和，会对估测精度造成一定的影响。同时，这两种方法，都是在像元尺度上对生物量进行反演计算，无法针对单株植物计算生物量，精度较低。

中国专利文献cn107449400a记录了一种森林地上生物量的测量系统及测量方法，通过无人机航拍影像，结合地面实测胸径数据建立模型，利用胸径-树冠面积回归方程和地上生物量-树冠面积回归方程获得林木的地上生物量。中国专利文献cn108921885a记录了一种综合三类数据联合反演森林地上生物量的方法，对高分辨率ccd影像和高光谱影像进行预处理，对激光雷达点云数据进行滤波，插值生成数字地形模型，点云数据归一化处理。然后分别基于预处理后的三个数据源提取纹理特征、光谱特征和点云结构特征，结合地面实测数据和提取的植被指数和纹理特征变量构建模型以预测森林地上生物量。中国专利文献cn104361592a记录了一种地上生物量估测装置，主要通过合成孔径雷达sar影像进行处理，基于散射数据估测树高，采用对象分割获得多边形林分，基于像元反射率建立多元线性模型估测地上生物量。然而以上发明方法中，都将目标区域中乔灌木按照同一种模型建模来推算地上生物量，并且没有将草地和林木区分开建模估测生物量，给地上生物量的测算带来误差。专利cn107449400中要求预测森林样地林木种类和建模样地相同，这在实际运用和操作中不太现实。而专利cn108921885a需要利用激光雷达测量树高，数据量极大，并且也将所有乔灌木混为一体，建立一个统一模型，这无形中也要求待测样地和建模样地林木种类一致。这些方法都没有以不同植物物种的角度建模，在识别物种的前提下进行地上生物量的估算。同时，城市绿地不同于城市森林，绿地中其他建设用地将会干扰地上生物量的测算，准确将植被区域提取后再计算显得尤为重要。

技术实现要素：

本发明所要解决的技术问题是提供一种城市绿地地上生物量估测方法，基于识别植物轮廓和物种前提，面向对象的城市绿地的地上生物量遥感反演方法，能快速和精确的计测大面积范围内的城市绿地地上生物量和单位面积生物量，并精确计算得到每一种植物、每一株植物的地上生物量，对于城市绿地资源监测和环境因子调查具有重要意义。

为解决上述技术问题，本发明所采用的技术方案是：

城市绿地地上生物量估测方法，包括以下步骤：

步骤一、获取目标区域的全色影像、高光谱影像数据，并对数据进行预处理；

步骤二、将高光谱图像进行分割，筛选出影像数据中的植被区域，构成目标区域植被影像数据，包括目标区域植被中的所有乔木、灌木和草地；

步骤三、利用全色影像，对植物区域的影像进行边缘检测，对植物物种的轮廓进行识别；

利用高光谱影像可见光、红边波段和近红外波段进行组合，定位、识别和聚集数据中最纯的像元，获取植物区域影像中纯净的端元波谱曲线，用于提取目标区域绿地中不同植物物种的端元波谱特征曲线并录入库；

步骤四、实测目标区域各种植物种类的地上生物量，建立目标区域绿地不同种类乔木和灌木的冠幅及其地上生物量的关系模型，草地建立面积和生物量的关系模型；

步骤五、根据步骤三中构建的高光谱影像下不同植物物种的端元波谱特征曲线，将影像中的所有数据进行分类，然后与全色影像植物物种的冠层轮廓监测结果进行叠加，识别定义每个轮廓的物种类别；

步骤六、根据步骤四中建立的各种植物对应的地上生物量的关系模型，结合步骤五中的识别结果，通过计算，获取目标区域绿地影像数据的地上生物量和每公顷地上生物量。

上述步骤一中数据的预处理包括：分别对全色影像和高光谱影像进行辐射定标、大气校正、几何校正、镶嵌和拼接。

上述的步骤二中图像分割的具体过程为：

利用最大类间方差算法otsu对归一化植被指数ndvi进行处理，将高光谱图像进行分割，即将影像数据中植被的区域和非植被的区域进行划分，获取目标区域影像数据，影像数据包括目标区域中的所有乔木、灌木和草地。

ndvi的公式为：

ndvi＝(ρnir-ρred)/(ρnir+ρred)；

式中，ρnir为近红波段的反射率，ρred为红波段的反射率。

otsu的公式为：

设ndvi影像数据的灰度级为l(g＝1,2,3,…,l)，处在灰度级i的像元个数用ni表示，ndvi的总像元个数n为：

直方图归一化后为：

式中，p(i)表示图像中灰度级为i的像元出现的概率。

用t表示阈值，将归一化图像中的像元按灰度级t划分为两类：c1≤t和c2＞t，两类出现的概率为：

假设两类像元的均值分别为m1和m2，ndvi图像全局像元均值为mg，那么：

w1*m1+w2*m2＝mg

w1+w2＝1

此时方差为：

δ²＝p1(m1-mg)²+p2(m2-mg)²

即：

δ²＝p1p2(m1-m2)²

式中，δ²为类间方差，当类间方差为最大值时，对应的阈值t为ndvi影像数据中植被和非植被区分的最佳分割阈值。

上述的步骤三的具体过程为：

基于索贝尔sobel算子进行边缘检测，然后进行植被轮廓识别，步骤为：

首先，设定sobel算子，公式为：

其次，利用全色高分辨率影像的像元集合z与a和b做平面卷积，分别实现经横向及纵向边缘检测的图像灰度值gx和gy的计算，公式如下：

最后，可计算全色影像的每一个像元的横向及纵向灰度大小，公式为：

使用近似值：

|f(gx，gy)|＝|gx|+|gy|

设置阈值为h，像元值f(gx，gy)大于h，即被标记为边缘点，这些边缘点的集合，即为轮廓识别的结果；

在全色影像轮廓识别的基础上，将每个轮廓内植物物种的高光谱影像作为数据集，通过高通滤波器模板对高光谱数据集g进行分离，得到噪声gn和信号gs，则有：

g＝gn+gs

通过最大化变换后数据的信噪比，即最大化信号协方差与噪声协方差的比值，得到变换矩阵v：

式中：yn和ys分别为变换后数据的噪声和信号，cov(·)表示协方差，cn和cs分别为数据的噪声和信号的协方差，cn＝cov(gn)，cs＝cov(gs)，上式优化问题可等价为：

式中，c表示数据的总体协方差，c＝cn+cs，根据拉格朗日乘子法，上式的最优解为：

cv＝λcnv

根据上式从大到小排列特征值，取前d个特征值对应的特征向量，可得到转换矩阵：

v＝[v1,v2,…,vd]

运用排除了噪声波段的mnf结果，通过ppi提取高光谱影像中的纯净像元，生成一副“像元纯度图像”，通过设定阈值，将纯净的像元提取，进行下一步计算；

将不同植物物种的端元波谱曲线提取出来，生成一条完整的端元波谱曲线并录入端元波谱数据库。

上述的步骤四的具体过程为：

调查记录城市绿地不同植物物种树冠投影面积lj、胸径dj和地上生物量yj，h、k、b和a为常数，建立不同植物物种的乔木和灌木yj、lj与dj的关系模型：

lj＝k*dj+b

yj＝a*dj^h

由上式可得出乔木和灌木的地上生物量为：

设草地的面积为rarea，c为常数，草地的地上生物量为yg:

yg＝c*rarea

上述的步骤五中数据分类采用波谱角分类方法，波谱角分类是在n维空间中把未知谱线与样本波谱曲线(步骤三中的端元波谱曲线)进行比较，选取最相似的曲线作为识别依据，从而进行归类的分类方法，即广义余弦夹角法。该方法将光谱曲线作为光谱空间的向量，通过计算两个向量之间的夹角，即计算待识别光谱与参考光谱之间的“角度”，进而确定二者的相似性，波谱角计算公式如下：

式中：α为未知波谱空间向量t与样本波谱空间向量r之间的夹角对比结果，ti和ri分别代表第i个波段上未知波谱t和样本波谱r的数值，参考波普库中的端元光谱，n为波段数。计算结果值为0～π/2，值越接近0，表示测试像元与参考光谱越接近。对于一个像元光谱，与所有的参考光谱都计算光谱角，其所属地物类别即为所有计算结果中α最小参考光谱所代表的地物类别。

上述的步骤六的具体过程为：

根据步骤五中的识别结果得到不同乔木和灌木植物物种的冠幅和物种类别，以及草地的面积和种类，结合步骤四中建立的不同物种的冠幅和生物量的关系模型，对目标区域的面积范围内的地上生物量进行反演，得到目标区域范围内的地上生物量(ubio)：

式中，ubio为目标区域内的生物量，yij为第j植物物种中第i株的生物量，yg为草地的生物量；

目标区域范围内的单位面积地上生物量(uavr)：

式中，uavr为目标区域内单位面积地上生物量，ubio为目标区域内的生物量，uarea为目标区域内总面积。

上述的生物量估算方法，还包括以下步骤：

步骤七、建立反演模型对目标区域绿地地上生物量进行反演，评价模型拟合的效果和精度。

上述的反演模型采取决定系数r²、均方根误差rmse和相对均方根误差rrmse评价模型，计算公式如下：

式中，xi为目标区域内地上生物量实测值；为目标区域内地上生物量实测平均值；为目标区域内地上生物量的模型估测值；n为样地的数量；i为目标区域样地编号。

本发明提供的一种城市森林地上生物量估测方法，基于高分辨率全色影像(dn值、纹理特征和空间几何特征)和高光谱影像(光谱特征)，用于识别植物轮廓和物种前提下，面向对象的城市绿地地上生物量遥感反演方法，以高光谱数据、高分辨率全色遥感影像和实测数据为数据源，构建了一套以植被冠层轮廓识别、高光谱端元波谱曲线特征提取和基于波谱角比对的物种识别为前提，结合不同物种冠层和地上生物量模型，进行城市绿地地上生物量反演的方法，本发明率先提出了面向对象的地上生物量遥感反演方法，能快速和精确的计测大面积范围内的城市绿地地上生物量和单位面积生物量，并精确计算得到每一种植物，每一株植物的地上生物量，计算结果效果理想、精度高，比起大多数基于像元的反演模型，更为快速、精确、全面。本专利提出的方法在估算大面积城市森林地上生物量有理想的效果。

附图说明

下面结合附图和实施例对本发明作进一步说明：

图1为本发明中预处理后的全色影像和高光谱影像；

图2为ndvi计算结果及otsu阈值分割图像；

图3为sobel算子的植物轮廓识别结果示意图；

图4为经过mnf和ppi后不同植被类型的端元波谱数据示意图；

图5为基于波谱角和轮廓检测的植物物种识别结果示意图；

图6为植物物种识别结果及基本信息。

具体实施方式

如图1-6中所示，城市绿地地上生物量估测方法，包括以下步骤：

步骤一、获取目标城市区域的高分辨率全色影像、高光谱影像数据，并对数据进行预处理；

步骤二、将高光谱图像进行分割，筛选出影像数据中的植被区域，构成目标城市区域影像数据，包括目标区域中的所有乔木、灌木和草地；

步骤三、利用高分辨率全色影像，对植物区域的影像进行边缘检测，对植物物种的冠层轮廓进行识别；

利用高光谱影像可见光、红边波段和近红外波段进行组合，定位、识别和聚集数据中最纯的像元，获取植物区域影像中纯净的端元波谱曲线，用于提取目标城市区域中不同植物物种的端元波谱特征曲线并录入库；

步骤四、实测目标城市区域各种植物物种的地上生物量，建立城市绿地区域不同乔木和灌木物种冠幅和地上生物量的关系模型，对草地建立面积和地上生物量的关系模型；

步骤五、根据步骤三中构建的高光谱影像下不同植物物种的端元波谱特征曲线，将影像中的所有数据进行分类，然后与全色影像植物物种的冠层轮廓监测结果进行叠加，识别定义每个轮廓的物种类别；

步骤六、根据步骤四中建立的各种植物对应的地上生物量关系模型，结合步骤五中的识别结果，通过计算，获取目标区域内的地上生物量和每公顷地上生物量。

上述步骤一中数据的预处理包括：分别对高分辨率全色影像和高光谱影像进行辐射定标、大气校正、几何校正、镶嵌和拼接。

上述的步骤二中图像分割的具体过程为：

利用最大类间方差算法otsu对归一化植被指数ndvi进行处理，将高光谱图像进行分割，即将影像数据中植被的区域和非植被的区域进行划分。

ndvi的公式为：

ndvi＝(ρnir-ρred)/(ρnir+ρred)；

式中，ρnir为近红波段的反射率，ρred为红波段的反射率。

otsu的公式为：

设ndvi影像数据的灰度级为l(g＝1,2,3,…,l)，处在灰度级i的像元个数用ni表示，ndvi的总像元个数n为：

直方图归一化后为：

式中，p(i)表示图像中灰度级为i的像元出现的概率。

用t表示阈值，将归一化图像中的像元按灰度级t划分为两类：c1≤t和c2＞t，两类出现的概率为：

假设两类像元的均值分别为m1和m2，ndvi图像全局像元均值为mg，那么：

w1*m1+w2*m2＝mg

w1+w2＝1

此时方差为：

δ²＝p1(m1-mg)²+p2(m2-mg)²

即：

δ²＝p1p2(m1-m2)²

式中，δ²为类间方差，当类间方差为最大值时，对应的阈值t为ndvi影像数据中植被和非植被区分的最佳分割阈值。

上述的步骤三的具体过程为：

基于索贝尔sobel算子进行边缘检测，然后进行植被轮廓识别，步骤为：

首先，设定sobel算子，公式为：

其次，利用全色高分辨率影像的像元集合z与a和b做平面卷积，分别实现经横向及纵向边缘检测的图像灰度值gx和gy的计算，公式如下：

最后，可计算高分辨率全色影像的每一个像元的横向及纵向灰度大小，公式为：

使用近似值：

|f(gx，gy)|＝|gx|+|gy|

设置阈值为h，像元值f(gx，gy)大于h，即被标记为边缘点，这些边缘点的集合，即为轮廓识别的结果；

在高分辨率全色影像轮廓识别的基础上，将每个轮廓内植物物种的高光谱影像作为数据集，通过高通滤波器模板对高光谱数据集g进行分离，得到噪声gn和信号gs，则有：

g＝gn+gs

通过最大化变换后数据的信噪比，即最大化信号协方差与噪声协方差的比值，得到变换矩阵v：

式中：yn和ys分别为变换后数据的噪声和信号，cov(·)表示协方差，cn和cs分别为数据的噪声和信号的协方差，cn＝cov(gn)，cs＝cov(gs)，上式优化问题可等价为：

式中，c表示数据的总体协方差，c＝cn+cs，根据拉格朗日乘子法，上式的最优解为：

cv＝λcnv

根据上式从大到小排列特征值，取前d个特征值对应的特征向量，可得到转换矩阵：

v＝[v1,v2,…,vd]

运用排除了噪声波段的mnf结果，通过ppi提取高光谱影像中的纯净像元，生成一副“像元纯度图像”，通过设定阈值，将纯净的像元提取，进行下一步计算；

将不同植物物种的端元波谱曲线提取出来，生成一条完整的端元波谱曲线并录入端元波谱数据库。

上述的步骤四的具体过程为：

调查记录目标区域城市绿地不同植物物种树冠投影面积lj、胸径dj和地上生物量yj，h、k、b和a为常数，建立不同植物物种的乔木和灌木yj、lj与dj的关系模型：

lj＝k*dj+b

yj＝a*dj^h

由上式可得出乔木和灌木的地上生物量为：

设草地的面积为rarea，c为常数，草地的地上生物量为yg:

yg＝c*rarea

上述的步骤五中数据分类采用波谱角分类方法，波谱角分类是在n维空间中把未知谱线与样本波谱曲线(步骤三中的端元波谱曲线)进行比较，选取最相似的曲线作为识别依据，从而进行归类的分类方法，即广义余弦夹角法。该方法将光谱曲线作为光谱空间的向量，通过计算两个向量之间的夹角，即计算待识别光谱与参考光谱之间的“角度”，进而确定二者的相似性，波谱角计算公式如下：

式中：α为未知波谱空间向量t与样本波谱空间向量r之间的夹角对比结果，ti和ri分别代表第i个波段上未知波谱t和样本波谱r的数值，参考波谱库中的端元光谱，n为波段数。计算结果值为0～π/2，值越接近0，表示测试像元与参考光谱越接近。对于一个像元光谱，与所有的参考光谱都计算光谱角，其所属地物类别即为所有计算结果中α最小参考光谱所代表的地物类别。

上述的步骤六的具体过程为：

根据步骤五中的识别结果得到不同乔木和灌木植物物种的冠幅和物种类别，以及草地的面积和种类，结合步骤四中建立的不同物种的冠幅和生物量的关系模型，对目标城市区域的面积范围内的地上生物量进行反演，得到目标区域范围内的地上生物量(ubio)：

式中，ubio为目标区域的生物量，yij为第j物种中第i株的生物量，yg为草地的生物量；

目标城市区域范围内的单位面积地上生物量(uaqr)：

式中，uavr为目标区域内单位面积地上生物量，ubio为目标区域内的生物量，uarea为目标区域内总面积。

上述的生物量估算方法，还包括以下步骤：

步骤七、建立反演模型对目标区域地上生物量进行反演，评价模型拟合的效果和精度。

上述的反演模型采取决定系数r²、均方根误差rmse和相对均方根误差rrmse评价模型，计算公式如下：

式中，xi为目标区域内地上生物量实测值；为目标区域内地上生物量实测平均值；为目标区域内地上生物量的模型估测值；n为样地的数量；i为目标区域样地编号。

实施例1：

1、本实施例试验区位于上海市奉贤区邬桥苗圃基地(121.41°e，30.96°n)，面积约3.83hm²，其植被类型丰富。其中植物包括：垂丝海棠(malushallianakoehne)、红花檵木(loropetalumchinense(r.br.)olivervar.rubrumyieh)、胡颓子(elaeagnuspungensthunb.)、鸡爪槭(acerpalmatumthunb.)、荚蒾(viburnumdilatatumthunb.)、金边黄杨(euonymusjaponicusthunb.var.aurea-marginatushort.)、金边马褂木(liriodendrontulipifera‘ureomarginatum’)、金叶皂荚(gleditsiatriacaanthos‘sunburst’)、木槿(hibiscussyriacuslinn.)、七叶树(aesculuschinensisbunge)、鹅掌楸(liriodendronchinense(hemsl.)sarg.)、香茶藨子(ribesodoratumwendland)、香樟(cinnamomumcamphora(linn)presl)、樱花(prunusserrulata)、玉兰(magnoliadenudatadesr.)、圆柏(sabinachinensis(linn.)ant.)、紫薇(lagerstroemiaindicalinn.)、榉树(zelkovaserrata(thunb.)makino)等。

2、借助无人机采集高分辨率影像和多光谱数据；

3、在研究区范围内设置28个方形样地(10×10m²)，样地中心点坐标使用gps(trimblegeoxh6000)测定，gps通过接收广域差分信号定位，精度优于0.5m。并在样地中记录树种并计数，同时测量每木的冠幅、胸径和树高，地上生物量通过方程结合胸径和树高进行计算，根据调查数据汇总，建立不同树种树冠和地上生物量的关系模型，建立草地和地上生物量的关系，如表1中所示：

表1

注：其中，树冠投影面积为lj，草地面积为rarea。

4、如图1中所示，数据预处理时，对高分辨率全色影像进行拼接，生成一幅完整的研究区高分辨率影像。同时结合地面实测控制点数据，采用二次多项式模型对影像进行几何精校正。借助传感器辐射定标数据对原始高光谱数据进行辐射定标，并利用经验线性模型结合地面实测标靶光谱数据进行大气校正。同时，利用地面实测控制点对高光谱影像进行几何精校正；

5、如图2中所示，利用经过预处理的高光谱影像，进行ndvi的计算，得到ndvi分布图，利用otsu法，求得阈值，本实例中为0.3，对ndvi进行分割，得到分布图中的植被区域；

6、利用经过预处理的高分辨率全色影像，根据步骤5中的植被区域进行裁剪，得到全色影像的植被分布区。

7、如图3所示，对全色影像的植被分布区，进行sobel边缘检测，识别不同植物物种的树冠轮廓以及草地轮廓；

8、如图4中所示，在识别植物物种的轮廓基础上，对实地调查的样地中已经标记的植物种类进行标定，对标定轮廓内的高光谱影像进行处理，包括最小噪声分离和纯净像元提取，提取不同典型物种的端元波谱曲线，并录入库；

9、如图5中所示，对其他为标定轮廓内的高光谱影像进行提取，将未知植物的波谱曲线与录入库的端元波谱曲线进行比对，利用波谱角比对的方法，对未知物种进行识别；

10、如图5和6中所示，计算植物物种的轮廓面积，结合识别结果，获得乔木和灌木的树冠面积，获得草地的分布面积，本例中共有不同植物物种23种，共计树冠面积2028个，草地125块。草地面积1019.19m²，乔木和灌木面积12735.29m²；

11、如图6中所示，带入步骤3中的23个关系模型，进行计算，获得影像区域的城市绿地地上生物量以及单位面积生物量，获得每一种植物，每一株植物的地上生物量；

12、采用采用决定系数(r²)、均方根误差(rmse)和相对均方根误差(rrmse)评价模型拟合的效果及估测精度，如表2中所示：

表2地上生物量各模型估算精度汇总表

本实施例中通过联合同时相获取的高分辨率全色影像和高光谱影像数据，然后对植被区域进行提取，识别植物物种轮廓，包括乔木和灌木的树冠轮廓和草地的轮廓，基于提取的特征变量建立城市绿地地上生物量估算模型，数据源分别从光谱和空间的角度记录了绿地结构信息，其数据之间互为补充，并且本方案是面向对象的地上生物量估测，脱离了以往基于像元的反演方法，且能识别不同植物物种，针对不同植物物种建立模型，测算地上生物量。故本方法增强了城市绿地地上生物量反演的能力和精度；结合表1和表2的数据验证结果表明，通过本实施例对城市绿地地上生物量进行提取，与使用其他相近遥感方法进行地上生物量估算结果相比，其相对均方根误差降低了5％以上。