CN106295232A - 一种基于灰色关联分析的测土配方施肥方法 - Google Patents

Abstract

本发明公开了一种基于灰色关联分析的测土配方施肥方法，所述基于灰色关联分析的测土配方施肥方法采用多因子肥料效应函数法对产量与施肥量进行估计，使用灰色关联分析模型对肥料效应函数进行选择；与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小。本发明既满足了精度的要求，也加强了肥效模型的通用性，无需再针对某一区域在进行实验，充分挖掘了田间试验数据与土壤养分数据之间的隐含联系，解决了长期以来多因子肥效模型区域性限制问题，大大降低了测土配方施肥所需的前期工作量，满足了农业生产实践的需求，为测土配方施肥技术提供了一种新的解决方案。

Description

一种基于灰色关联分析的测土配方施肥方法

技术领域

本发明属于测土配方施肥技术领域，尤其涉及一种基于灰色关联分析的测土配方施肥方法。

背景技术

测土配方施肥技术是以土壤测试和肥料田间试验为基础，根据作物需肥规律、土壤供肥性能和肥料效应，在合理施用有机肥料的基础上，提出氮、磷、钾及中、微量元素等肥料的施用数量、施肥时期和施用方法。这其中，肥料效应模型的选择对测土配方施肥的科学性有着重要的作用，其计算结果的精度及可靠性将直接影响到测土配方施肥技术是否能达到提高肥料利用率和减少用量，提高作物产量，改善农产品品质，节省劳力，节支增收的目的。

现有测土配方施肥技术是由土肥专家根据当地田间试验的结果结合当地的环境气候特点与土壤养分数据，因地制宜的研究与构建出养分丰缺指标体系来指导当地施肥配肥。这种方式需要在当地进行过较多次的田间试验，积累不同年度的资料数据，并且对土肥专家的水平有着一定的要求，而且由于地区的限制，导致基于这些田间试验建立的肥效模型无法做到地区间的通用。由于现有技术有一方面是根据土肥专家以往种植经验来确定肥效模型的，模型受制于专家水平与主观因素的影响，计算结果的精度有待商榷，且每个地区都有其各自的施肥模型而无法通用，如若对实验数据暂缺区域应用测土配方施肥技术，一方面直接套用其他地区的肥效模型会导致结果的准确程度很差，而另一方面再对该区域进行相应的田间试验将耗费大量的时间与人力，这显然是不必要的。基于田间实验数据所建立的多因子肥效模型可以在产量与施肥量估计时达到相当的精度，但由于其具有一定的地区局限性，其适用范围往往限于当地已进行过肥效实验的地块而无法运用在其他地区。

发明内容

本发明的目的在于提供一种基于灰色关联分析的测土配方施肥方法，旨在解决现有测土配方施肥技术存在无法做到地区间的通用，计算结果不准确，耗费大量时间与人力的问题。

本发明是这样实现的，一种基于灰色关联分析的测土配方施肥方法，所述基于灰色关联分析的测土配方施肥方法采用多因子肥料效应函数法对产量与施肥量进行估计，使用灰色关联分析模型对肥料效应函数进行选择；与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小；采用氮、磷、钾肥建立施肥量与产量的关系函数如下：

$\hat{y}=b_0+b_{1}N+b_{2}P+b_{3}K+b_{12}NP+b_{13}NK+b_{23}PK+b_{11}{N}^2+b_{22}{P}^2+b_{33}{K}^2;$

式中：b_i为系数，N、P、K分别为氮磷钾的施用量，为产量；

所述灰色关联模型为邓氏关联度模型，进行过相应实验已探明施肥量与产量关系的地区和未进行过相关实验的地区，要素则代表反映这两个区域环境特征的量化指标，要素分别有如下量化的因子组成：

系统1：{0.035,0.215,0.325,0.475,1.475,2.225,3.225}

系统2：{0.045,0.315,0.451,0.451,1.201,1.201,2.201}

系统3：{0.141,0.555,0.829,1.221,1.721,1.721,2.721}

系统1与系统2的关联度大小为0.815，系统1与系统3的关联度大小为0.741。

进一步，所述与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小具体包括:

(1)根据目标作物选择模型库中所有的对应模型；

(2)计算所有被选中比较序列中各序列值与参考序列相应序列值的绝对差Δ_i(k)；序列：Y＝{1,6.1,9.3,13.6,35,63.6,92.1}；X＝{1,7,10,10,26.7,26.7,48.9}；

绝对差有如下计算方式：

Δ＝|Y(k)-X(k)|＝{0,0.86,0.74,3.5,8.3,36.9,43.2}

式中：Δ为绝对差，k为元素在序列中的位置；

(3)找到所有绝对差中的最小绝对差a与最大绝对差b；绝对差的结果：{0,0.86,0.74,3.5,8.3,36.9,43.2}；最小绝对差为0，最大绝对差为43.2；

(4)按以下公式计算所有选中的比较序列的各序列值与参考序列相应序列值的关联数ξ_i(k)：

ξ_i(k)＝(a+0.5b)/(Δ_i(k)+0.5b)

(5)经过计算得到各比较序列与参考序列的关联数序列，计算各序列的平均值作为各比较序列与参考序列的关联度，按以下公式进行计算：

$r_i=\frac{1}{n}\Σ{\mathrm{\ξ}}_i\left(k\right),(i=1,2,...,n);$

式中：r_i为参考序列与比较序列的关联度，n为序列中的元素数，Σξ_i(k)表示第i个序列的关联数之和；两个比较序列的关联数：

ξ₁＝{1,0.979，0.982,0.918,0.828,0.520,0.480}

ξ₂＝{1,0.947,0.921,0.890,0.636,0.440,0.354}

那么根据公式得：

$r_1=\frac{1}{7}(1+...+0.48)=0.815$

$r_2=\frac{1}{7}(1+...+0.354)=\mathrm{0.741.}$

进一步，所述基于灰色关联分析的测土配方施肥方法包括以下步骤：

步骤一，利用田间实验数据建立效应数学模型库，模型库中除模型的各项系数外，还包括作物名称、前作、前作产量及施肥量、海拔、坡度、气候区、土壤类型、土壤pH值以及土壤中的有机质、碱解氮、速效钾、有效磷的含量，并数据进行标准化处理；

步骤二，利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量；

步骤三，将模型库中标准化后的数据累加，生成各模型比较序列；

步骤四，在对目标地块进行产量和施肥量的计算时，通过调用土壤信息数据库和直接输入的方法获得地块的前作情况、土壤养分状况、酸碱度、地块基本情况，将数据进行量化和标准化构成参考序列，与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小；调用其中关联度最大的模型所对应的最大产量及施肥量作为该地块某作物的施肥方案；

步骤五，调用上述模型进行最经济产量的计算，即按照报酬递减规律，根据当前肥料价格及作物单价，利用边际产量等于边际产值时利润最大的原理，求施肥模型中三个因子的一阶偏导数等于边际成本时的解作为最经济施肥量，对应产量为最经济产量，以此作为另一施肥方案。

进一步，利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量包括：

对于两因素试验配置的二元一次方程组，利用边际产量为0时的解求解函数极值点作为最大产量和施肥量；对于三因素及以上的模型，采用数值迭代的方法，给定一个合理的迭代区间求解模型函数在该区间的极大值点，将其作为最大产量与施肥量保存至模型库中；

二元二次肥料效应回归方程式有如下形式：

y＝B₀+B₁X₁+B₂X₁ ²+B₃X₂+B₄X₂ ²+B₅X₁X₂；

式中：y为产量，B_i为系数，X₁、X₂分别为两种肥料的用量，根据方程可求出产量y对施肥量X₁、X₂的偏导数即边际产量：

$\frac{\∂y}{\∂X_1}=B_1+2B_2{X}_1+B_5{X}_2$

$\frac{\∂y}{\∂X_2}=B_3+2B_4{X}_2+B_5{X}_1$

又当且时，可知该效应函数的拟合曲面为凸形，函数一定有极大值点，并满足时，对应的施肥量为最高产量施肥量，也就是边际产量均为0时，获得最大的产量。

进一步，对于三因素及其以上的方程，求解其在某一区间的极大值利用求约束条件下n维极值的复形调优法目标函数为：

J＝-f(x₀+x₁+x₂)；

式中：J为所要求解最大产量的相反数，f(x₀+x₁+x₂)为模型库中的多因子肥效函数，x_i分别为氮、磷、钾三种肥料的施用量；

常量约束条件为：

a_i＜x_i＜b_i；

式中：a_i为多因子肥效函数对应实验数据的0水平施肥量，b_i为多因子肥效函数对应实验数据的3水平施肥量；

函数约束条件为：

0＜f(x₀+x₁+x₂)；

由约束条件利用复形调优法求解J的极小值即f(x₀+x₁+x₂)的极大值的过程如下所示

复形共有2n个顶点，设给定初始复形中的第一个顶点坐标：

X₍₀₎＝(x₀₀,x₁₀,···,x_n-1,0)；

且此顶点坐标满足所有的常数约束条件和函数约束条件；

(1)在n维变量空间中在确定出初始复形的其余2n-1个顶点，其方法如下：利用伪随机数按常量约束条件产生第j个顶点X_(j)＝(x_0j,x_1j,···,x_n-1,j)(j＝1,2,···,2n-1)中的各分量x_ij(i＝1,2,···,2n-1)，即

x_ij＝a_i+r(b_i-a_i)；

式中：为r是区间[0,1]之间的一个伪随机数；

在检查是否符合函数约束条件，如果不符合，则需要作调整，直到全部顶点均符合常量约束和函数约束条件为止；调整的原则为：

前j个顶点以满足所有的约束条件，而第j+1个顶点不满足约束条件，则做如下调整变换(j＝1,2,···,2n-1)：

X_(j+1)＝(X_(j+1)+T)/2；

其中：

$T=\frac{1}{j}\underset{k=1}{\overset{j}{\Σ}}{X}_{\left(k\right)};$

初始复形的2n个顶点确定以后，计算各顶点处的目标函数值：

J_(j)＝-f(X_(j)),j＝0,···,2n-1

(2)确定：

$J_{\left(R\right)}=-f\left(X_{\left(R\right)}\right)=\underset{0\≤i\≤2n-1}{max}(-f(i\left)\right);$

$J_{\left(G\right)}=-f\left(X_{\left(G\right)}\right)=\underset{0\≤i\≤2n-1;i\≠R}{max}(-f(i\left)\right);$

其中：X_(R)被称为最坏点；

(3)计算最坏点的对称点

X_T＝(1+α)X_F-αX_(R)；

式中：

$X_F=\frac{1}{2n-1}\underset{i\≠R}{\overset{2n-1}{\underset{i=0}{\Σ}}}{X}_{\left(i\right)};$

α称为反射系数，取1.3；

(4)确定一个新的顶点替代最坏点X_(R)以构成新的复形，其方法如下：

如果J(X_T)＞J(X_(G))，则用下式修改X_T：

X_T＝(X_F+X_T)/2；

直到J(X_T)≤J(X_(G))为止；

然后检查X_T是否满足所有约束条件，如果对于某个分量X_T(j)不满足常量约束条件，即如果X_T(j)＜a_j或者X_T(j)＞b_j；

则令：

X_T(j)＝a_j+δ或X_T(j)＝b_j-δ；

式中：δ在本发明中取10-⁶，重复步骤(4)；

如果X_T不满足函数约束条件，则用下式修改X_T：

X_T＝(X_F+X_T)/2；

重复(4)；

直到-f(X_T)≤-f(X_(G))且满足所有约束条件为止，令：

X_(R)＝X_T,f(X_(R))＝f(X_T)；

重复(2)～(4)，直到复形中各顶点的距离小于预先给定的精度要求为止，也就代表迭代满足了原先设定的精度要求，搜索到了极值点。

进一步，标准化是将模型库中除拟合函数系数外，其余数据均通过极值化或均值化以消除各类型数据不同量纲的影响，累加是将标准化后的数据依次相加。

进一步，标准化是通过对所收集的数据进行处理，通过各种数据变换消除其量纲；

对数据序列X＝(x(1),x(2),···,x(n))变换得到Y＝(y(1),y(2),···,y(n))，其中

$y\left(k\right)=\frac{x\left(k\right)}{\underset{k}{\max x\left(k\right)}},k=1,2,...,n;$

则称由序列X到序列Y的变换为极值化处理；

变换为以下形式：

$y\left(k\right)=\frac{x\left(k\right)}{\stackrel{\‾}{X}},k=1,2,\mathrm{...},n;\stackrel{\‾}{X}=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}x\left(k\right);$

X序列到Y序列的变换为均值化处理；

土壤中碱解氮含量所构成的序列如下所示：

{171，160，97，170，290}；

则经过极值化处理后转化为序列：

{0.590,0.552,0.334,0.586,1}

经过均值化处理后转化为序列:

{0.963,0.901,0.546,0.957,1.633}。

本发明提供的基于灰色关联分析的测土配方施肥方法，采用多因子肥料效应函数法对产量与施肥量进行估计，使用灰色关联分析模型对肥效函数进行科学的选择，既满足了精度的要求，也加强了肥效模型的通用性，无需再针对某一区域在进行实验，充分挖掘了田间试验数据与土壤养分数据之间的隐含联系，解决了长期以来多因子肥效模型区域性限制问题，大大降低了测土配方施肥所需的前期工作量，满足了农业生产实践的需求。应用灰色关联分析对实验地块所处的环境进行定量描述，确定各肥效模型的适用条件，在实际运用中只需要获得目标地块相应的数据，即可关联到与之最相似的实验地块并调用模型进行产量估计，解决了多因子肥料效应模型的区域性限制问题，为测土配方施肥技术提供了一种新的解决方案。

本发明将田间试验数据与土壤养分调查数据进行了有效整合，拓展了其用途；采用在农业上用途较为广泛的灰色关联分析模型进行环境特征的识别，有效判断了施肥模型的适用条件；采用数值迭代而不是求导的方式求解多因子肥效模型的极值，保证了解的准确性和合理性。本发明计算精度较高，得出的结果与生产实践贴近；支持作物种类较多；适用性较强，对于任意区域，只要其与模型所对应区域的关联度达到一定阈值，即可从模型库中调用模型进行计算，无需再针对某一区域进行相应的肥效实验；可扩展性强，随着实验数据的不断增加，本发明求得结果的精确程度以及支持作物种类的数量也都会随之增加。

附图说明

图1是本发明实施例提供的基于灰色关联分析的测土配方施肥方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明采用多因子肥料效应函数法对产量与施肥量进行估计，使用灰色关联分析模型对肥效函数进行科学的选择，既满足了精度的要求，也加强了肥效模型的通用性，无需再针对某一区域在进行实验，充分挖掘了田间试验数据与土壤养分数据之间的隐含联系，解决了长期以来多因子肥效模型区域性限制问题，大大降低了测土配方施肥所需的前期工作量，满足了农业生产实践的需求。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例的基于灰色关联分析的测土配方施肥方法包括以下步骤：

S101：利用田间实验数据建立效应数学模型库；

S102：利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量；

S103：将模型库中标准化后的数据累加，生成各模型比较序列；

S104：在对目标地块进行产量和施肥量的计算时，通过调用土壤信息数据库和直接输入的方法获得地块的前作情况、土壤养分状况、酸碱度、地块基本情况，将这些数据进行量化和标准化构成参考序列，与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小；

S105：调用模型进行最经济产量的计算，即按照报酬递减规律，根据当前肥料价格及作物单价，利用边际产量等于边际产值时利润最大的原理，求施肥模型中三个因子的一阶偏导数等于边际成本时的解作为最经济施肥量，对应产量为最经济产量。

在本发明的实施例中：

1、多因子肥料效应函数法是通过肥力效应试验数据建立的能够体现施肥量与产量之间数量关系的函数，多因子是指函数中包含多种肥料及其组合，在本发明中，采用氮、磷、钾肥建立施肥量与产量的关系函数，函数形式如下所示：

$\hat{y}=b_0+b_{1}N+b_{2}P+b_{3}K+b_{12}NP+b_{13}NK+b_{23}PK+b_{11}{N}^2+b_{22}{P}^2+b_{33}{K}^2;$

式中：b_i为系数，N、P、K分别为氮磷钾的施用量(纯质，千克/亩)，为产量(千克/亩)。

2、本发明所利用的灰色关联模型为邓氏关联度模型，它的意义在于，对于两系统之间的要素(在本发明中两系统分别表示进行过相应实验已探明施肥量与产量关系的地区和未进行过相关实验的地区，要素则代表反映这两个区域环境特征的量化指标，如土壤中有机质的含量等)，其随时间或不同对象而变化的关联度大小的量度，对于一个系统发展变化态势提供了量化的量度(在本发明中关联度的大小就是两地块相似的程度)，例如有3个系统，其要素分别有如下量化的因子组成：

系统1：{0.035,0.215,0.325,0.475,1.475,2.225,3.225}

系统2：{0.045,0.315,0.451,0.451,1.201,1.201,2.201}

系统3：{0.141,0.555,0.829,1.221,1.721,1.721,2.721}

如何量化的判断系统1与其余系统的关联度(元素的相似性及系统的发展态势)就是灰色关联度要解决的主要问题，经过计算可知系统1与系统2的关联度大小为0.815，系统1与系统3的关联度大小为0.741，由此可以判断系统1与系统2的因子及发展态势更为相似。

3、在使用本发明进行产量与施肥量计算时，应首先获知所要预测作物的名称，例如对某一区域的马铃薯产量及施肥量进行计算，那么目标作物即为马铃薯。模型库即为存储作物施肥模型的数据库，其中每个模型都对应了一个特定的作物，在得知了目标作物以后，就需要在模型库中选择出所有的该作物对应的模型进行下一步的计算。

4、例如有2系统序列如下所示：

Y＝{1,6.1,9.3,13.6,35,63.6,92.1}

X＝{1,7,10,10,26.7,26.7,48.9}

那么这两个系统的绝对差有如下计算方式：

Δ＝|Y(k)-X(k)|＝{0,0.86,0.74,3.5,8.3,36.9,43.2}

式中：Δ为绝对差，k为元素在序列中的位置。

在具体运用时这一步通过电脑编程来实现。

5、依旧以上述两系统为例，由绝对差的结果：

{0,0.86,0.74,3.5,8.3,36.9,43.2}

可以看出该例中最小绝对差为0，最大绝对差为43.2，对于由多个绝对差序列组成的数据，则找出在所有绝对差中最小差和最大差。具体运用中通过电脑编程来实现。

6、按以下公式进行计算：

式中：r_i为参考序列与比较序列的关联度，n为序列中的元素数，∑ξ_i(k)表示第i个序列的关联数之和。

例如有两个比较序列的关联数：

ξ₁＝{1,0.979，0.982,0.918,0.828,0.520,0.480}

ξ₂＝{1,0.947,0.921,0.890,0.636,0.440,0.354}

那么根据公式得：

$r_1=\frac{1}{7}(1+...+0.48)=0.815$

$r_2=\frac{1}{7}(1+...+0.354)=0.741$

7、二元二次肥料效应回归方程式有如下形式：

y＝B₀+B₁X₁+B₂X₁ ²+B₃X₂+B₄X₂ ²+B₅X₁X₂；

式中：y为产量，B_i为系数，X₁、X₂分别为两种肥料的用量。那么根据此方程可求出产量y对施肥量X₁、X₂的偏导数即边际产量：

$\frac{\∂y}{\∂X_1}=B_1+2B_2{X}_1+B_5{X}_2$

$\frac{\∂y}{\∂X_2}=B_3+2B_4{X}_2+B_5{X}_1$

又当且时，可知该效应函数的拟合曲面为凸形，函数一定有极大值点，并满足时，对应的施肥量为最高产量施肥量，也就是边际产量均为0时，可以获得最大的产量，这就是对于二元肥效方程组的处理过程。

8、对于三因素及其以上的方程，本发明求解其在某一区间的极大值主要是利用求约束条件下n维极值的复形调优法，在本发明中，目标函数为：

J＝-f(x₀+x₁+x₂)；

式中：J为所要求解最大产量的相反数，f(x₀+x₁+x₂)为模型库中的多因子肥效函数，x_i分别为氮、磷、钾三种肥料的施用量。

常量约束条件为：

a_i＜x_i＜b_i；

式中：a_i为多因子肥效函数对应实验数据的0水平施肥量，b_i为多因子肥效函数对应实验数据的3水平施肥量。

函数约束条件为：

0＜f(x₀+x₁+x₂)；

由这些约束条件利用复形调优法求解J的极小值即f(x₀+x₁+x₂)的极大值的过程如下所示：

复形共有2n(在本发明中n为3)个顶点，设给定初始复形中的第一个顶点坐标：

X₍₀₎＝(x₀₀,x₁₀,···,x_n-1,0)；

且此顶点坐标满足所有的常数约束条件和函数约束条件。

(1)在n维变量空间中在确定出初始复形的其余2n-1个顶点，其方法如下：利用伪随机数按常量约束条件产生第j个顶点X_(j)＝(x_0j,x_1j,···,x_n-1,j)(j＝1,2,···,2n-1)中的各分量x_ij(i＝1,2,···,2n-1)，即

x_ij＝a_i+r(b_i-a_i)；

式中：为r是区间[0,1]之间的一个伪随机数。

显然，由上述方法产生的初始复形的各顶点满足常量约束条件。然后在检查它们是否符合函数约束条件，如果不符合，则需要作调整，直到全部顶点均符合常量约束和函数约束条件为止。调整的原则为：

假设前j个顶点以满足所有的约束条件，而第j+1个顶点不满足约束条件，则做如下调整变换(j＝1,2,···,2n-1)：

X_(j+1)＝(X_(j+1)+T)/2；

其中：

$T=\frac{1}{j}\underset{k=1}{\overset{j}{\Σ}}{X}_{\left(k\right)}$

这个过程一直做到满足所有约束条件为止。

初始复形的2n个顶点确定以后，计算各顶点处的目标函数值：

J_(j)＝-f(X_(j)),j＝0,···,2n-1

(2)确定：

$J_{\left(R\right)}=-f\left(X_{\left(R\right)}\right)=\underset{0\≤i\≤2n-1}{max}(-f(i\left)\right)$

$J_{\left(G\right)}=-f\left(X_{\left(G\right)}\right)=\underset{0\≤i\≤2n-1;i\≠R}{max}(-f(i\left)\right)$

其中：X_(R)被称为最坏点。

(3)计算最坏点的对称点

X_T＝(1+α)X_F-αX_(R)

式中：

$X_F=\frac{1}{2n-1}\underset{i\≠R}{\overset{2n-1}{\underset{i=0}{\Σ}}}{X}_{\left(i\right)}$

α称为反射系数，在本发明中取1.3。

(4)确定一个新的顶点替代最坏点X_(R)以构成新的复形，其方法如下：

如果J(X_T)＞J(X_(G))，则用下式修改X_T：

X_T＝(X_F+X_T)/2

直到J(X_T)≤J(X_(G))为止。

然后检查X_T是否满足所有约束条件，如果对于某个分量X_T(j)不满足常量约束条件,即如果X_T(j)＜a_j或者X_T(j)＞b_j

则令：

X_T(j)＝a_j+δ或X_T(j)＝b_j-δ；

式中：δ在本发明中取10-⁶。重复步骤(4)。

如果X_T不满足函数约束条件，则用下式修改X_T：

X_T＝(X_F+X_T)/2

重复(4)。

直到-f(X_T)≤-f(X_(G))且满足所有约束条件为止。此时令：

X_(R)＝X_T,f(X_(R))＝f(X_T)

重复(2)～(4)，直到复形中各顶点的距离小于预先给定的精度要求为止，也就代表迭代满足了原先设定的精度要求，搜索到了极值点。

在实际运用中，以上步骤主要依靠计算机来完成。

9、标准化是通过对所收集的数据进行处理，通过各种数据变换消除其量纲，使其具有可比性，以保证建模质量和系统分析的正确结果。其中极值化处理和均值化处理是常见的标准化处理手段。

对数据序列X＝(x(1),x(2),···,x(n))变换得到Y＝(y(1),y(2),···,y(n))，其中

$y\left(k\right)=\frac{x\left(k\right)}{\underset{k}{\max x\left(k\right)}},k=1,2,...,n$

则称由序列X到序列Y的变换为极值化处理。

又若以上的变换为以下形式：

$y\left(k\right)=\frac{x\left(k\right)}{\stackrel{\‾}{X}},k=1,2,\mathrm{...},n;\stackrel{\‾}{X}=\frac{1}{N}\underset{k=1}{\overset{n}{\Σ}}x\left(k\right)$

那么就称X序列到Y序列的变换为均值化处理。

例如本发明中由土壤中碱解氮含量所构成的序列如下所示：

{171，160，97，170，290}

则经过极值化处理后转化为序列：

{0.590,0.552,0.334,0.586,1}

经过均值化处理后转化为序列:

{0.963,0.901,0.546,0.957,1.633}

下面结合具体实施例对本发明的应用原理作详细的描述。

1、利用田间实验数据建立效应数学模型库。模型库中除模型的各项系数外，还包括作物名称、前作、前作产量及施肥量、海拔、坡度、气候区、土壤类型、土壤pH值以及土壤中的有机质、碱解氮、速效钾、有效磷的含量，并将这些数据进行标准化处理。

2、利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量。采用数值迭代的方法，给定一个合理的迭代区间求解模型函数在该区间的极大值点，将其作为最大产量与施肥量保存至模型库中。迭代区间应选取作物2水平施肥量的邻近区间，以此保证迭代结果的合理性。

3、将模型库中标准化后的数据累加，生成各模型比较序列。标准化是将模型库中除拟合函数系数外，其余数据均通过极值化(除以该类型数据的极大值)或均值化(除以该类型数据的平均值)以消除各类型数据不同量纲的影响，累加是将标准化后的数据依次相加，以加强数据变化的规律性。

4、在对目标地块进行产量和施肥量的计算时，通过调用土壤信息数据库和直接输入的方法获得地块的前作情况、土壤养分状况、酸碱度、地块基本情况，将这些数据进行量化和标准化构成参考序列，与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小，具体方法是:

(1)根据目标作物选择模型库中所有的对应模型；

(2)计算所有被选中比较序列中各序列值与参考序列相应序列值的绝对差Δ_i(k)；

(3)找到所有绝对差中的最小绝对差a与最大绝对差b；

(4)按以下公式计算所有选中的比较序列的各序列值与参考序列相应序列值的关联数ξ_i(k)：

ξ_i(k)＝(a+0.5b)/(Δ_i(k)+0.5b)；

(5)经过上步的计算可得到各比较序列与参考序列的关联数序列，计算各序列的平均值作为各比较序列与参考序列的关联度。

调用其中关联度最大的模型所对应的最大产量及施肥量作为该地块某作物的推荐施肥方案一。

5、调用上述模型进行最经济产量的计算，即按照报酬递减规律，根据当前肥料价格及作物单价，利用边际产量等于边际产值时利润最大的原理，求施肥模型中三个因子的一阶偏导数等于边际成本时的解作为最经济施肥量，对应产量为最经济产量，以此作为推荐施肥方案二。

下面结合实验对本发明的应用效果作详细的描述。

表1是应用本发明提供的测土配方施肥技术的某实验田与该地区常规施肥的试验田的对比(模型库中无此区域的实验数据与施肥模型)：

由该表可以看出，基于灰色关联模型的测土配方施肥技术即本发明对于作物有着明显增产效果，同时其预报产量与实际产量的误差较小，实验证明，该技术可以应用于实际的农业生产中。

另一方面，使用现有技术对该实验田进行产量预测，则要求在该地区进行

对此实验，累计多年的实验数据，由专家根据实验结果和种植经验建立施肥指标体系后，才可以进行作物产量及施肥量的预测。

应用现有技术对以上试验田进行预测，由于缺少该地区的施肥指标体系，只能调用与之相近区域的施肥模型进行估计，按照上述实验最高产量509.33千克/亩使用现有技术对其施肥量进行估计，得应施纯N 11.89千克/亩、纯P₂O₅1.6千克/亩、纯K₂O 11.5千克，对比上表结果发现由现有技术预测的施肥量与实际施肥量存在较为明显的误差。

造成这种较大误差的原因是因为该地区缺少相应的实验数据，无法因地制宜为其配置施肥指标体系，而使用其他地区的施肥模型则会造成预测结果的失真。虽然该地配置肥效实验可以解决该问题，但需要在当地配置多点试验，累积多不同年度的数据资料，不利于测土配方施肥技术的快速推广。

根据模型库中所收录到的实验数据来看，本发明目前支持的作物品种有：白菜，茶叶，大豆，大麦，甘蓝，甘蔗，辣椒，马铃薯，荞麦，水稻，莴笋，西兰花，小麦，油菜，玉米。随着实验数据的逐步补充，本发明所支持的作物种类将进一步增加。

具体的应用实施例如下所示：

假设所要估计产量的地块的各项数据如下表所示：

又已知在该地块上进行过的田间试验分析显示该地块上的水稻最高产量为720千克/亩，且模型库中没有这条试验数据。对该地块应用本发明提供的技术，得到的结果如下表所示：

由此可以看出即使在模型库中没有该实验数据的前提下，使用本发明会自动关联到生产能力相似的地块对产量及相应施肥量进行估计，从而得到一个较为合理的结果。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于灰色关联分析的测土配方施肥方法，其特征在于，所述基于灰色关联分析的测土配方施肥方法采用多因子肥料效应函数法对产量与施肥量进行估计，使用灰色关联分析模型对肥料效应函数进行选择；与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小；采用氮、磷、钾肥建立施肥量与产量的关系函数如下：

$\hat{y}=b_0+b_{1}N+b_{2}P+b_{3}K+b_{12}NP+b_{13}NK+b_{23}PK+b_{11}{N}^2+b_{22}{P}^2+b_{33}{K}^2;$

式中：b_i为系数，N、P、K分别为氮磷钾的施用量，为产量；

所述灰色关联模型为邓氏关联度模型，进行过相应实验已探明施肥量与产量关系的地区和未进行过相关实验的地区，其要素代表反映这两个区域环境特征的量化指标，要素分别有如下量化的因子组成：

系统1：{0.035,0.215,0.325,0.475,1.475,2.225,3.225}

系统2：{0.045,0.315,0.451,0.451,1.201,1.201,2.201}

系统3：{0.141,0.555,0.829,1.221,1.721,1.721,2.721}

系统1与系统2的关联度大小为0.815，系统1与系统3的关联度大小为0.741。

2.如权利要求1所述的基于灰色关联分析的测土配方施肥方法，其特征在于，所述与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小具体包括:

(1)根据目标作物选择模型库中所有的对应模型；

(2)计算所有被选中比较序列中各序列值与参考序列相应序列值的绝对差Δ_i(k)；序列：Y＝{1,6.1,9.3,13.6,35,63.6,92.1}；X＝{1,7,10,10,26.7,26.7,48.9}；

绝对差有如下计算方式：

Δ＝|Y(k)-X(k)|＝{0,0.86,0.74,3.5,8.3,36.9,43.2}

式中：Δ为绝对差，k为元素在序列中的位置；

(3)找到所有绝对差中的最小绝对差a与最大绝对差b；绝对差的结果：{0,0.86,0.74,3.5,8.3,36.9,43.2}；最小绝对差为0，最大绝对差为43.2；

(4)按以下公式计算所有选中的比较序列的各序列值与参考序列相应序列值的关联数ξ_i(k)：

ξ_i(k)＝(a+0.5b)/(Δ_i(k)+0.5b)

(5)经过计算得到各比较序列与参考序列的关联数序列，计算各序列的平均值作为各比较序列与参考序列的关联度，按以下公式进行计算：

$r_i=\frac{1}{n}{\mathrm{\Σ\ξ}}_i\left(k\right),(i=1,2,...,n);$

式中：r_i为参考序列与比较序列的关联度，n为序列中的元素数，∑ξ_i(k)表示第i个序列的关联数之和；两个比较序列的关联数：

ξ₁＝{1,0.979，0.982,0.918,0.828,0.520,0.480}

ξ₂＝{1,0.947,0.921,0.890,0.636,0.440,0.354}

那么根据公式得：

$r_1=\frac{1}{7}(1+...+0.48)=0.815$

$r_2=\frac{1}{7}(1+...+0.354)=\mathrm{0.741.}$

3.如权利要求1所述的基于灰色关联分析的测土配方施肥方法，其特征在于，所述基于灰色关联分析的测土配方施肥方法包括以下步骤：

步骤一，利用田间实验数据建立效应数学模型库，模型库中除模型的各项系数外，还包括作物名称、前作、前作产量及施肥量、海拔、坡度、气候区、土壤类型、土壤pH值以及土壤中的有机质、碱解氮、速效钾、有效磷的含量，并数据进行标准化处理；

步骤二，利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量；

步骤三，将模型库中标准化后的数据累加，生成各模型比较序列；

步骤四，在对目标地块进行产量和施肥量的计算时，通过调用土壤信息数据库和直接输入的方法获得地块的前作情况、土壤养分状况、酸碱度、地块基本情况，将数据进行量化和标准化构成参考序列，与模型库中相应数据所构成比较序列进行几何相似性的计算即计算其与各比较序列的灰色关联度大小；调用其中关联度最大的模型所对应的最大产量及施肥量作为该地块某作物的施肥方案；

步骤五，调用上述模型进行最经济产量的计算，即按照报酬递减规律，根据当前肥料价格及作物单价，利用边际产量等于边际产值时利润最大的原理，求施肥模型中三个因子的一阶偏导数等于边际成本时的解作为最经济施肥量，对应产量为最经济产量。

4.如权利要求3所述的基于灰色关联分析的测土配方施肥方法，其特征在于，利用数值迭代逼近的方法求取模型库中各模型的区间极值与对应的施肥量包括：

对于两因素试验配置的二元一次方程组，利用边际产量为0时的解求解函数极值点作为最大产量和施肥量；对于三因素及以上的模型，采用数值迭代的方法，给定一个合理的迭代区间求解模型函数在该区间的极大值点，将其作为最大产量与施肥量保存至模型库中；

二元二次肥料效应回归方程式有如下形式：

y＝B₀+B₁X₁+B₂X₁ ²+B₃X₂+B₄X₂ ²+B₅X₁X₂；

式中：y为产量，B_i为系数，X₁、X₂分别为两种肥料的用量，根据方程可求出产量y对施肥量X₁、X₂的偏导数即边际产量：

$\frac{\∂y}{\∂X_1}=B_1+2B_2{X}_1+B_5{X}_2$

$\frac{\∂y}{\∂X_2}=B_3+2B_4{X}_2+B_5{X}_1$

又当且时，可知该效应函数的拟合曲面为凸形，函数一定有极大值点，并满足时，对应的施肥量为最高产量施肥量，也就是边际产量均为0时，获得最大的产量。

5.如权利要求4所述的基于灰色关联分析的测土配方施肥方法，其特征在于，对于三因素及其以上的方程，求解其在某一区间的极大值利用求约束条件下n维极值的复形调优法目标函数为：

J＝-f(x₀+x₁+x₂)；

式中：J为所要求解最大产量的相反数，f(x₀+x₁+x₂)为模型库中的多因子肥效函数，x_i分别为氮、磷、钾三种肥料的施用量；

常量约束条件为：

a_i＜x_i＜b_i；

式中：a_i为多因子肥效函数对应实验数据的0水平施肥量，b_i为多因子肥效函数对应实验数据的3水平施肥量；

函数约束条件为：

0＜f(x₀+x₁+x₂)；

由约束条件利用复形调优法求解J的极小值即f(x₀+x₁+x₂)的极大值的过程如下所示

复形共有2n个顶点，设给定初始复形中的第一个顶点坐标：

X₍₀₎＝(x₀₀,x₁₀,…,x_n-1,0)；

且此顶点坐标满足所有的常数约束条件和函数约束条件；

(1)在n维变量空间中在确定出初始复形的其余2n-1个顶点，其方法如下：利用伪随机数按常量约束条件产生第j个顶点X_(j)＝(x_0j,x_1j,…,x_n-1,j)(j＝1,2,…,2n-1)中的各分量x_ij(i＝1,2,…,2n-1)，即

x_ij＝a_i+r(b_i-a_i)；

式中：为r是区间[0,1]之间的一个伪随机数；

在检查是否符合函数约束条件，如果不符合，则需要作调整，直到全部顶点均符合常量约束和函数约束条件为止；调整的原则为：

前j个顶点以满足所有的约束条件，而第j+1个顶点不满足约束条件，则做如下调整变换(j＝1,2,…,2n-1)：

X_(j+1)＝(X_(j+1)+T)/2；

其中：

$T=\frac{1}{j}\underset{k=1}{\overset{j}{\Σ}}{X}_{\left(k\right)};$

初始复形的2n个顶点确定以后，计算各顶点处的目标函数值：

J_(j)＝-f(X_(j)),j＝0,…,2n-1

(2)确定：

$J_{\left(R\right)}=-f\left(X_{\left(R\right)}\right)=\underset{0\≤i\≤2n-1}{max}(-f(i\left)\right);$

$J_{\left(G\right)}=-f\left(X_{\left(G\right)}\right)=\underset{0\≤i\≤2n-1;i\≠R}{max}(-f(i\left)\right);$

其中：X_(R)被称为最坏点；

(3)计算最坏点的对称点

X_T＝(1+α)X_F-αX_(R)；

式中：

$X_F=\frac{1}{2n-1}\underset{i\≠R}{\overset{2n-1}{\underset{i=0}{\mathrm{\Σ}}}}{X}_{\left(i\right)};$

α称为反射系数，取1.3；

(4)确定一个新的顶点替代最坏点X_(R)以构成新的复形，其方法如下：

如果J(X_T)＞J(X_(G))，则用下式修改X_T：

X_T＝(X_F+X_T)/2；

直到J(X_T)≤J(X_(G))为止；

然后检查X_T是否满足所有约束条件，如果对于某个分量X_T(j)不满足常量约束条件，即如果X_T(j)＜a_j或者X_T(j)＞b_j；

则令：

X_T(j)＝a_j+δ或X_T(j)＝b_j-δ；

式中：δ在本发明中取10^-6，重复步骤(4)；

如果X_T不满足函数约束条件，则用下式修改X_T：

X_T＝(X_F+X_T)/2；

重复(4)；

直到-f(X_T)≤-f(X_(G))且满足所有约束条件为止，令：

X_(R)＝X_T,f(X_(R))＝f(X_T)；

重复(2)～(4)，直到复形中各顶点的距离小于预先给定的精度要求为止，也就代表迭代满足了原先设定的精度要求，搜索到了极值点。

6.如权利要求3所述的基于灰色关联分析的测土配方施肥方法，其特征在于，标准化是将模型库中除拟合函数系数外，其余数据均通过极值化或均值化以消除各类型数据不同量纲的影响，累加是将标准化后的数据依次相加。

7.如权利要求6所述的基于灰色关联分析的测土配方施肥方法，其特征在于，标准化是通过对所收集的数据进行处理，通过各种数据变换消除其量纲；

对数据序列X＝(x(1),x(2),…,x(n))变换得到Y＝(y(1),y(2),…,y(n))，其中

$y\left(k\right)=\frac{x\left(k\right)}{\underset{k}{\max x\left(k\right)}},k=1,2,...,n;$

则称由序列X到序列Y的变换为极值化处理；

变换为以下形式：

$y\left(k\right)=\frac{x\left(k\right)}{\stackrel{\‾}{X}},k=1,2,...,n;\stackrel{\‾}{X}=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}x\left(k\right);$

X序列到Y序列的变换为均值化处理；

土壤中碱解氮含量所构成的序列如下所示：

{171，160，97，170，290}；

则经过极值化处理后转化为序列：

{0.590,0.552,0.334,0.586,1}

经过均值化处理后转化为序列:

{0.963,0.901,0.546,0.957,1.633}。