一种基于卫星遥感影像的叶龄空间分布检测方法

技术领域

本发明属于农业遥感应用领域，具体涉及一种水稻叶龄诊断方法。

背景技术

水稻生育叶龄诊断技术是在水稻“器官同伸理论”与“叶龄模式理论”的基础上发展起来的水稻栽培新技术，使传统的种管收流程式栽培技术发展成为按叶龄诊断、预测、调控的栽培技术体系。根据水稻某时的叶龄进程，可推测以后一段时间的叶龄进程，从而推测出幼穗分化、拔节、减数分裂、抽穗等关键时期。传统的水稻叶龄诊断方法有：叶龄跟踪方法、种谷方向法、叶脉偏向法、伸长叶枕距法、变形叶鞘法、最长叶法、茎数调查法等。

已有叶龄诊断专利包括：一种用于记录水稻叶龄的标记尺(CN217032268U)，通过标记尺来测量最上端的叶片与下面的叶片长度是否一样；一种水稻叶龄诊断方法(CN114913217A)，获取所述水稻当前叶的长度和前一完全叶的长度以及获取所述水稻所处环境的温度，计算所述水稻当前叶的长度和前一叶的完全长度的比值，根据所述环境的温度和所述比值，确定所述水稻当前叶的叶龄。

以上发明虽然可以进行水稻叶龄的检测，但需要人力实地进行观测，耗费人力较大。

现有技术中还有一种水稻叶龄指标管理数字化诊断无人车(CN114080889A)，无人车能够对叶龄进行数据采集，能够进行叶片采样收集；一种叶龄识别方法、装置、电子设备及可读存储介质(CN112800845A)，获得并保存初始状态下目标子区域内植株的初始叶龄，之后每间隔预设时长获得目标子区域内植株图像，并对植株图像进行心叶识别，获得心叶数量(植株主茎上)，根据心叶数量对已保存的距离当前最近的目标历史叶龄进行更新，获得植株的当前叶龄。

以上发明虽然可以进行水稻叶龄的检测，但只适用于小范围地块。

发明内容

本发明的目的是提供一种基于卫星遥感影像的叶龄空间分布检测方法，借助卫星影像计算稻田的叶龄空间分布，为后续按照叶龄进行水稻田间管理提供数据基础。改善现有技术只适用于小范围地块，且耗时耗力的问题。

本发明所述的一种基于卫星遥感影像诊断叶龄的方法包括如下步骤：

S1，对卫星遥感影像进行大气校正、去云预处理；

S2，根据S1中经过大气校正，去云预处理的卫星遥感影像计算归一化植被指数NDVI；

S3，提取叶龄待测区域对应的NDVI时间序列；

S4,对S3中获得的NDVI时间序列进行S-G滤波处理，降低噪声的影响；

S5,建立S4中的降噪后的NDVI和对应的实测区域叶龄的回归方程y＝ax+b，叶龄x代表NDVI，y代表叶龄，a和b与水稻品种、水稻最大叶数等有关；

S6，基于S5获得的回归方程，根据实时采集的卫星遥感影像，获得所述影像中待检测叶龄的区域内的叶龄空间分布。

优选实施方式，所述卫星遥感影像是待检测叶龄的水稻种植区域的卫星遥感影像。

优选实施方式，所述步骤S1包括如下步骤：

S11辐射定标：将卫星遥感影像中卫星传感器记录下的数字信号值转化为地物反射太阳光线的辐射亮度值或是反射率，所述辐射亮度值为：

L

式中，L

所述反射率为：

式中，ρ

S12大气校正：根据大气辐射传输模型进行遥感影像的大气校正；

S13去云：使用卫星遥感影像中的QA波段去云。

优选实施方式，所述步骤S2包括如下步骤：

S21通过归一化植被指数NDVI量化植被；归一化植被指数NDVI公式如下，

其中，NIR是指该影像像元的近红外波段值，RED是指该影像像元的红波段值。

优选实施方式，所述S21包括如下步骤：计算每一个影像像元的归一化植被指数NDVI。

优选实施方式，所述S-G滤波处理是基于最小二乘原理的多项式平滑算法，是应用加权滑动平均窗口进行滤波，在滑动窗口中其权重值取决于在一个窗口内做最小二乘拟合的多项式系数。

本发明还要求保护一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行本发明所述的任意一种基于卫星遥感影像的叶龄空间分布检测方法。

本发明还包括一种基于卫星遥感影像的叶龄空间分布检测装置，所述装置包括处理器和存储介质，所述存储介质中存入有计算机程序，所述计算机程序被处理器运行时执行本发明所述的任意一种基于卫星遥感影像的叶龄空间分布检测方法。

本发明的有益效果：

卫星遥感具有探测范围广，获取资料的速度快、周期短，获取的信息量大，受地面条件限制少等优势，是目前很多现有技术采用的一种图像获取技术手段。但是，卫星遥感获得的图像信息也存在的分辨率过低的缺陷，导致其无法应用于需要图像分辨率较高的技术领域中。

现有技术中，田间拍摄水稻照片多为厘米及毫米级尺度，在该尺度下，叶长、叶宽、叶子表面的虫洞都可以观察到，因此可以基于计算机视觉方法，综合利用颜色、纹理等信息进行叶龄诊断，即：该技术对所述水稻照片的图片分辨率有一定要求，必须能够清晰的识别叶片的具体形状、颜色和纹理。而采用卫星遥感技术获得的图像数据，即使是亚米级高分辨率数据，分辨率也在30cm以上，无法识别每一株水稻的叶长、叶宽、叶片表面特征。

鉴于上述情况，由于水稻叶龄检测技术对水稻叶片图像的分辨率的较高的要求，导致没有人会去思考采用分辨率较低的卫星影像去实现对水稻叶龄的检测。本发明，采用了卫星影像实现了水稻叶龄的检测克服了由于其分辨率过低导致无法实现叶龄检测的缺陷，通过对卫星遥感影像的处理、并构建叶龄和NDVI的关系模型，快速地获得大范围稻田的叶龄空间分布，实现了叶龄诊断，所述测量方法与现有技术相比，在节省时间和人力的前提下，能够快速获得大范围的叶龄检测结果。

本发明所述方法与现有方法相比较，节约了大量的人力，降低了人工成本，且适用于大范围的区域性叶龄分布的检测，该方法地应用于将有利于推动使用叶龄指导水稻生产管理，确保各项生产技术措施及时性、准确性，对水稻生产的优质、高产、高效有着重要的现实意义。

本发明所述的方法应用于大面积的水稻种植管理中。还可以推广到玉米、马铃薯、小麦等作物的叶龄分布检测中。

附图说明

图1是实施方式六所示的一种基于卫星遥感影像的叶龄空间分布检测方法的流程图。

图2是实施方式五中提及的S-G滤波处理的原理示意图。

图3是实施方式六中提及的，针对闫家岗农场卫星NDVI与叶龄数据关系图。图中圆点表示原始数据；

点画线代表决定系数R

短横划线代表决定系数R

实线代表决定系数R

长横划线代表决定系数R

图4是实施方式六所述的，查哈阳农场卫星NDVI与叶龄数据关系图。图中圆点表示原始数据；

点画线代表决定系数为R

短横划线代表决定系数R

实线代表决定系数R

长横划线代表决定系数R

具体实施方式

下面将结合附图对本发明做进一步详细说明。

实施方式一.本实施方式所述的一种基于卫星遥感影像的叶龄空间分布检测的方法包括如下步骤：

S1，对卫星遥感影像进行大气校正、去云预处理；

S2，根据S1中经过大气校正，去云预处理的卫星遥感影像计算归一化植被指数NDVI；

S3，提取叶龄待测区域对应的NDVI时间序列；

S4,对S3中获得的NDVI时间序列进行S-G滤波处理，降低噪声的影响；

S5,建立S4中的降噪后的NDVI和对应的实测叶龄数据的回归方程y＝ax+b，叶龄x代表NDVI，y代表叶龄，a和b与水稻品种、水稻最大叶数等有关；

S6，基于S5获得的回归方程，根据实时采集的卫星遥感影像，获得所述影像中待检测叶龄的区域内的叶龄空间分布。

本实施方式中，所述的实测叶龄数据，是提前获得的数据，可以采用任何一种现有技术进行对所述区域的水稻的叶龄进行实际测量获得，测量的是上一个水稻生长周期数据，叶龄与时间、地理位置相对应。所述叶龄数据还可以通过读取之前的叶龄数据库获得。

本实施方式所述的卫星遥感影像是包含有待检测叶龄的水稻种植区域的卫星遥感影像。

实施方式二.本实施方式是对实施方式一所述的一种基于卫星遥感影像的叶龄空间分布检测方法的优化，本实施方式中，所述步骤S1包括如下步骤：

S11辐射定标：将卫星遥感影像中卫星传感器记录下的数字信号值转化为地物反射太阳光线的辐射亮度值或是反射率，

L

式中，L

式中，ρ

S12大气校正：根据大气辐射传输模型进行遥感影像的大气校正；

使用大气校正方法Sen2cor进行进行遥感影像的大气校正。所述大气校正方法Sen2cor的算法核心是大气辐射传输模型libRadtran。所述大气辐射传输模型libRadtran用于计算地球大气辐射传输，辐射方程可以表示为，

/>

式中，L为位置x,y,z的辐射亮度，β为体积消光系数,ds为光线传播方向上通过的λ

厚度，J为源函数，可定义为，

λ

式中，ω为单次散射反照率，p(Ω,Ω')为从方向Ω'到Ω的相位函数，B(T)为普朗克函数。

所述大气校正方法Sen2cor是将大气辐射传输模型libRadtran代码生成的大气校正查找表集成在Sen2cor程序中。大气校正查找表存储的反射率的光谱分辨率为0.6nm，利用滤波函数对光谱进行重采样，得到特定通道的反射率。

S13去云：使用卫星遥感影像中的QA波段去云。

实施方式三.本实施方式是对实施方式一所述的一种基于卫星遥感影像的叶龄空间分布检测方法的优化，本实施方式中，所述步骤S2包括如下步骤：

S21通过归一化植被指数NDVI量化植被；归一化植被指数NDVI公式如下，

其中，NIR是指该影像像元的近红外波段值，RED是指该影像像元的红波段值。

实施方式四.本实施方式是对实施方式三所述的一种基于卫星遥感影像的叶龄空间分布检测方法的优化，本实施方式中，所述S21包括如下步骤：计算每一个影像像元的归一化植被指数NDVI。

实施方式五.本实施方式是对实施方式一所述的一种基于卫星遥感影像的叶龄空间分布检测方法的优化，本实施方式中，参见图2所示，所述S-G滤波处理是基于最小二乘原理的多项式平滑算法，是应用加权滑动平均窗口进行滤波，在滑动窗口中其权重值取决于在一个窗口内做最小二乘拟合的多项式系数。

实施方式六.参见图1说明本实施方式。本实施方式是采用本发明所述方法应用于实际场景的叶龄空间分布检测方法：

确定某四地农场的水稻叶龄实测点的位置，从移栽后直到孕穗期以前，每隔两天人工测量一次目标点水稻的叶龄，汇总整理实测的叶龄值。同时要收集品种、插秧时间、采集时间等信息。

选择叶龄实测点位置处、叶龄监测时期的哨兵2号影像。哨兵2号卫星5天即可实现重复观测，一景哨兵2号影像覆盖范围是10000平方公里，可见光和近红外波段的影像像素单元是10m，一景L2A级别产品占用空间大约是800MB，影像数据量大，数据预处理和后续处理分析的计算量也大，借助Google earth engine云平台进行遥感影像的处理会非常便利。GEE将地理数据分析算法封装成函数接口的形式提供给用户，通过互联网可以使用户快速地访问海量的遥感影像以及其他空间数据，在云端实现数据处理分析的全部流程。

处理卫星影像得到地物的光谱反射率信息之前，需要进行辐射定标、大气校正、去云等预处理。辐射定标是将卫星传感器记录下的数字信号值转化为地物反射太阳光线的辐射亮度值或是反射率，其中，辐照亮度为：

L

反射率为，

大气校正是为消除气溶胶、水汽、气体分子等对太阳光线散射和吸收的影响，使卫星接收的辐射亮度全部来自于地物。Sen2cor是欧洲航天局(ESA)官方提供的哨兵2号卫星影像的大气校正方法，算法的核心是大气辐射传输模型libRadtran。libRadtran是用于计算地球大气辐射传输的软件包，该算法运用离散纵标法计算地球大气、太阳以及热辐射传输，适用于各种大气条件。辐射方程可以表示为，

式中，L为位置(x,y,z)的辐亮度，β为体积消光系数,ds为光线传播方向上通过的λ

厚度，J为源函数，可定义为，

λ

式中，ω为单次散射反照率，p(Ω,Ω')为从方向Ω'到Ω的相位函数，B(T)为普朗克函数。

Sen2cor将libRadtran代码生成的大气校正查找表集成在程序中，缩短了程序运行的时间并简化了操作流程。该查找表包含了各种大气条件、太阳角度、地形、4种基本气溶胶类型(农村、城市、沙漠和海洋)以及6种不同的臭氧含量等，存储的反射率的光谱分辨率为0.6nm，随后利用滤波函数对光谱进行重采样，得到特定通道的反射率。

在光学遥感影像上，云遮盖了地物，影像不能反映地物的信息。使用哨兵2号影像的QA60波段或者哨兵2号Cloud Probability数据集去云。QA60可以识别不透明的云和卷云，不透明的云在蓝波段(B2)具有高反射率，为了避免雪的影响，还使用了云的反射率很高而雪的反射率较低的SWIR波段(B11和B12)；B10对应高大气吸收带，只能检测到高空云，且卷云是半透明的，在B2蓝色波段无法检测到，因此B2反射率低、B10波段反射率高的像素很有可能是卷云。QA60波段二进制数第10位数代表不透明的云，若0则没有不透明的云，若1存在不透明的云，第11位数代表卷云，若0则无卷云，若1存在卷云。哨兵2号CloudProbability数据集是使用sentinel2-cloud-detector库(使用LightGBM)创建的，在应用梯度提升基础算法之前，使用双线性插值将所有波段上采样到10m分辨率，较高的数值更可能是云或高反射表面(例如屋顶或雪)。

植被指数是利用遥感影像进行波段计算的结果，目的是为了突出植被的信息。NDVI通过测量近红外波段(植被强烈反射)和红波段(植被吸收)之间的差异来量化植被。NDVI公式如下，

其中，NDVI是指该影像像元的归一化植被指数，NIR是指该影像像元的近红外波段值，RED是指该像元的红波段值。

NDVI的范围始终为-1至+1，但每种地覆盖并没有明确的界限。例如，当NDVI值为负数时，可能是水；NDVI值接近1，可能是茂密的绿叶；当NDVI接近于零时，没有绿叶，可能是城市化区域。

NDVI是公认表征植被变化最有效的参数之一，可较好地反映植被绿度变化，指示植被生长状态及植被覆盖度。使用不同时间的NDVI数据来指示水稻的生长发育状态，而水稻生育进程可以用叶龄来定量表征，基于此尝试建立NDVI与水稻叶龄的关系。

受云、冰雪、气溶胶等因素的干扰，NDVI时间序列数据包含大量噪声，使用S-G滤波处理时间序列数据，降低噪声的影响。S-G滤波是基于最小二乘原理的多项式平滑算法，是应用加权滑动平均窗口进行滤波，在滑动窗口中其权重值取决于在一个窗口内做最小二乘拟合的多项式系数。S-G滤波的滤波效果主要由两个参数控制，分别是滤波窗口大小和平滑多项式的次数，滤波窗口的宽度越大，多项式的次数越低，拟合的NDVI时间序列越平滑。

S-G滤波平滑过程如图2，设定滤波窗口包括五个点，使用x和其前后两个点的多项式拟合值来取代x，然后依次移动，直到把把所有的点遍历完。

S-G滤波平滑过程矩阵算子的求解如下，假设滤波窗口的宽度为n＝2m+1,各测量点x＝(-m+0,-m+1,-m+2..,m-1,m)，采用k-1次多项式对窗口内的数据点进行拟合，

y＝a

得到n个这样的方程，组成k元线性方程组，要使方程组有解，应该满足n大于等于k，一般选择n>k，通过最小二乘程法确定参数解A。由此得到，

用矩阵表示为，

Y

A的最小二乘解

Y的滤波值为

B＝X·(X

在滤波平滑过程中只更新第m个点的坐标，故在得到B矩阵后，只取B矩阵的第m行,并与窗口内的2m+1个值使用多项式拟合，从而更新窗口内的第m个点，然后窗口移动依次遍历每一个点。

建立NDVI和实测叶龄的一元线性回归方程(y＝ax+b)，基于回归方程计算稻田叶龄的空间分布。

闫家岗农场共有10个无人机地块叶龄采样点，采样时间为5月20日到7月23日，水稻品种是育品16(最大叶数13叶)。共收集到50个叶龄采样值对应的NDVI数据。由闫家岗农场卫星NDVI与叶龄数据关系图(图3)可知，当水稻品种是育品16时(最大叶数13叶)，得到的一元线性回归方程为y＝13.043x+4.0361，决定系数R2＝0.89，均方根误差RMSE＝0.92。

查哈阳农场共有10个无人机地块叶龄采样点，采样时间为6月19日到7月17日，水稻品种是寒稻608(最大叶数11叶)。共收集到27个叶龄采样值对应的NDVI数据，并进行了S-G滤波降噪处理。由查哈阳农场卫星NDVI与叶龄数据关系图(图4)可知，当水稻品种是寒稻608时(最大叶数11叶)，得到的一元线性回归方程为y＝14.838x+1.9854，决定系数R2＝0.72，均方根误差RMSE＝0.84。

叶龄采样点数据量增加后，还可以将数据集分为训练集和测试集，避免回归模型过拟合的问题。

上述实施方式是对本发明所述基于卫星遥感影像的叶龄空间分布检测的方法的举例说明，但不限制本发明的保护范围。