CN112577906A - 城市绿地土壤含水率检测方法 - Google Patents

Abstract

城市绿地土壤含水率检测方法，对高光谱数据和全色影像进行预处理；设置样地，并对样品进行含水率测定；通过影像分割和轮廓识别得到植被区域和非植被区域，选取非植被区域的裸土样本获取裸土的端元波谱曲线；提取裸露土壤分布区域；利用热红外数据针对植物覆盖区和裸露土壤区分别建立反演模型并结合实测值，计算反演模型中的常量值，反演区域土壤含水率。本发明结合图像分割、轮廓识别和波谱角分类，区分了植被覆盖区域和裸露土壤区域，建立了两种土壤含水率模型，能够间接和直接的检测城市绿地的土壤含水率，采用的方法便捷实用，可实现大范围面积的土壤含水率检测，比起大多数的土壤含水率测定方法，更为快速、精确、全面。

Description

城市绿地土壤含水率检测方法

技术领域

本发明涉及土壤含水率测算领域，具体涉及一种城市绿地土壤含水率检测方法。

背景技术

土壤含水率是一个重要的土壤物理参数，它是生态系统水循环、植物生长、土壤承载能力等研究领域中必不可少的基本参数。首先，土壤含水率影响着地表与大气的水分和能量交换。另外，水是植被生长发育的基本条件，土壤含水率是植被健康生长的保障，也是灌溉管理和城市绿地管理预报中的重要参数。因此，土壤含水率是评价土地资源的重要指标，是精准林业中极为关键的参数，土壤含水率的监测有着十分重要的意义。土壤含水率的监测方法和手段也成为人们一直关注的热点问题之一。

目前，土壤含水率的方法主要是：1)土壤含水率直接监测方法主要有：取土烘干称质量法、中子仪法、张力计法、频域反射仪法和时域反射仪法等。2)土壤含水率间接测定方法主要是基于土壤水分反射光谱法，该方法在测量时依据土壤反射率光谱曲线与土壤含水量的变化规律，通过提取相应的特征波段，基于数学统计方法，进行拟合分析而确定出土壤含水量与土壤反射率相关性模型，并对土壤含水量进行定量反演。以上两种方法都能测定土壤含水率，但都存在问题。第一种方法，能比较准确的获得土体剖面含水率，可测定土层较多，但一般都依赖于密集的野外采样，并且成本偏高、调查周期较长、容易受制于采样时间和采样范围。另外，这些方法以点测量为基础，代表性差，体现区域土壤含水率空间变异性尤为困难，难以实现大尺度土壤含水率变化的实时监测和快速更新；第二种方法，通过建立遥感影像特征波段与土壤含水率之间的关系，能反演土壤的含水率。但是，现阶段的土壤含水率遥感反演，一种是基于多光谱卫星或者航拍影像进行估算反演，一种是利用地物光谱仪测定土壤含水率并建模反演。但是，多光谱影像在区分土壤与其他地物信息上有难度，常常产生混合像元，降低反演精度。而多数是利用地物光谱仪来进行测定，测量的范围有限，很难在大面积范围内快速的计测土壤含水率。

中国专利文献CN 104697937A记录了一种土壤属性高光谱识别技术方法，利用遥感数据，获取不同时间的土壤高光谱图像，预处理后，经过监督分类得到裸露土壤，提取所述裸露土壤的地表反射率，依据所述裸露土壤的地表反射率建立裸露土壤地表反射率反演模型，室内测定土壤可蚀性数据并计算土壤K值；根据所述土壤K值和所述地表反射率反演模型中的光谱数据，建立影响可蚀性K的土壤属性的高光谱模型，解决了高光谱遥感技术不能用于测定土壤可蚀性的问题。中国专利文献CN 106918566A记录了基于高光谱数据的土壤水分、盐分信息联合提取方法，对数据进行预处理后根据反射率求归一化反射率、表观吸收率一阶导数和二阶导数。建立基于主成分分析的评价系统，利用敏感波段，使用逐步回归方法，建立模型，对土壤水分进行反演计算。以上两种方式，都只能测算裸露土壤的含水率，对于被植被遮挡的区域，无法直接测算出植被遮挡下土壤的含水率，且没有良好的方法将高光谱影像中的植被区域和非植被区域区分开，对土壤含水率的测算带来了误差。

发明内容

本发明所要解决的技术问题是提供一种城市绿地土壤含水率检测方法，用于识别城市绿地植被和裸露土壤的前提下，面向城市绿地的土壤含水率快速测定方法。能快速和精确的计测大面积范围内的城市绿地中土壤的含水率，对于城市绿地土壤物理性质监测和指导灌溉具有重要意义。

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

城市绿地土壤含水率检测方法，包括以下步骤：

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

步骤二、在地面设置样地，样地包括植被覆盖区域和裸土区域，在植被覆盖区域中采集植被冠层叶片及其所生长处的表层土壤混合样品，在裸土区域采集表层土壤混合样品，并对样品进行土壤含水率测定；

步骤三、将高光谱图像进行分割，区分出影像数据中植被覆盖的区域和非植被区域，非植被区域包括裸露土壤区域和其他区域，植被覆盖区域包括城市绿地中的所有乔木、灌木和草地；

步骤四、利用全色影像，对目标区域的影像进行边缘检测，对植被区域和非植被区域的轮廓进行识别，结合高光谱图像分割结果和全色影像的识别结果，作为最终植被区域和非植被区域的影像数据集；

利用高光谱影像不同波段进行组合，结合地面数据，选取典型裸露土壤样本区，通过设定阈值，定位、识别和聚集数据中纯净的像元，获取裸露土壤数据中纯净的端元波谱曲线，用于非植被区域影像数据集中提取裸露土壤分布区域；

步骤五、利用步骤五中构建的高光谱影像下裸露土壤的端元波谱特征曲线，将非植被区域影像中数据进行分类，提取裸露土壤分布范围；

步骤六、利用热红外数据，将水分胁迫指数CWSI、冠层相对温度CRTD和地表相对温度SRTD作为土壤含水率反演模型的参量；

步骤七、对植被覆盖区和裸露土壤区域的土壤含水率分别建立反演模型，结合样地实测值，计算出模型中常量值，反演区域土壤含水率。

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

上述的步骤三中图像分割的具体过程为利用归一化植被指数NDVI和大津Otsu算法将植被区域和非植被区域进行初次划分：

利用归一化植被指数NDVI和大津0tsu算法对图像进行分割。

首先，计算NDVI，公式为：

NDVI＝(ρ_NIR-ρ_RED)/(ρ_NIR+ρ_RED)

式中，ρ_NIR为近红波段的反射率，ρ_RED为红波段的反射率；

其次，利用大津Otsu算法对图像进行分割，步骤如下：

设NDVI图像的灰度级为L(G＝1，2，3，…，L)，处在灰度级i的像元个数用n_i表示，NDVI的总像元个数为：

直方图归一化后为：

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

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

假设两类像元的均值分别为m₁和m₂，NDVI图像全局像元均值为m_g，那么：

w₁*m₁+w₂*m₂＝m_g

w₁+w₂＝1

此时方差为：

δ²＝p₁(m₁-m_g)²+p₂(m₂-m_g)²

即：

δ²＝p₁p₂(m₁-m₂)²

式中，δ²为类间方差，类间方差为最大值时，对应的阈值T为NDVI图像的最佳分割阈值。

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

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

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

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

通过下式，可计算全色高分辨率影像的每一个像元的横向及纵向灰度大小：

使用近似值：

|f(G_x，G_y)|＝|G_x|+|G_y|

设置阈值为h，像元值f(G_x，G_y)大于h，即被标记为边缘点，这些边缘点的集合，即为轮廓识别的结果；然后在全色影像轮廓识别的基础上，参照样地数据，将非植被区域中裸露土壤部分的高光谱影像作为数据集，通过高通滤波器模板对高光谱数据集G进行分离，得到噪声G_n和信号G_s，则有：

G＝G_n+G_S

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

式中：Y_N和Y_S分别为变换后数据的噪声和信号，Cov(·)表示协方差，C_N和C_S分别为数据的噪声和信号的协方差，C_N＝Cov(GN)，C_S＝Cov(G_S)，上式优化问题可等价为：

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

CV＝λC_NV

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

V＝[V₁，V₂，…，V_d]

运用排除了噪声波段的MNF结果，通过PPI提取高光谱影像中的裸露土壤的纯净像元，通过设定阈值，将裸露土壤的纯净像元进行提取，并绘制纯净像元的波谱曲线。

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

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

上述的步骤六中，反演的具体参量公式为：

式中，CWSI为水分胁迫指数，T_l为植被覆盖区域的植被冠层温度，T_wet和T_dry分别代表在相同的气象环境中冠层温度的下限和上限，温度数据直接由热红外数据获得；

式中，CRTD为冠层相对温度，SRTD为地表相对温差，T_vmax为植被覆盖区的植被冠层温度最大值，T_vmin为植被覆盖区的植被冠层温度最小值，T_smax为地表裸露土壤区域温度最大值，T_smin为地表裸露土壤区域温度最小值，温度数据直接由热红外数据获得。

上述的步骤七中反演模型分为：

Y_c＝αCWSI+βCRTD+a

Y_s＝ySRTD+b

式中：Y_c为植被覆盖区域土壤含水率，Y_c为裸露土壤区域土壤含水率，CWSI为水分胁迫指数，CRTD为冠层相对温度，SRTD为地表相对温差，α、β、γ、a和b为常数；

利用室外采集到的样品，在室内测定冠层叶片含水率及其对应生长处的土壤含水率，测定裸露土壤含水率，带入上述公式，求解出α、β、γ、a和b，作为整个待测目标区域反演的参数。

上述的城市绿地土壤含水率检测方法，还包括：

步骤八、采取决定系数R²和相对均方根误差rRMSE评价模型拟合的效果和精度，计算公式如下：

式中，xi为目标区域土壤含水率实测值；

为目标区域土壤含水率平均值；

为目标区域土壤含水率的模型估测值；n为样地的数量；i为目标区域单一样地。

本发明提供的一种城市绿地土壤含水率检测方法，基于高分辨率全色影像、高光谱和热红外数据，利用光谱特性计算植被指数，结合轮廓识别和波谱角分类，区分了植被覆盖区域和裸露土壤区域，结合热红外数据，计算相应指数，建立了两种土壤含水率模型，采用的方法便捷实用，可实现大范围面积的土壤含水率检测，比起大多数的土壤含水率测定方法，更为快速、精确、全面，本专利提出的方法在估算大面积城市绿地土壤含水率有理想的效果。

附图说明

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

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

图2为NDVI计算结果及Otsu阈值分割图像得到植被区域示意图；

图3为经过MNF和PPI后裸露土壤的端元波谱数据示意图；

图4为经过波谱角识别的土壤分布特征示意图；

图5为土壤率反演结果示意图。

具体实施方式

如图1-5中所示，城市绿地土壤含水率检测方法，包括以下步骤：

步骤一、借助卫星、有人机或无人机采集全色影像、高光谱影像；

步骤二、在地面设置样地，并在样地中采集冠层叶片和表层(0～20cm)土壤混合样品，用密闭容器带回实验室进行含水率的测定；

步骤三、对高分辨率全色影像、高光谱影像和热红外影像数据预处理为：借助传感器辐射定标数据对原始高光谱据进行辐射定标，并利用经验线性模型结合地面实测标靶光谱数据进行大气校正；然后，利用地面实测控制点，采用二次多项式模型对高光谱影像进行几何矫正，然后进行镶嵌和拼接；

步骤四、利用经过预处理的高光谱影像，进行NDVI的计算，得到NDVI分布图，利用Otsu法，求得阈值为0.3，对NDVI进行分割，得到分布图中的植被区域；步骤五、利用全色影像，对目标区域的影像进行边缘检测，对植被区域和非植被区域的轮廓进行识别，结合高光谱图像分割结果和全色影像的识别结果，作为最终植被区域和非植被区域的影像数据集；

在高光谱影像上，对应实地调查的样地中裸露土壤区域进行标定，对标定区域的高光谱影像进行处理，包括MNF噪声分离和PPI纯净像元提取，提取裸露土壤的端元波谱曲线，并录入库；

步骤六、利用波谱角分类的方法，设定阈值为0.1，对影像中裸露土壤部分区域进行识别；

步骤七、利用三个关于水分和温度的公式(CWSI、CRTD和SRTD)，建立模型，并带入实测数据和热红外数据，计算出相应的常数，确定模型。

步骤八、利用模型，分别对植被覆盖区域和裸露土壤区域进行土壤含水率的反演检测。

步骤九、利用决定系数R²和相对均方根误差rRMSE评价模型拟合的效果和精度。

上述的步骤四中图像分割的具体过程为利用归一化植被指数NDVI和大津0tsu算法将植被区域和非植被区域进行初次划分：

利用归一化植被指数NDVI和大津Otsu算法对图像进行分割。

首先，计算NDVI，公式为：

NDVI＝(p_NIR-ρ_RED)/(ρ_NIR+ρ_RED)

式中，ρ_NIR为近红波段的反射率，ρ_RED为红波段的反射率；

其次，利用大津Otsu算法对图像进行分割，步骤如下：

设NDVI图像的灰度级为L(G＝1，2，3，…，L)，处在灰度级i的像元个数用n_i表示，NDVI的总像元个数为：

直方图归一化后为：

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

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

假设两类像元的均值分别为m₁和m₂，NDVI图像全局像元均值为m_g，那么：

w₁*m₁+w₂*m₂＝m_g

w₁+w₂＝1

此时方差为：

δ²＝p₁(m₁-mg)²+p₂(m₂-m_g)²

即：

δ²＝p₁p₂(m₁-m₂)²

式中，δ²为类间方差，类间方差为最大值时，对应的阈值T为NDVI图像的最佳分割阈值。

上述的步骤五中识别的具体过程为：

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

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

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

通过下式，可计算全色高分辨率影像的每一个像元的横向及纵向灰度大小：

使用近似值：

|f(G_x，G_y)|＝|G_x|+|G_y|

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

然后在全色影像轮廓识别的基础上，参照样地数据，将非植被区域中裸露土壤部分的高光谱影像作为数据集，，通过高通滤波器模板对高光谱数据集G进行分离，得到噪声G_n和信号G_s，则有：

G＝G_n+G_S

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

式中：Y_N和Y_S分别为变换后数据的噪声和信号，Cov(·)表示协方差，C_N和C_S分别为数据的噪声和信号的协方差，C_N＝Cov(G_N)，C_S＝Cov(G_S)，上式优化问题可等价为：

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

CV＝λC_NV

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

V＝[V₁，V₂，…，V_d]

运用排除了噪声波段的MNF结果，通过PPI提取高光谱影像中的裸露土壤的纯净像元，通过设定阈值，将裸露土壤的纯净像元进行提取，并绘制纯净像元的波谱曲线。

上述的步骤六中数据分类以裸露土壤纯净像元的波谱曲线为参照，采用波谱角分类方法识别非植被区域内的裸露土壤区域，波谱角分类是在n维空间中把未知谱线与步骤四中的样本波谱曲线作为样本进行比较，选取最相似的曲线作为识别依据，从而进行归类的分类方法，即广义余弦夹角法，该方法将光谱曲线作为光谱空间的向量，通过计算两个向量之间的夹角，即计算待识别光谱与参考光谱之间的“角度”，进而确定二者的相似性，波谱角计算公式如下：

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

上述的步骤七中，反演的具体参量公式为：

式中，CWSI为水分胁迫指数，T_l为植被覆盖区的植被冠层温度，T_wet和T_dry分别代表在相同的气象环境中冠层温度的下限和上限，温度数据直接由热红外数据获得；

式中，CRTD为冠层相对温度，SRTD为地表相对温差，T_vmax为植被覆盖区的植被冠层温度最大值，T_vmin为植被覆盖区的植被冠层温度最小值，T_smax为地表裸露土壤区域温度最大值，T_smin为地表裸露土壤区域温度最小值，温度数据直接由热红外数据获得。

上述的步骤九中反演模型分为：

Y_c＝αCWSI+βCRTD+a

Y_s＝γSRTD+b

式中：Y_c为植被覆盖区域土壤含水率，Y_c为裸露土壤区域土壤含水率，CWSI为水分胁迫指数，CRTD为冠层相对温度，SRTD为地表相对温差，α、β、γ、a和b为常数；

利用室外采集到的样品，在室内测定冠层叶片含水率及其对应生长处的土壤含水率，测定裸露土壤含水率，带入上述公式，求解出α、β、γ、a和b，作为整个待测目标区域反演的参数。

上述步骤九具体过程为：

采取决定系数R²和相对均方根误差rRMSE评价模型拟合的效果和精度，计算公式如下：

式中，xi为目标区域土壤含水率实测值；

为目标区域土壤含水率平均值；

为目标区域土壤含水率的模型估测值；n为样地的数量；i为目标区域单一样地。

实施例1：

1、本实施例试验区位于上海市奉贤区邬桥苗圃基地(121.41°E，30.96°N)，面积约3.83hm²。其植被类型丰富。其中植物包括：植物包括：垂丝海棠(Malus halliana Koehne)、红花橙木(Loropetalum chinense(R.Br.)Oliver var.rubrum Yieh)、胡颓子(Elaeagnuspungens Thunb.)、鸡爪槭(Acer palmatum Thunb.)、荚蒾(Viburnum dilatatum Thunb.)、金边黄杨(Euonymus japonicus Thunb.var.aurea-marginatus Hort.)、金边马褂木(Liriodendron tulipifera‘ureomarginatum’)、金叶皂荚(Gleditsia triacaanthos‘sunburst’)、木槿(Hibiscus syriacus Linn.)、七叶树(Aesculus chinensis Bunge)、鹅掌楸(Liriodendron chinense(Hemsl.)Sarg.)、香茶藨子(Ribes odoratum Wendland)、香樟(Cinnamomum camphora(Linn)Presl)、樱花(Prunus serrulata)、玉兰(Magnoliadenudata Desr)、圆柏(Sabina chinensis(Linn.)Ant.)、紫薇(Lagerstroemia indicaLinn.)、榉树(Zelkova serrata(Thunb.)Makino)等；

2、借助卫星、有人机或无人机采集高分辨率全色影像、高光谱影像和热红外影像数据；

3、在研究区范围内设置28个方形样地(10×10m²)。样地中心点坐标使用GPS(Trimble GeoXH6000)测定，GPS通过接收广域差分信号定位，精度优于0.5m。并在样地中采集冠层叶片和土壤表层(0～20cm)的土样进行混合，带入实验室内，利用烘干称质量法对叶片和土壤的含水率进行检测；

4、如图1中所示，数据预处理时，对高分辨率全色影像、高光谱影像和热红外影像数据预处理为：借助传感器辐射定标数据对原始高光谱据进行辐射定标，并利用经验线性模型结合地面实测标靶光谱数据进行大气校正；然后，利用地面实测控制点，采用二次多项式模型对高光谱影像进行几何矫正，然后进行镶嵌和拼接；

5、利用经过预处理的高光谱影像，进行NDVI的计算，得到NDVI分布图，利用Otsu法，求得阈值，本实例中为0.3，对NDVI进行分割，初次识别目标区域中的植被区域和非植被区域；

6、如图2中所示，利用Soble算子对植被区域和非植被区域进行轮廓识别，结合NDVI的识别结果，确定目标区域中的植被区域和非植被区域的分布特征；

7、如图3中，在高光谱影像上，对应实地调查的样地中裸露土壤区域进行标定，对标定区域的高光谱影像进行处理，包括MNF噪声分离和PPI纯净像元提取，提取裸露土壤的端元波谱曲线，并录入库；

8、如图4中所示，利用波谱角分类的方法，设定阈值，本实例中为0.1，对影像中裸露土壤部分区域进行识别；

9、利用三个关于水分和温度的公式(CWSI、CRTD和SRTD)，建立模型，并带入实测数据和热红外数据，计算出相应的常数，确定模型；

Y_c＝-0.0474CWSI-0.0472CRTD+0.0980

Y_s＝-0.0616SRTD+0.1016

10、如图5中所示，利用模型，分别对植被覆盖区域和裸露土壤区域进行土壤含水率的反演检测；

11、采用采用决定系数(R²)和相对均方根误差(rRMSE)评价模型拟合的效果及估测精度。分别为植被区域(0.75和12.52％)和裸露土壤区域(0.84和11.12％)，效果理想；

12、本实施例中通过联合同时相获取的高分辨率全色影像、高光谱影像数据和热红外数据，然后基于提取植被覆盖区和裸露土壤区域。分别建立反演模型，用于检测大范围面积内城市绿地的土壤含水率，检测效果理想，方便快捷。

Claims (8)

1.城市绿地土壤含水率检测方法，其特征在于，包括以下步骤：

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

步骤二、在地面设置样地，样地包括植被覆盖区域和裸土区域，在植被覆盖区域中采集植被冠层叶片及其所生长处的表层土壤混合样品，在裸土区域采集表层土壤混合样品，并对样品进行土壤含水率测定；

步骤三、将高光谱图像进行分割，区分出影像数据中植被覆盖的区域和非植被区域，非植被区域包括裸露土壤区域和其他区域，植被覆盖区域包括城市绿地中的所有乔木、灌木和草地；

步骤四、利用全色影像，对目标区域的影像进行边缘检测，对植被区域和非植被区域的轮廓进行识别，结合高光谱图像分割结果和全色影像的识别结果，作为最终植被区域和非植被区域的影像数据集；

利用高光谱影像不同波段进行组合，结合地面数据，选取典型裸露土壤样本区，通过设定阈值，定位、识别和聚集数据中纯净的像元，获取裸露土壤数据中纯净的端元波谱曲线，用于非植被区域影像数据集中提取裸露土壤分布区域；

步骤五、利用步骤五中构建的高光谱影像下裸露土壤的端元波谱特征曲线，将非植被区域影像中数据进行分类，提取裸露土壤分布范围；

步骤六、利用热红外数据，将水分胁迫指数CWSI、冠层相对温度CRTD和地表相对温度SRTD作为土壤含水率反演模型的参量；

步骤七、对植被覆盖区和裸露土壤区域的土壤含水率分别建立反演模型，结合样地实测值，计算出模型中常量值，反演区域土壤含水率。

2.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述步骤一中数据的预处理包括：分别对全色影像、高光谱影像和热红外影像数据进行辐射定标、大气校正、几何校正、镶嵌和拼接。

3.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述的步骤三中图像分割的具体过程为利用归一化植被指数NDVI和大津0tsu算法将植被区域和非植被区域进行初次划分：

利用归一化植被指数NDVI和大津Otsu算法对图像进行分割；

首先，计算NDVI，公式为：

NDVI＝(ρ_NIR-ρ_RED)/(ρ_NIR+ρ_RED)

式中，ρ_NIR为近红波段的反射率，ρ_RED为红波段的反射率；

其次，利用大津Otsu算法对图像进行分割，步骤如下：

设NDVI图像的灰度级为L(G＝1，2，3，…，L)，处在灰度级i的像元个数用n_i表示，NDVI的总像元个数为：

直方图归一化后为：

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

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

假设两类像元的均值分别为m₁和m₂，NDVI图像全局像元均值为m_g，那么：

w₁*m₁+w₂*m₂＝m_g

w₁+w₂＝1

此时方差为：

δ²＝p₁(m₁-m_g)²+p₂(m₂-m_g)²

即：

δ²＝p₁p₂(m₁-m₂)²

式中，δ²为类间方差，类间方差为最大值时，对应的阈值T为NDVI图像的最佳分割阈值。

4.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述的步骤四的具体过程为：

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

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

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

通过下式，可计算全色高分辨率影像的每一个像元的横向及纵向灰度大小：

使用近似值：

|f(G_x，G_y)|＝|G_x|+|G_y|

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

然后在全色影像轮廓识别的基础上，参照样地数据，将非植被区域中裸露土壤部分的高光谱影像作为数据集，通过高通滤波器模板对高光谱数据集G进行分离，得到噪声G_n和信号G_s，则有：

G＝G_n+G_S

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

式中：Y_N和Y_S分别为变换后数据的噪声和信号，Cov(·)表示协方差，C_N和C_S分别为数据的噪声和信号的协方差，C_N＝Cov(G_N)，C_S＝Cov(G_S)，上式优化问题可等价为：

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

CV＝λC_NV

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

V＝[V₁，V₂，…，V_d]

运用排除了噪声波段的MNF结果，通过PPI提取高光谱影像中的裸露土壤的纯净像元，通过设定阈值，将裸露土壤的纯净像元进行提取，并绘制纯净像元的波谱曲线。

5.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述的步骤五中数据分类以裸露土壤纯净像元的波谱曲线为参照，采用波谱角分类方法识别非植被区域内的裸露土壤区域。

6.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述的步骤六中，反演的具体参量公式为：

式中，CWSI为水分胁迫指数，T_l为植被覆盖区的植被冠层温度，T_wet和T_dry分别代表在相同的气象环境中冠层温度的下限和上限，温度数据直接由热红外数据获得；

式中，CRTD为冠层相对温度，SRTD为地表相对温差，T_vmax为植被覆盖区的植被冠层温度最大值，T_vmin为植被覆盖区的植被冠层温度最小值，T_smax为地表裸露土壤区域温度最大值，T_smin为地表裸露土壤区域温度最小值，温度数据直接由热红外数据获得。

7.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，所述的步骤七中反演模型分为：

Y_c＝αCWSI+βCRTD+a

Y_s＝γSRTD+b

式中：Y_c为植被覆盖区域土壤含水率，Y_c为裸露土壤区域土壤含水率，CWSI为水分胁迫指数，CRTD为冠层相对温度，SRTD为地表相对温差，α、β、γ、a和b为常数；

利用室外采集到的样品，在室内测定冠层叶片含水率及其对应生长处的土壤含水率，测定裸露土壤含水率，带入上述公式，求解出α、β、γ、a和b，作为整个待测目标区域反演的参数。

8.根据权利要求1所述的城市绿地土壤含水率检测方法，其特征在于，还包括：步骤八、采取决定系数R²和相对均方根误差rRMSE评价模型拟合的效果和精度，计算公式如下：

式中，xi为目标区域土壤含水率实测值；

为目标区域土壤含水率平均值；

为目标区域土壤含水率的模型估测值；n为样地的数量；i为目标区域单一样地。

Images