本发明属于遥感数据分析技术领域，涉及一种基于遥感影像植被指数的土壤生产力等级识别方法。本发明通过遥感的方法解决区域尺度上大范围土壤生产力等级识别与可视化难题，通过卫星遥感影像的植被指数与作物产量建立遥感估产模型，通过作物产量与作物无肥区产量建立产量关系模型，通过作物无肥区产量与土壤性质建立土壤生产力等级模型，其中，土壤生产力等级根据无肥区产量划分；融合三种模型进一步建立卫星遥感影像的植被指数与土壤生产力等级和土壤性质的关联，从而实现通过卫星遥感影像来大范围识别土壤生产力等级以及仿真相应的
1.一种基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述方法包括如下步骤：S1、构建研究区数据库；
获取研究区内目标作物的物候信息、研究区内目标作物关键生育期的多种卫星的遥感影像数据、研究区内目标作物的施肥区产量数据、无肥区产量数据和对应的土壤数据；所述土壤数据包含土壤全氮、有机质、有效磷、有效钾和pH；
S2.1、分别提取步骤S1获得的研究区内目标作物关键生育期的各卫星遥感影像中的光谱特征，获得光谱特征参数图层，所述光谱特征参数图层包括近红外波段反射率ρNIR参数图层、红光波段反射率ρR参数图层、绿光波段反射率ρG参数图层和蓝光波段反射率ρB参数图层；
S2.2、根据步骤S2.1的各光谱特征参数图层，计算植被指数特征，获得植被指数特征图层，所述植被指数特征图层包括归一化植被指数NDVI图层、差值植被指数DVI图层、垂直植被指数PVI图层、土壤调节植被指数SAVI图层和植被覆盖度VFC图层；
S2.3、根据步骤S2.1的各光谱特征参数图层，任意选取两个不同特征时相的同一波段的光谱特征进行差值计算，获得光谱差值特征图层，所述光谱差值特征图层包括近红外波段反射率差值图层、红光波段反射率差值图层、绿光波段反射率差值图层和蓝光波段反射率差值图层；
S2.4、根据步骤S2.2的各植被指数特征图层，选取与步骤S2.3相同的两个特征时相的同一植被指数特征进行差值计算，获得植被指数差值特征图层，所述植被指数差值特征图层包括归一化植被指数差值图层、差值植被指数差值图层、垂直植被指数差值图层、土壤调节植被指数差值图层和植被覆盖度差值图层；
将步骤S2获得的所有图层合并为一个图层；并将该图层与已知的目标作物样本作为输入集，运行随机森林算法得到研究区目标作物种植区域；其中，所述目标作物样本即确定种植目标作物的样本的坐标信息；
通过步骤S2获得的研究区目标作物关键生育期的各卫星遥感影像中的归一化植被指数NDVI与步骤S1获得的目标作物施肥区产量数据分别构建线性模型、多项式模型、指数模型、对数模型和乘幂模型，并通过模型评价指标进行最优模型筛选，所述模型评价指标包括
决定系数R、相对误差R和均方根误差RMSE，依据决定系数R越大且相对误差R和均方根误差RMSE越小的原则确定最优的研究区目标作物的遥感估产模型；根据最优的研究区目标作物的遥感估产模型确定预测目标作物产量最佳的卫星遥感影像数据和最佳生育期；
S5.1、根据步骤S1获得的目标作物施肥区产量数据和目标作物无肥区产量数据建立产量关系模型；
S5.2、根据步骤S1获得的目标作物无肥区产量数据划分土壤基础生产力，并与土壤数据建立土壤生产力等级模型；
S5.3、通过步骤S4确定的最优的研究区目标作物的遥感估产模型、步骤S5.1建立的产量关系模型和步骤S5.2建立的土壤生产力等级模型建立研究区目标作物关键生育期卫星遥感影像中的归一化植被指数NDVI与土壤生产力等级的关系模型，进而建立土壤性质仿真；
S5.4、根据步骤S4确定的最佳的卫星遥感影像数据和最佳生育期在步骤S3确定的研究区目标作物种植区域内提取任一点的归一化植被指数NDVI，借助步骤S5.3获得的研究区目标作物关键生育期卫星遥感影像中的归一化植被指数NDVI与土壤生产力等级的关系模型确定该点的土壤生产力等级并仿真出该点的土壤性质。
2.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S1中，卫星包括MODIS、SPOT和MERIS卫星。
3.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S1中，所述研究区内目标作物关键生育期的多种卫星的遥感影像数据根据研究区气象数据筛选无云层遮挡的数据。
4.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S1中，若研究区面积较大，则通过影像拼接的方式拼接遥感影像，并以研究区边界矢量图层为掩膜进行裁剪。
5.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S2.2中，通过公式1～公式5计算植被指数特征；
公式1～公式5中，NDVI为归一化植被指数，ρNIR为近红外波段反射率，ρR为红光波段反射率，DVI为差值植被指数，PVI为垂直植被指数，SR为红光波段土壤反射率，VR为红光波段植被反射率，SNIR为近红外波段土壤反射率，VNIR为近红外波段植被反射率，SAVI为土壤调节植被指数，L为土壤调节系数，取值范围0‑1，L＝0时表示植被覆盖率为零，L＝1时，表示土壤背景的影像为零，即植被覆盖度非常高；VFC为植被覆盖度，NDVIsoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIveg代表完全被植被所覆盖的NDVI值，即纯植被像元的NDVI值。
6.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S4中，线性模型：y＝a×NDVI+b
多项式模型：y＝a×(NDVI) +b×NDVI+c
其中，y代表目标作物产量，NDVI代表目标作物某一关键生育期的归一化植被指数值；
a、b、c代表方程的回归参数，并且不同方程之间的参数值不同。
7.根据权利要求1所述的基于遥感影像植被指数的土壤生产力等级识别方法，其特征在于，所述步骤S4中，2
其中，R为决定系数， Yi、 ESS是回归平方和，TSS是残差平方和；Yi表示参与建模的每个NDVI样本对应的实际产量，表示所有样本的对应实际产量的平均值，表示根据选定模型计算得到的每个样本对应的预测产量；
其中，RMSE为均方根误差；n代表测量次数，Xobs,i代表真实值，Xmodel,i表示模型预测值。
一种基于遥感影像植被指数的土壤生产力等级识别方法
技术领域
[0001] 本发明属于遥感数据分析技术领域，涉及一种基于遥感影像植被指数的土壤生产力等级识别方法。
背景技术
[0002] 土壤是作物生产的基础，作物产量潜力和水肥调控作用的持续稳定发挥依赖于良好的土壤条件。土壤生产力是指一定的资源投入下土壤能够产出人类所需的植物产品的能力，是土壤接纳和储存养分水分，提供作物根系发育环境能力的具体表现。
[0003] 定量化评价土壤生产力，及其在生产中的大面积识别仍然是一个科学难题和挑战。目前关于土壤生产力的评价主要是基于土壤属性，如土壤性质或通过作物农学属性，如作物产量。这些方法主要通过定点采样进行作物产量或土壤养分的测定，适合于静态(单一时相)以及小尺度(田块)评价；而区域土壤生产力及土壤性质表现出很大的异质性；因此，定点采样方法耗时耗力，对于大尺度评价及识别存在诸多困难。
[0004] 卫星遥感具有高时空分辨率、成本低、易获得的特点，是区域大尺度植被动态反演的重要途径。植被指数已被广泛用来定性和定量评价植被覆盖度及其生长状态。由于植被光谱反映的是植被、土壤亮度、环境影响、阴影、土壤颜色和湿度等的混合效应，且受大气空间和时相变化的影响，因此不同的植被指数往往可以表征不同的农学特征。例如，比值植被指数又称为绿度，为二通道反射率之比，能较好地反映植被覆盖度和生长状况的差异，特别适用于植被生长旺盛、具有高覆盖度的植被监测。蓝光、红光和近红外通道的组合可大大消除大气中气溶胶对植被指数的干扰，所组成的抗大气植被指数可大大提高植被长势监测和作物估产精度。最为常见的归一化植被指数(NDVI)为两个通道反射率之差除以它们的和。在植被处于中、低覆盖度时，该指数随覆盖度的增加而迅速增大，当达到一定覆盖度后增长缓慢，所以适用于植被早、中期生长阶段的动态监测，常用于作物的长势监测和产量估计。
目前，国内学者分别利用NOAA‑NDVI与MODIS‑NDVI遥感数据，对四川、湖南、广西、江西和黑龙江建立以省为单位的水稻估产模型和冬小麦产量估算模型，并达到一定的估产精度。美国学者提出一种普适的冬小麦产量回归预测模型，最初模型建立于美国堪萨斯州。随后，该模型在没有经过任何修改和标定的情况下直接应用于欧洲乌克兰地区。在冬小麦成熟期前
6周预测单产，误差在15％左右，预测总产误差10％左右。总之，卫星遥感技术的快速发展，为区域大尺度产量精确估测提供了可能。
发明内容
[0005] 传统土壤生产力等级识别方法主要依赖于采用土壤取样或作物产量实地实测。这种方法虽然识别结果较为准确，但人力物力耗费大、分析周期长，不适用于区域尺度上大范围的土壤生产力等级识别；这对于土壤智慧型的作物精准化管理造成根本性限制。本发明的目的是提供一种基于遥感影像植被指数的土壤生产力等级识别方法，通过遥感的方法解决区域尺度上大范围土壤生产力等级识别与可视化难题，通过卫星遥感影像的植被指数与作物产量建立遥感估产模型，通过作物产量与作物无肥区产量建立产量关系模型，通过作物无肥区产量与土壤性质建立土壤生产力等级模型，其中，土壤生产力等级根据无肥区产量划分；融合三种模型进一步建立卫星遥感影像的植被指数与土壤生产力等级和土壤性质的关联，从而实现通过卫星遥感影像来大范围识别土壤生产力等级以及仿真相应的土壤性质。
[0006] 为了实现上述目的，本发明提供了如下技术方案：
[0007] 一种基于遥感影像植被指数的土壤生产力等级识别方法，包括如下步骤：
[0008] S1、构建研究区数据库；
[0009] 获取研究区内目标作物的物候信息、研究区内目标作物关键生育期的多种卫星的遥感影像数据、研究区内目标作物的施肥区产量数据、无肥区产量数据和对应的土壤数据；所述土壤数据包含土壤全氮、有机质、有效磷、有效钾和pH；
[0010] S2、基于卫星遥感影像的植被指数计算；
[0011] S2.1、分别提取步骤S1获得的研究区内目标作物关键生育期的各卫星遥感影像中的光谱特征，获得光谱特征参数图层，所述光谱特征参数图层包括近红外波段反射率ρNIR参数图层、红光波段反射率ρR参数图层、绿光波段反射率ρG参数图层和蓝光波段反射率ρB参数图层；
[0012] S2.2、根据步骤S2.1的各光谱特征参数图层，计算植被指数特征，获得植被指数特征图层，所述植被指数特征图层包括归一化植被指数 NDVI图层、差值植被指数DVI图层、垂直植被指数PVI图层、土壤调节植被指数SAVI图层和植被覆盖度VFC图层；
[0013] S2.3、根据步骤S2.1的各光谱特征参数图层，任意选取两个不同特征时相的同一波段的光谱特征进行差值计算，获得光谱差值特征图层，所述光谱差值特征图层包括近红外波段反射率差值图层、红光波段反射率差值图层、绿光波段反射率差值图层和蓝光波段反射率差值图层；
[0014] S2.4、根据步骤S2.2的各植被指数特征图层，选取与步骤S2.3相同的两个特征时相的同一植被指数特征进行差值计算，获得植被指数差值特征图层，所述植被指数差值特征图层包括归一化植被指数差值图层、差值植被指数差值图层、垂直植被指数差值图层、土壤调节植被指数差值图层和植被覆盖度差值图层；
[0015] S3、提取研究区目标作物种植区域；
[0016] 将步骤S2获得的所有图层合并为一个图层；并将该图层与已知的目标作物样本作为输入集，运行随机森林算法得到研究区目标作物种植区域；其中，所述目标作物样本即确定种植目标作物的样本的坐标信息；
[0017] S4、建立研究区目标作物的遥感估产模型；
[0018] 通过步骤S2获得的研究区目标作物关键生育期的各卫星遥感影像中的归一化植被指数NDVI与步骤S1获得的目标作物施肥区产量数据分别构建线性模型、多项式模型、指数模型、对数模型和乘幂模型，并通过模型评价指标进行最优模型筛选，所述模型评价指标2 2
包括决定系数R、相对误差R和均方根误差RMSE，依据决定系数R越大且相对误差R和均方根误差RMSE越小的原则确定最优的研究区目标作物的遥感估产模型；根据最优的研究区目标作物的遥感估产模型确定预测目标作物产量最佳的卫星遥感影像数据和最佳生育期；
[0019] S5、研究区土壤生产力等级的识别；
[0020] S5.1、根据步骤S1获得的目标作物施肥区产量数据和目标作物无肥区产量数据建立产量关系模型；
[0021] S5.2、根据步骤S1获得的目标作物无肥区产量数据划分土壤基础生产力，并与土壤数据建立土壤生产力等级模型；
[0022] S5.3、通过步骤S4确定的最优的研究区目标作物的遥感估产模型、步骤S5.1建立的产量关系模型和步骤S5.2建立的土壤生产力等级模型建立研究区目标作物关键生育期卫星遥感影像中的归一化植被指数 NDVI与土壤生产力等级的关系模型，进而建立土壤性质仿真；
[0023] S5.4、根据步骤S4确定的最佳的卫星遥感影像数据和最佳生育期在步骤S3确定的研究区目标作物种植区域内提取任一点的归一化植被指数NDVI，借助步骤S5.3获得的研究区目标作物关键生育期卫星遥感影像中的归一化植被指数NDVI与土壤生产力等级的关系模型确定该点的土壤生产力等级并仿真出该点的土壤性质。
[0024] 所述步骤S1中，卫星包括MODIS、SPOT和MERIS卫星。
[0025] 所述步骤S1中，所述研究区内目标作物关键生育期的多种卫星的遥感影像数据根据研究区气象数据筛选无云层遮挡的数据。
[0026] 所述步骤S1中，若研究区面积较大，则通过影像拼接的方式拼接遥感影像，并以研究区边界矢量图层为掩膜进行裁剪。
[0027] 所述步骤S2.2中，通过公式1～公式5计算植被指数特征；
[0028]
[0029] DVI＝ρNIR‑ρR 公式2
[0030]
[0031]
[0032]
[0033] 公式1～公式5中，NDVI为归一化植被指数，ρNIR为近红外波段反射率，ρR为红光波段反射率，DVI为差值植被指数，PVI为垂直植被指数，SR为红光波段土壤反射率，VR为红光波段植被反射率，bNIR为近红外波段土壤反射率，VNIR为近红外波段植被反射率，SAVI为土壤调节植被指数，L为土壤调节系数，取值范围0‑1，L＝0时表示植被覆盖率为零，L＝1时，表示土壤背景的影像为零，即植被覆盖度非常高；VFC 为植被覆盖度，NDVIsoil为完全是裸土或无植被覆盖区域的NDVI值， NDVIveg代表完全被植被所覆盖的NDVI值，即纯植被像元的NDVI值；
[0034] 所述步骤S4中，
[0035] 线性模型：y＝a×NDVI+b
[0036] 多项式模型：y＝a×(NDVI)2+b×NDVI+c
[0037] 指数模型：y＝a×eb×NDVI
[0038] 对数模型：y＝a×ln(NDVI)+b
[0039] 乘幂模型：y＝a×(NDVI)b
[0040] 其中，y代表目标作物产量，NDVI代表目标作物某一关键生育期的归一化植被指数值；a、b、c代表方程的回归参数，并且不同方程之间的参数值不同。
[0041] 所述步骤S4中，
[0042]2
[0043] 其中，R为决定系数， Yi、 ESS是回归平方和，TSS是残差平方和；Yi表示参与建模的每个NDVI样本对应的实际产量，表示所有样本的对应实际产量的平均值，表示根据选定模型计算得到的每个样本对应的预测产量；
[0044]
[0045] 其中，R为相对误差，Δ为绝对误差，L为实际真值；
[0046]
[0047] 其中，RMSE为均方根误差；n代表测量次数，Xobs,i代表真实值，Xmodel,i表示模型预测值。
[0048] 与现有技术相比，本发明的有益效果在于：
[0049] 本发明基于作物无肥区产量和土壤性质结合的方法来划分土壤生产力等级。根据作物无肥区产量和土壤性质(包含土壤全氮、有机质、有效磷、速效钾和pH)等指标的值域范围，将其划分为高、中、低和极低4类，分别对应土壤生产力等级的4个级别。
[0050] 本发明基于作物关键生育期卫星遥感影像的植被指数识别作物种植区域，以及建立作物产量的遥感反演技术。依据研究区域的作物物候信息确定作物的关键生育期，并在MODIS、SPOT和MERIS等不同卫星的遥感影像中选择处于作物关键生育期内(或时间接近)的遥感影像，分别计算多种植被指数并根据不同地物植被指数的差异识别作物种植区域；通过建立卫星遥感影像的植被指数与作物产量的遥感估产模型，实现作物产量的遥感反演(即通过少量样本建立遥感估产模型后，可实现通过计算卫星遥感影像的植被指数反演整个作物种植区域内的作物产量)。
[0051] 本发明融合作物产量的遥感反演技术和土壤生产力等级划分的区域大范围土壤生产力等级识别和可视化技术。根据作物无肥区产量和作物产量间的产量关系模型，进一步融合遥感估产技术和土壤生产力等级划分方法，实现作物种植区土壤生产力等级识别及土壤性质仿真。
[0052] 综上，本发明建立了卫星遥感影像的植被指数与作物产量、作物产量与作物无肥区产量、作物无肥区产量与土壤性质间的关联，并进一步构建了卫星遥感影像的植被指数与土壤生产力等级和土壤性质间的关联，从而实现了通过卫星遥感反演进行土壤生产力等级识别划分和土壤性质的仿真。
附图说明
[0053] 图1为本发明实施例的华北平原土壤生产力等级识别方法的流程示意图；
[0054] 图2为本发明实施例的华北平原分辨率为1km*1km的遥感影像；
[0055] 图3为本发明实施例的植被指数计算和种植区域提取的流程示意图；
[0056] 图4为本发明实施例提取的华北平原冬小麦种植区域示意图；
[0057] 图5为本发明实施例的冬小麦产量数据(施肥区)和冬小麦无肥区产量数据建立产量关系示意图；
[0058] 图6为本发明实施例建立的华北平原冬小麦关键生育期卫星遥感影像中的归一化植被指数NDVI与土壤生产力等级的关系示意图。
具体实施方式
[0059] 下面结合附图和实施例对本发明进行进一步说明。
[0060] 以华北平原冬小麦为例。
[0061] 如图1所示，华北平原土壤生产力等级识别方法，包括如下步骤：
[0062] S1、华北平原数据库构建；
[0063] S1.1、获取2005‑2009年华北平原(包含河北、河南、山东)冬小麦物候信息，如表1所示，一般情况下数据量越大，后期建模效果越好，本实施例选择了2005‑2009年；
[0064] S1.2、下载华北平原边界矢量图层；从2005‑2009年MODIS、SPOT 和MERIS卫星1km×1km分辨率的遥感影像数据中挑选覆盖华北平原区域且处于华北平原冬小麦关键生育期的遥感影像数据下载备用，如图2所示，同时结合当年当地的气象情况，进一步筛选云量较少，质量较好的遥感影像数据，运用ENVI软件对上述遥感影像进行拼接，使之成为涵盖整个华北平原的遥感影像数据，而后以华北平原边界矢量图层为掩膜裁剪该遥感影像数据，得到最终的华北平原冬小麦关键生育期各卫星遥感影像数据。由于华北平原区域较广，故本实施例采用了 2005‑2009年的1km×1km空间分辨率的SPOT‑VGT 7d、 ENVISAT MERIS 10d和MODIS 16d的卫星遥感影像数据。
[0065] S1.3、获取2005‑2009年华北平原(河北、河南、山东)冬小麦产量数据(施肥区)、冬小麦无肥区产量数据和对应的土壤数据；所述土壤数据包含土壤全氮、有机质、有效磷、有效钾和pH。
[0066] 表1 2005‑2009年华北平原(河北、河南、山东)三省冬小麦物候期时间表[0067]
[0068] (备注：表1所记录的物候期为大概日期，后期参考表1用于卫星遥感影像数据产品的选择时无需严格依据该时间段，若无物候期内的遥感数据，日期接近即可)[0069] S2、基于卫星遥感影像的植被指数计算，如图3所示；
[0070] S2.1、运用ENVI软件分别提取步骤S1获得的华北平原冬小麦关键生育期的各卫星遥感影像中的光谱特征，获得光谱特征参数图层，所述光谱特征参数图层包括近红外波段反射率ρNIR参数图层、红光波段反射率ρR参数图层、绿光波段反射率ρG参数图层和蓝光波段反射率ρB参数图层；所述冬小麦关键生育期包括返青期、拔节期、抽穗期和灌浆期；
[0071] S2.2、根据步骤S2.1的各光谱特征参数图层，通过公式1～公式5 计算植被指数特征，获得植被指数特征图层，所述植被指数特征图层包括归一化植被指数NDVI图层、差值植被指数DVI图层、垂直植被指数PVI图层、土壤调节植被指数SAVI图层和植被覆盖度VFC图层；
[0072]
[0073] DVI＝ρNIR‑ρR 公式2
[0074]
[0075]
[0076]
[0077] 公式1～公式5中，NDVI为归一化植被指数；ρNIR为近红外波段反射率；ρR为红光波段反射率；DVI为差值植被指数；PVI为垂直植被指数；SR为红光波段土壤反射率；VR为红光波段植被反射率；SNIR为近红外波段土壤反射率；VNIR为近红外波段植被反射率；SAVI为土壤调节植被指数；L为土壤调节系数，取值范围0‑1，L＝0时表示植被覆盖率为零，L＝1时，表示土壤背景的影像为零，即植被覆盖度非常高； VFC为植被覆盖度；NDVIsoil为完全是裸土或无植被覆盖区域的NDVI 值；NDVIveg代表完全被植被所覆盖的NDVI值，即纯植被像元的NDVI 值；
[0078] S2.3、根据步骤S2.1的各光谱特征参数图层，任意选取两个不同特征时相的同一波段的光谱特征进行差值计算，获得光谱差值特征图层，所述光谱差值特征图层包括近红外波段反射率差值图层、红光波段反射率差值图层、绿光波段反射率差值图层和蓝光波段反射率差值图层；
[0079] 本实施例选取2008‑03‑06和2008‑06‑10两个特征时相，分别计算各波段的光谱差值特征，如表2所示，其中，
[0080]
[0081]
[0082]
[0083]
[0084] S2.4、根据步骤S2.2的各植被指数特征图层，选取与步骤S2.3相同的两个特征时相的同一植被指数特征进行差值计算，获得植被指数差值特征图层，所述植被指数差值特征图层包括归一化植被指数差值图层、差值植被指数差值图层、垂直植被指数差值图层、土壤调节植被指数差值图层和植被覆盖度差值图层；
[0085] 本实施例选取2008‑03‑06和2008‑06‑10两个特征时相，分别计算各植被指数差值特征，如表2所示，其中，
[0086] NDVI6_NDVI3＝NDVI2008‑06‑1‑NDVI2008‑03‑06
[0087] DVI6_DVI3＝DVI2008‑06‑1‑DVI2008‑03‑0
[0088] PVI6_PVI3＝PVI2008‑06‑10‑PVI2008‑03‑06
[0089] SAVI6_SAVI3＝SAVI2008‑06‑1‑SAVI2008‑03‑0
[0090] VFC6_VFC3＝VFC2008‑06‑10‑VFC2008‑03‑06
[0091] 表2实施例遥感影像特征信息列表
[0092]
[0093]
[0094] S3、提取华北平原冬小麦种植区域，如图3所示；
[0095] 运用ENVI软件的合并功能将步骤S2获得的所有图层合并为一个图层；并将该图层与已知的冬小麦样本作为输入集，运行随机森林算法得到得华北平原冬小麦种植区域，如图4所示；其中，所述冬小麦样本即确定种植冬小麦的样本的坐标信息；
[0096] S4、建立华北平原冬小麦的遥感估产模型；
[0097] 通过步骤S2获得的华北平原冬小麦关键生育期的各卫星遥感影像中的归一化植被指数NDVI与步骤S1获得的冬小麦产量数据分别构建线性模型、多项式模型、指数模型、对数模型和乘幂模型，并通过模型评价指标进行最优模型筛选，所述模型评价指标包括决定2 2
系数R、相对误差R和均方根误差RMSE，依据决定系数R越大且相对误差R和均方根误差RMSE越小的原则确定最优的华北平原冬小麦的遥感估产模型；根据最优的华北平原冬小麦的遥感估产模型确定预测冬小麦产量最佳的卫星遥感影像数据和最佳生育期；
[0098] 线性模型：y＝a×NDVI+b
[0099] 多项式模型：y＝a×(NDVI)2+b×NDVI+c
[0100] 指数模型：y＝a×eb×NDVI
[0101] 对数模型：y＝a×ln(NDVI)+b
[0102] 乘幂模型：y＝a×(NDVI)b
[0103] 其中，y代表冬小麦产量，NDVI代表冬小麦某一关键生育期的归一化植被指数值；a、b、c代表方程的回归参数，并且不同方程之间的参数值不同。
[0104]
[0105] 其中，R2为决定系数， Yi、ESS是回归平方和，TSS是残差平方和；Yi表示参与建模的每个NDVI样本对应的实际产量，表示所有样本的对应实际产量的平均值，表示根据选定模型计算得到的每个样本对应的预测产量；
[0106]
[0107] 其中，R为相对误差，Δ为绝对误差，L为实际真值；
[0108]
[0109] 其中，RMSE为均方根误差；n代表测量次数，Xobs，i代表真实值 (观察值)，Xmodel，i表示模型预测值；
[0110] 本发明实施例在获取冬小麦产量y和归一化植被指数NDVI数据后运用Excel建模，2
模型参数a、b、c可由Excel自动计算，无需在意取值。而后依据R越大且RMSE和R越小原则筛选出最优的华北平原冬小麦的遥感估产模型。本发明实施例的最优的华北平原冬小麦的遥感估产模型为线性模型，结论是华北平原冬小麦估产的最佳的卫星遥感影像数据是SPOT卫星遥感影像数据，最佳生育期是拔节期。
[0111] S5、华北平原土壤生产力等级的识别；
[0112] S5.1、根据步骤S1获得的冬小麦产量数据和冬小麦无肥区产量数据建立产量关系模型，如图5所示；
[0113] S5.2、根据步骤S1获得的冬小麦无肥区产量数据划分土壤基础生产力，并与土壤数据建立土壤生产力等级模型，如表3所示；
[0114] 表3不同土壤生产力等级对应的土壤性质
[0115]
[0116] S5.3、通过步骤S4确定的最优的华北平原冬小麦的遥感估产模型、步骤S5.1建立的产量关系模型和步骤S5.2建立的土壤生产力等级模型建立华北平原冬小麦关键生育期卫星遥感影像中的归一化植被指数 NDVI与土壤生产力等级的关系模型，如图6以及表4所示，进而建立土壤性质仿真；
[0117] S5.4、根据步骤S4确定的最佳的(SPOT)卫星遥感影像数据和最佳生育期(拔节期)在步骤S3确定的华北平原冬小麦种植区域内提取任一点的归一化植被指数NDVI，借助步骤S5.3获得的华北平原冬小麦关键生育期卫星遥感影像中的归一化植被指数NDVI与土壤生产力等级的关系模型确定该点的土壤生产力等级并仿真出该点的土壤性质。
[0118] 表4华北平原冬小麦拔节期NDVI与产量的关系
[0119]
[0120]
