本发明涉及一种小麦氮浓度反演方法，尤其涉及一种基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，属于遥感反演技术领域。
背景技术：
遥感技术常用来实现小麦植株氮浓度的快速、无损探测。长期以来遥感信息的获取手段主要依托地面设备、有人机、卫星三种平台。基于它们获得的数据，学者们开展了大量的研究，具体可分为经验方法、半经验方法和机理模型法三类。
经验方法指基于植被原始光谱或其导数光谱，利用偏最小二乘回归分析、人工神经网络分析及支持向量机分析等提取信息能力强的数学方法，直接建立光谱信息与目标参数的关系模型，利用模型对目标参数进行预测的方法。这种方法往往在建模时能获得较好的结果，但是由于缺乏机理性，经常只适合研究者当时的研究情况，很难在不同地点甚至同一地点不同年份推广应用。
半经验方法指光谱指数法，它通过光谱指数来建立对目标参数的定量关系模型，然后利用模型来对目标参数进行预测，是目前使用最多的方法。光谱指数由植被反射光谱曲线上某些特定波段反射率组合而成，波段组合的选取参照一定的物理基础，能部分消除环境背景的影响，提高对目标参数的敏感性，具有一定的机理性解释。
机理模型法借助叶片光学模型与冠层结构模拟模型的结合来实现作物生理、生化参数的反演。但实际上，目前并没有以植株氮浓度为输入的机理模型，研究者只是基于氮与叶绿素之间的经验关系，将叶绿素对光的吸收转化为氮对光的吸收，从而对机理模型进行改进，来实现对氮素的反演。
在过去的研究中，由于相对于作物个体大小，获得的大部分遥感数据存在着空间分辨率不足的问题，因此，以上三种方法都是只基于遥感影像的光谱信息来反演植株氮浓度的。单一依赖光谱特征的方法，容易受“同物异谱、异物同谱”现象的影响；并且，在小麦生长后期，田间覆盖较高、茂密较大情况下，容易出现“饱和”现象。
低空无人机获得的影像具有超高的空间分辨率，其分辨率能达到厘米乃至毫米级，富含丰富的纹理信息。但是目前，学者们在基于无人机影像反演小麦植株氮浓度时，还主要应用光谱信息，罕有关于纹理信息应用于小麦植株氮浓度反演的报道。
技术实现要素：
为了解决上述技术所存在的不足之处，本发明提供了一种基于无人机影像光谱与纹理特征的小麦氮浓度反演方法。
为了解决以上技术问题，本发明采用的技术方案是：一种基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，方法包括以下步骤：
ⅰ、无人机观测系统搭建与影像获取；无人机影像在10:00～14:00之间晴朗无云时进行获取；影像至少具有722nm、702nm、670nm、550nm、450nm五个波段光谱信息；影像地面分辨率为0.8cm；
ⅱ、影像中土壤信息的自动剔除；利用影像的670nm、550nm、450nm三个波段来计算过绿指数exg和过红指数exr的差值exg-exr；其中，exg、exr的计算公式，如公式①和公式②所示；分析exg-exr的灰度影像，基于otsu方法自动确定影像分割阙值，取大于阙值的像元作为小麦区域，小于阙值的像元作为非小麦区域；去除非小麦区域的影像，仅保留小麦区域影像；
exg＝2×r550-r670-r450公式①
式中，exg为过绿指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率；r450为450nm处的波段反射率；
exr＝1.4×r670-r550公式②
式中，exr为过红指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率；
ⅲ、诊断植株氮浓度的纹理指数构建；光谱指数中的双峰冠层氮指数具有最好反演小麦作物植株氮浓度的能力，因此应用双峰冠层氮指数的纹理特征构建诊断小麦植株氮浓度的纹理指数；
双峰冠层氮指数图像基于步骤ⅱ中剔除非小麦区域的无人机影像计算获得，其计算方法如公式③所示；然后应用灰度共生矩阵法计算双峰冠层氮指数图像的纹理特征；基于获得的灰度共生矩阵，计算纹理特征的同质性和对比度；构建诊断植株氮浓度的纹理指数nti，如公式④所示；
dcni＝(r722-r702)/(r702-r670)/(r722-r670 0.03)公式③
其中，dcni表示双峰冠层氮指数，r722为722nm处的波段反射率；r702为702nm处的波段反射率；r670为670nm处的波段反射率；
其中，nti为诊断植株氮浓度的纹理指数；homogeneitydcni为基于dcni图像计算的同质性；contrastdcni为基于dcni图像计算的对比度；小麦氮素营养越缺乏，dcni图像的同质性会变小，对比度会变大，nti值变小；因此nti与植株氮浓度呈正相关；
ⅳ、耦合纹理指数与光谱指数的植株氮浓度诊断模型构建；耦合诊断植株氮浓度的纹理指数nti和已有诊断植株氮浓度的光谱指数双峰冠层氮指数dcni，采用多元线性回归法建立植株氮浓度的诊断模型，用于小麦植株氮浓度预测，如公式⑤所示；
n％＝a×nti b×dcni c公式⑤
其中，n表示植株氮浓度的诊断模型；a，b，c均为模型系数，由地区决定。
进一步地，步骤ⅲ中灰度共生矩阵描述的是在θ方向上，相隔d像元距离的一对像元分别具有灰度级i和j出现的概率，如公式⑥所示：
式中，glcm为灰度共生矩阵；aij为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数；
计算灰度共生矩阵时，选择的窗口大小为整个目标区域；灰度分级数设为10，即将目标区域dcni图像最大值和最小值之间10等分，并按值大小命名为从1到10共十级，像元dcni值落入哪一等分内，即为哪一级灰度值；搜索的方向设定为0°、45°、90°以及135°四个方向，步长为1个像元。
进一步地，步骤ⅲ中计算对比度的方法如公式⑦所示：
contrastdcni＝∑i,j|i-j|2ai,j公式⑦
式中，contrastdcni为dcni图像的对比度；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
进一步地，步骤ⅲ中计算同质性的方法如公式⑧所示：
式中，homogeneitydcni为dcni图像的同质性；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
本发明摒弃了传统小麦植株氮浓度反演中仅利用光谱信息的缺陷，提出可基于无人机高空间分辨率影像纹理特征来诊断植株氮浓度，并结合已有诊断植株氮浓度的光谱指数，提出用于诊断植株氮浓度的纹理指数；本发明提出了耦合影像光谱、纹理特征开展小麦植株氮浓度反演的技术方法。与现有方法相比，本发明充分利用了低空无人机影像具有高空间分辨率的特征，引入影像纹理特征，耦合光谱和纹理开展小麦植株氮浓度反演，减弱“同物异谱，异物同谱”、“饱和”等现象带来的反演误差，提高植株氮浓度反演精度。
附图说明
图1为本发明的整体流程示意图。
具体实施方式
下面结合具体实施方式对本发明作进一步详细的说明。
图1所示的一种基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，方法包括以下步骤：
ⅰ、无人机观测系统搭建与影像获取；基于大疆m600pro搭载德国cubert公司生产的s185组成无人机观测系统。无人机影像在10:00～14:00之间晴朗无云时进行获取；无人机的飞行高度为30m，航向、旁向的重叠率为75％；影像对应地面分辨率为0.8cm；光谱范围为450～950nm，采样间隔为4nm，光谱分辨率为8nm，共具有125个波段；经检验，s185获取影像的质量具有较高的精度。其中，影像的722nm、702nm、670nm、550nm以及450nm五个光谱波段应用到本发明中。
ⅱ、影像中土壤信息的自动剔除；为了消除土壤背景信息的干扰，用来计算纹理特征的影像需要去除土壤背景。本发明在去除影像中土壤信息的时候采用了otsu自动阙值的方法。首先利用影像的670nm(红光波段)、550nm(绿光波段)、450nm(蓝光波段)三个波段来计算过绿指数exg和过红指数exr的差值exg-exr；其中，exg、exr的计算公式，如公式①和公式②所示；分析exg-exr的灰度影像，基于otsu方法自动确定影像分割阙值，取大于阙值的像元作为小麦区域，小于阙值的像元作为非小麦区域；去除非小麦区域的影像，仅保留小麦区域影像；
exg＝2×r550-r670-r450公式①
式中，exg为过绿指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率；r450为450nm处的波段反射率；
exr＝1.4×r670-r550公式②
式中，exr为过红指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率。
ⅲ、诊断植株氮浓度的纹理指数构建；从叶片尺度看，植物缺氮表现为叶片边缘发黄，整片叶子黄绿相间，叶子变黄比例随植株缺氮程度的增加而增大；从冠层尺度看，植物缺氮表现为从下部往上部叶片逐渐表现为缺氮症状(失绿)，因此从上方俯视看冠层，随着植株缺氮症状的加剧，冠层颜色的对比度逐渐加大，同质性减弱。
此外，光谱指数中的双峰冠层氮指数具有最好反演小麦作物植株氮浓度的能力，因此应用双峰冠层氮指数的纹理特征构建诊断小麦植株氮浓度的纹理指数。
双峰冠层氮指数图像基于步骤ⅱ中剔除非小麦区域的无人机影像计算获得，其计算方法如公式③所示。
dcni＝(r722-r702)/(r702-r670)/(r722-r670 0.03)公式③
其中，dcni表示双峰冠层氮指数，r722为722nm处的波段反射率；r702为702nm处的波段反射率；r670为670nm处的波段反射率。
然后应用灰度共生矩阵法计算双峰冠层氮指数图像的纹理特征；灰度共生矩阵描述的是在θ方向上，相隔d像元距离的一对像元分别具有灰度级i和j出现的概率，如公式⑥所示；
式中，glcm为灰度共生矩阵；aij为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数；
从灰度共生矩阵计算的常用纹理信息包括对比度、能量值、同质性、平均值、熵、方差、二阶距、相关性等。本发明计算灰度共生矩阵时，选择的窗口大小为整个目标区域(如：某一变量施肥单元区)；灰度分级数设为10，即将目标区域dcni图像最大值和最小值之间10等分，并按值大小命名为从1到10共十级，像元dcni值落入哪一等分内，即为哪一级灰度值；搜索的方向设定为0°、45°、90°以及135°四个方向，步长为1个像元。
基于获得的灰度共生矩阵，计算纹理特征的同质性和对比度；选择这两个参数的原因是根据前面所述叶片、冠层在缺氮时其对比度会加大，同质性会变弱。
计算对比度的方法如公式⑦所示：
contrastdcni＝∑i,j|i-j|2ai,j公式⑦
式中，contrastdcni为dcni图像的对比度；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
步骤ⅲ中计算同质性的方法如公式⑧所示：
式中，homogeneitydcni为dcni图像的同质性；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
构建诊断植株氮浓度的纹理指数nti，如公式④所示；
其中，nti为诊断植株氮浓度的纹理指数；homogeneitydcni为基于dcni图像计算的同质性；contrastdcni为基于dcni图像计算的对比度；小麦氮素营养越缺乏，dcni图像的同质性会变小，对比度会变大，nti值变小；因此nti与植株氮浓度呈正相关。
ⅳ、耦合纹理指数与光谱指数的植株氮浓度诊断模型构建；耦合诊断植株氮浓度的纹理指数nti和已有诊断植株氮浓度的光谱指数双峰冠层氮指数dcni，采用多元线性回归法建立植株氮浓度的诊断模型，用于小麦植株氮浓度预测，如公式⑤所示；考虑到不同地区小麦品种、气候的差异，不同地区公式中的模型系数会有所不同，其值由不同区域采集的训练样本来确定。
n％＝a×nti b×dcni c公式⑤
其中，n表示植株氮浓度的诊断模型；a，b，c均为模型系数，由地区决定。
上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。
1.一种基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，其特征在于：所述方法包括以下步骤：
ⅰ、无人机观测系统搭建与影像获取；无人机影像在10:00～14:00之间晴朗无云时进行获取；影像至少具有722nm、702nm、670nm、550nm、450nm五个波段光谱信息；影像地面分辨率为0.8cm；
ⅱ、影像中土壤信息的自动剔除；利用影像的670nm、550nm、450nm三个波段来计算过绿指数exg和过红指数exr的差值exg-exr；其中，exg、exr的计算公式，如公式①和公式②所示；分析exg-exr的灰度影像，基于otsu方法自动确定影像分割阙值，取大于阙值的像元作为小麦区域，小于阙值的像元作为非小麦区域；去除非小麦区域的影像，仅保留小麦区域影像；
exg＝2×r550-r670-r450公式①
式中，exg为过绿指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率；r450为450nm处的波段反射率；
exr＝1.4×r670-r550公式②
式中，exr为过红指数；r550为550nm处的波段反射率；r670为670nm处的波段反射率；
ⅲ、诊断植株氮浓度的纹理指数构建；光谱指数中的双峰冠层氮指数具有最好反演小麦植株氮浓度的能力，因此应用双峰冠层氮指数的纹理特征构建诊断小麦植株氮浓度的纹理指数；
双峰冠层氮指数图像基于步骤ⅱ中剔除非小麦区域的无人机影像计算获得，其计算方法如公式③所示；然后应用灰度共生矩阵法计算双峰冠层氮指数图像的纹理特征；基于获得的灰度共生矩阵，计算纹理特征的同质性和对比度；构建诊断植株氮浓度的纹理指数nti，如公式④所示；
dcni＝(r722-r702)/(r702-r670)/(r722-r670 0.03)公式③
其中，dcni表示双峰冠层氮指数，r722为722nm处的波段反射率；r702为702nm处的波段反射率；r670为670nm处的波段反射率；
其中，nti为诊断植株氮浓度的纹理指数；homogeneitydcni为基于dcni图像计算的同质性；contrastdcni为基于dcni图像计算的对比度；小麦氮素营养越缺乏，dcni图像的同质性会变小，对比度会变大，nti值变小；因此nti与植株氮浓度呈正相关；
ⅳ、耦合纹理指数与光谱指数的植株氮浓度诊断模型构建；耦合诊断植株氮浓度的纹理指数nti和已有诊断植株氮浓度的光谱指数双峰冠层氮指数dcni，采用多元线性回归法建立植株氮浓度的诊断模型，用于小麦植株氮浓度预测，如公式⑤所示；
n％＝a×nti b×dcni c公式⑤
其中，n表示植株氮浓度的诊断模型；a，b，c均为模型系数，由地区决定。
2.根据权利要求1所述的基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，其特征在于：所述步骤ⅲ中灰度共生矩阵描述的是在θ方向上，相隔d像元距离的一对像元分别具有灰度级i和j出现的概率，如公式⑥所示：
式中，glcm为灰度共生矩阵；aij为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数；
计算灰度共生矩阵时，选择的窗口大小为整个目标区域；灰度分级数设为10，即将目标区域dcni图像最大值和最小值之间10等分，并按值大小命名为从1到10共十级，像元dcni值落入哪一等分内，即为哪一级灰度值；搜索的方向设定为0°、45°、90°以及135°四个方向，步长为1个像元。
3.根据权利要求2所述的基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，其特征在于：所述步骤ⅲ中计算对比度的方法如公式⑦所示：
contrastdcni＝∑i,j|i-j|2ai,j公式⑦
式中，contrastdcni为dcni图像的对比度；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
4.根据权利要求2所述的基于无人机影像光谱与纹理特征的小麦氮浓度反演方法，其特征在于：所述步骤ⅲ中计算同质性的方法如公式⑧所示：
式中，homogeneitydcmi为dcni图像的同质性；ai,j为灰度级为i的像元与灰度级为j的像元同时出现的概率；i，j的取值为1到n，n为总的灰度分级数。
技术总结