三元非结构肥效模型提高水稻施肥推荐的可靠性

李 娟，章明清*，章赞德，许文江，姚宝全

（1 福建省农业科学院土壤肥料研究所，福州 350013；2 福建省大田县农田建设与土壤肥料技术推广站，福建 366100；3 福建省亚热带植物研究所，福建厦门 361000；4 福建省农田建设与土壤肥料技术推广总站，福州 350003）

当前，通过构建氮磷钾三元肥效模型确定最佳施肥量，是实现水稻计量施肥的重要技术途径之一[1-2]，其中，二次多项式函数是研究和应用最多的肥效模型种类[3-7]。但是，众多研究表明，一元和二元二次多项式肥效模型典型式的比例分别仅占60%左右和40.2%[8-9]，三元二次多项式肥效模型的典型式比例则更低至23.6%[10]。由于构建肥效模型时出现了大量非典型式，严重削弱了该法的计量精确性和实用价值。尽管国内外学者对此进行了深入研究，提出了许多有意义的改进措施[3]，但该问题至今仍然困扰着计量施肥研究和应用。纵观这些研究和改进措施中，对二次多项式等函数形式的肥效模型本身是否存在模型设定偏误及提出改进建议还鲜见研究报道。

研究表明，当前广泛应用的一元、二元、三元二次多项式肥效模型及其它类似的多项式肥效模型存在明显的设定偏误[3]。为此，章明清等[11]提出了一元非结构肥效模型，较好地克服了模型设定缺陷。与一元二次多项式肥效模型相比，新模型在拟合水稻氮磷钾单因素田间肥效试验结果时，具有较高的拟合精度、较宽的适用范围和推荐施肥量较低等优点。为此，在一元非结构肥效模型基础上，本研究利用氮磷钾三因素田间肥效试验结果，探讨三元非结构肥效模型的构建方法及其对田间肥效试验资料的拟合效果，旨在扩大三元肥效模型的适用性，为计量施肥研究和应用提供一种新的模型方法。

1 材料与方法

1.1 水稻氮磷钾田间肥效试验设计

近10年来，笔者参加福建省莆田市仙游县和漳州市平和县的水稻测土配方施肥工作。两个项目县均地处南亚热带海洋性季风气候带，气候条件尤其适合水稻生长发育。 稻田土壤类型主要有黄泥田、灰泥田及灰沙田等土属。为探讨三元非结构肥效模型的建模方法，本研究针对这两个项目县选择13个代表性水稻氮磷钾田间肥效试验结果作为研究案例。试验采用“3414”设计方案，即：1) N0P0K0；2)N0P2K2；3) N1P2K2；4) N2P0K2；5) N2P1K2；6)N2P2K2；7) N2P3K2；8) N2P2K0；9) N2P2K1；10)N2P2K3；11) N3P2K2；12) N1P1K2；13) N1P2K1；14)N2P1K1。其中，“2”水平为试验前当地氮磷钾推荐施肥量，“0”水平表示不施肥，“1”水平和“3”水平的施肥量分别为“2”水平的50%和150%。供试水稻品种采用当地大面积种植的良种，试验设计方案和田间管理措施与文献[12]相同。供试土壤主要理化性状采用常规方法[13]测定。基础土壤的主要理化性状和处理 6) 的施肥量以及各处理的试验产量结果，分别见表1和表2。

为客观地评价三元非结构肥效模型的科学性和实用性，截止2017年底，作者收集到近10年来福建早稻、晚稻和中稻的氮磷钾“3414”设计试验资料668个 (不包括表1和表2的相关试验点)。这些试验资料主要来自福建省水稻主产区完成的肥效试验结果 (表3)。田间试验设计方案与文献[12]相同或相似，不再赘述。

表1 早稻代表性试验点供试土壤理化性状及其处理氮磷钾施肥量Table 1 Physical and chemical properties of soils on the representative sites and fertilization rate

1.2 三元非结构肥效模型的构建

针对二次多项式肥效模型存在的模型设定偏误和多重共线性等问题[3,14]，章明清等[11]根据水稻氮磷钾单因素肥效试验结果，提出了一元非结构肥效模型：

式中：Y表示稻谷产量，X表示施肥量，s0为土壤供肥当量，计量单位均为kg/hm2。c为施肥增产效应系数；A表示施肥量X = 0时土壤肥力与稻谷产量之间的转换系数，综合反映了试验地的土壤生产力。在(1) 式模型中，当施肥量和土壤供肥当量都等于零时，作物产量必等于零。因此，根据植物营养元素功能不可相互替代的原理，三元非结构肥效模型可由如下最简化的数学形式来描述：

式中：A=AN×AP×AK。数学理论分析表明，(2) 式在一定施肥量范围内模型存在一个稻谷产量峰值，该峰值对应的施肥量即为最高产量施肥量。因此，根据微积分原理，令 (2) 式的水稻产量Y分别对N、P、K施肥量的导数等于零，得到最高产量施肥量计算式：

令 (2) 式的水稻产量Y分别对N、P、K施肥量的导数等于稻谷和肥料价格倒数比，得到经济产量施肥量计算式：

表2 早稻代表性试验点14个施肥处理的稻谷产量 (kg/hm2)Table 2 Yields of early rice in the 14 treatments in the representative experimental sites

表3 福建省各地区水稻氮磷钾田间肥效试验数(n)Table 3 Number of NPK fertilization experiments with “3414” design collected from cities in Fujian Province

1.3 三元非结构肥效模型的参数估计与统计检验

(2) 式模型是非线性模型，而且不能直接进行线性化处理，模型参数估计需采用非线性最小二乘法[15]。假设非线性肥效模型为Y =f(X,a)，为求得参数a的估计值，可求解最小二乘问题：

(2) 式模型的回归显著性检验与三元二次多项式肥效模型相似，但其回归自由度为6。本研究在计算机上的具体实现则使用MATLAB软件的nlinfit功能函数进行非线性肥效模型参数估计和统计检验，对三元二次多项式肥效模型回归分析则使用regress功能函数。文中图形采用MATLAB语言编程绘制，具体计算的数学原理和有关功能函数的使用方法可参阅相关专著[15-16]。

2 结果与分析

2.1 三元二次多项式肥效模型 (TPFM) 的拟合效果

常用的三元二次多项式肥效模型可用如下数学形式来表达：

式中：Y表示模型拟合产量；N、P、K分别表示N、P2O5、K2O施肥量 (kg/hm2)；b0至b9表示模型的肥效系数。根据表1中各试验点的氮磷钾施肥量及表2对应试验点各处理的稻谷产量，采用普通最小二乘法对 (6) 式进行回归建模 (表4)。结果表明，除1号和2号试验点外，其它各试验点建立的肥效模型均达到统计显著水平。

进一步对表4的三元二次多项式肥效模型进行典型性判别[10]，表明3号和4号试验点构建的肥效模型，虽然模型系数的代数符号满足要求，但模型不存在全局极大值，属于无最高产量点的非典型式；5号试验点构建的肥效模型，虽然模型系数的代数符号合理且模型存在最高产量点，但属于推荐施肥量外推的非典型式。6号至13号试验点建立的三元二次多项式肥效模型满足作物施肥效应的一般规律，均属于典型式。

2.2 三元非结构肥效模型 (TNFM) 的拟合效果

对表1和表2的同一批试验数据，采用 (2) 式进行非线性回归建模 (表5)。统计检验显示，1号和2号试验点的显著水平概率值P由 (6) 式模型未达显著水平的0.058和0.087提高到显著水平的0.017和0.010，回归模型通过了显著性检验；(6) 式模型能通过显著性检验的其它11个试验点，(2) 模型拟合结果同样能通过显著性检验，而且，除了10号试验点外，其它12个试验点的显著性概率值P均明显小于表4的 (6) 式模型相应指标。

表5的三元非结构肥效模型典型性判别[10]结果表明，由 (6) 式模型回归建模未能通过显著性检验的1号和2号试验点以及归属于非典型式的第3、4、5号试验点，均转化成了典型三元肥效模型，模型系数的代数符号、模型最高产量点和推荐施肥量等方面均满足了作物施肥效应的一般规律；(6) 式模型能得到典型式的6号至13号试验点，(2) 式模型的建模结果同样能得到典型式。

因此，统计显著性检验指标F值、拟合优度R2和显著水平概率值P以及典型式出现比例等方面都表明，(2) 式模型比 (6) 式模型具有更高的拟合精度和更宽的适用范围。

2.3 三元非结构肥效模型的推荐施肥量

根据表5中各试验点的三元非结构肥效模型的参数估计结果以及 (3) 式和 (4) 式，计算各试验点的氮、磷、钾最高施肥量，并以N 4.3元/kg、P2O55.0元/kg、K2O 5.0元/kg和稻谷2.0元/kg的平均市场价格计算经济施肥量。结果表明，1号至5号试验点的推荐施肥量均落在试验设计的施肥量范围内，相应的模型预测产量也都落在试验各处理的产量范围之内，无异常情况出现。

在根系养分吸收机理模型研究中，为评价机理模型养分吸收量预测值的可靠性，一般都采用在预测值与实测值之间进行一元线性回归分析[17]。本文采用相同方法考察（2）式模型在推荐施肥量上的可靠性，即：针对6号至13号试验点的（6）式模型参数估计结果，采用边际产量导数法计算相应试验点的氮、磷、钾的最高施肥量和经济施肥量，并与（2）式模型对应推荐施肥量绘制成图1。结果显示，两种肥效模型的氮、磷、钾的最高施肥之间以及经济施肥量之间都存在显著水平的线性正相关，说明在三元典型肥效模型前提下，新模型的推荐施肥量对二次多项式模式推荐施肥量具有继承性和可靠性。同时，线性回归模型的一次项系数分别为0.876和0.940，表明当三元二次多项式肥效模型推荐施肥量每增加1kg时，三元非结构肥效模型的最高施肥量和经济施肥量平均只增加了0.876kg和0.940kg，新模型较好地克服了二次多项式肥效模型推荐施肥量偏高的不足[2,18]。

图1 三元非结构肥效模型和三元二次多项式肥效模型推荐施肥量的相关性Fig. 1 Agreement of the recommended fertilizer rate between ternary non-structural and ternary quadratic polynomial fertilizer response models (TNFM and TPFM)

2.4 三元非结构肥效模型拟合效果评价

分别采用 (2) 式和 (6) 式对收集的 668个“3414”田间试验建三元肥效模型 (表6)。统计结果表明，目前常用的三元二次多项式肥效模型的典型式平均比例仅为19.5%，而三元非结构肥效模型的平均比例提高到39.1%，是前者的2.0倍，显著提升了建模成功率。

进一步分析还表明，受多重共线性危害的影响[3,19]，二次多项式肥效模型一次项或二次项的系数代数符号不合理的非典型式平均比例达到32.3%，而消除了多重共线性危害的三元非结构肥效模型的N0、P0、K0或c1、c2、c3等系数符号不合理的模型比例则平均为6.9%，大幅度降低了此类非典型式出现的机率。二次多项式肥效模型的一次项和二次项系数代数符号合理但模型无最高产量点的非典型式平均比例为14.4%；由于非结构肥效模型所具有的数学结构特点，模型系数代数符号合理则必有最高产量点。二次多项式模型系数代数符号合理且存在最高产量点但推荐施肥量外推的非典型式平均比例为4.0%，但非结构肥效模型的此类非典型式比例提高到30.7%，显示两种肥效模型在无最高产量点和推荐施肥量外推的非典型式类型出现比例方面有明显的差别。

3 讨论

3.1 三元二次多项式肥效模型的设定偏误

在肥料效应函数中，单位养分增产量与施肥量的函数关系的不同假设，就会得到数学形式各异和适用性不同的肥效模型[20]。假设增施单位量肥料的增产量和该养分最高产量施肥量与现有施肥量之差成正比，由此推导可得到一元二次多项式肥效模型。该模型被国内外的大量肥效试验尤其是氮肥试验结果所证实[20-23]。Cerrato等[24]、杨靖一等[25]和陈新平等[18]对蔬菜、冬小麦和夏玉米的多项式、线性或二次函数+平台等氮肥效应模型进行了比较研究，表明二次型模型具有较强的通用性。毛达如等[26]对2次、1.5次、0.75次和0.5次的氮磷二元多项式肥效模型比较也显示，在灌溉地上二次多项式能较好地反映冬小麦的肥效规律。可惜的是，668个水稻氮磷钾田间肥效试验表明，三元二次多项式肥效模型典型式的平均比例仅为19.5% (表6)，过低的建模成功率令人不得不怀疑二次多项式肥效模型设定本身的合理性。

对二次多项式的分析表明，该类模型是假设单位养分增产量与施肥量之间的函数关系为线性关系，结果导致最高施肥量之前和最高施肥量之后的施肥效应是对称关系[3]。这种模型设定忽略了高产耐肥新品种在过量施肥时因作物耐肥特性，使产量降低幅度得到较大程度缓解的当前生产实际。同时也忽略了土壤对养分的缓冲能力从而减轻了过量施肥对作物产量负效应的作用。此外，二次多项式模型的回归变量间存在强烈的多重共线性[3,14]，制约了普通最小二乘法回归建模的有效性和统计检验的可靠性。上述两个缺陷是导致当前常用的三元二次多项式肥效模型的典型式比例明显偏低的重要原因。

针对多项式统计模型的多重共线性问题，数理统计学家已经提出了诸如岭回归、主成分回归、偏最小二乘回归等有偏估计方法[19,27]，来消除或削弱这种多重共线性的危害。例如，采用主成分回归技术对福建171个早稻氮磷钾“3414”试验结果建立三元二次多项式肥效模型，其典型式的比例从普通最小二乘法的27.5%提高到43.3%[28]。但是，有偏估计只能消除或缓解多重共线性危害，不能有效地解决肥效模型设定的偏误问题。

表6 三元二次多项式肥效模型和三元非结构肥效模型拟合效果比较Table 6 Fitting effect of ternary non-structural fertilizer response model compared with that of ternary quadratic polynomial fertilizer response model

3.2 三元非典型式的类型和分类方法

在肥效模型研究中，常用的边际产量导数法计算推荐施肥量的方法仅适用于典型肥效模型，对非典型肥效模型的计算结果则会产生误导。因此，三元肥效模型典型性判别在推荐施肥中具有重要作用。总结福建省668个水稻氮磷钾田间肥效试验结果以及先前的相关田间肥效试验资料[10]，三元二次多项式肥效模型的非典型式有三种类型。在通过统计显著性检验前提下，从直观上看，如果肥效模型的一次项系数至少有一个为负数或二次项系数至少有一个为正数 (不考虑交互项系数符号)，因其不满足作物施肥效应的一般规律，此类肥效模型就属于模型系数不合理的非典型式，平均占到试验点总数的32.3% (表6)；第二种类型是模型一次项或二次项系数符号合理，但肥效模型不存在全局最高产量点，此类模型属于无最高产量点的非典型式，其平均比例占试验点总数的14.4% (表6)；第三种类型是模型一次项或二次项系数符号合理，而且肥效模型存在最高产量点，但推荐施肥量属于远外推结果，此类模型属于推荐施肥量外推的非典型式，其比例平均占试验点总数的4.0% (表6)。

与之对应，在表6的建模结果中，三元非结构肥效模型也会存在不同类型的非典型式。在通过统计显著性检验前提下，当模型参数N0、P0、K0、c1、c2、c3中有任何一个或一个以上的参数为负数时，因其违背了施肥效应的一般规律，该类肥效模型就属于参数符号不合理的非典型式，平均占试验点总数的6.9%，大幅度低于三元二次多项式肥效模型相同类型的比例。当模型参数N0、P0、K0、c1、c2、c3均为正数时，非结构肥效模型因数学结构特点必有最高产量点，因而不存在无最高产量点的非典型式类型。然而，模型参数代数符号合理但推荐施肥量属于远外推的非典型式类型平均达到试验点总数的30.7%，远高于三元二次多项式肥效模型相同类型的比例。可以想像，相信为降低远外推模型比例，非结构肥效模型对试验施肥量设计具有更高的要求。幸运的是，这一要求在试验设计中容易做到。

模型拟合过程还表明，“3414”试验采用多点分散不设重复的试验方法，无论是多项式肥效模型还是非结构肥效模型，都有相当大比例的试验点未能通过显著性检验，其中，二次多项式肥效模型的平均比例达到30.1%，非结构肥效模型则为23.4%。另外，未能通过统计显著性检验或者得到典型肥效模型的相关试验点经常会相对集中出现，说明建模成功与否还和田间试验质量密切相关。

3.3 三元非结构肥效模型的适用性

研究表明，一元二次多项式肥效模型是一元非结构肥效模型的简化式和特例[11]。针对三元非结构肥效模型的指数项，根据高等数学的泰勒展开式，即：其中，由于三元非结构肥效模型的参数c1、c2、c3均在10-3量级，若只取展开式的前两项，则 (2) 式模型可以转化为Y =A(N0+BN-c1N2) (P0+CP -c2P2) (K0+DK -c3K2) , 其中，B=1 -N0c1，C= 1 -P0c2，D= 1 -K0c3。对该式进行代数式展开，并忽略c1、c2、c3的两两乘积项以及c1c2c3的乘积项和NPK三因子交互项，则 (2) 式可以进一步转化为Y =A(N0P0K0+BP0K0N +CN0K0P +DN0P0K -c1P0K0N2-c2N0K0P2-c3N0P0K2+BCK0NP +BDP0NK +CDN0PK)，结果显示与 (6) 式具有相同的数学形式。由此可见，三元二次多项式肥效模型是三元非结构肥效模型的简化式和特例。当某些试验结果确实使上述忽略项对作物产量影响足够小时，(2) 式和 (6) 式模型都能得到很好的拟合效果；反之，三元二次多项式模型可能会因过分简化导致拟合效果较差，而三元非结构肥效模型因未进行这种简化则可能较好地拟合相关试验结果。

三元非结构肥效模型假设单位养分增产量与施肥量之间的函数关系为非线性关系，较好地克服了二次多项式肥效模型的设定偏误；模型本身不能直接进行线性化转换，较好地克服了多重共线性问题。在668个水稻氮磷钾田间肥效试验资料 (表6)中，三元非结构肥效模型的典型式平均比例达到39.1%，是三元二次多项式肥效模型典型式比例的2.0倍。图1分析表明，新模型推荐的最高施肥量和经济施肥量与三元二次多项式肥效模型的相应推荐施肥量具有显著水平的线性正相关，同时较好地克服了推荐施肥量偏高的问题。因此，新模型具有较宽的适用范围。

4 结论

三元非结构肥效模型较好地克服了三元二次多项式肥效模型的设定偏误和多重共线性危害，在水稻氮磷钾施肥效应建模中，具有更高的拟合精度和更宽的适用范围。