李 娟，章明清*，章贊德，許文江，姚寶全

（1 福建省農(nóng)業(yè)科學(xué)院土壤肥料研究所，福州 350013；2 福建省大田縣農(nóng)田建設(shè)與土壤肥料技術(shù)推廣站，福建 366100；3 福建省亞熱帶植物研究所，福建廈門 361000；4 福建省農(nóng)田建設(shè)與土壤肥料技術(shù)推廣總站，福州 350003）

當(dāng)前，通過構(gòu)建氮磷鉀三元肥效模型確定最佳施肥量，是實現(xiàn)水稻計量施肥的重要技術(shù)途徑之一[1-2]，其中，二次多項式函數(shù)是研究和應(yīng)用最多的肥效模型種類[3-7]。但是，眾多研究表明，一元和二元二次多項式肥效模型典型式的比例分別僅占60%左右和40.2%[8-9]，三元二次多項式肥效模型的典型式比例則更低至23.6%[10]。由于構(gòu)建肥效模型時出現(xiàn)了大量非典型式，嚴(yán)重削弱了該法的計量精確性和實用價值。盡管國內(nèi)外學(xué)者對此進(jìn)行了深入研究，提出了許多有意義的改進(jìn)措施[3]，但該問題至今仍然困擾著計量施肥研究和應(yīng)用?？v觀這些研究和改進(jìn)措施中，對二次多項式等函數(shù)形式的肥效模型本身是否存在模型設(shè)定偏誤及提出改進(jìn)建議還鮮見研究報道。

研究表明，當(dāng)前廣泛應(yīng)用的一元、二元、三元二次多項式肥效模型及其它類似的多項式肥效模型存在明顯的設(shè)定偏誤[3]。為此，章明清等[11]提出了一元非結(jié)構(gòu)肥效模型，較好地克服了模型設(shè)定缺陷。與一元二次多項式肥效模型相比，新模型在擬合水稻氮磷鉀單因素田間肥效試驗結(jié)果時，具有較高的擬合精度、較寬的適用范圍和推薦施肥量較低等優(yōu)點。為此，在一元非結(jié)構(gòu)肥效模型基礎(chǔ)上，本研究利用氮磷鉀三因素田間肥效試驗結(jié)果，探討三元非結(jié)構(gòu)肥效模型的構(gòu)建方法及其對田間肥效試驗資料的擬合效果，旨在擴大三元肥效模型的適用性，為計量施肥研究和應(yīng)用提供一種新的模型方法。

1 材料與方法

1.1 水稻氮磷鉀田間肥效試驗設(shè)計

近10年來，筆者參加福建省莆田市仙游縣和漳州市平和縣的水稻測土配方施肥工作。兩個項目縣均地處南亞熱帶海洋性季風(fēng)氣候帶，氣候條件尤其適合水稻生長發(fā)育。 稻田土壤類型主要有黃泥田、灰泥田及灰沙田等土屬。為探討三元非結(jié)構(gòu)肥效模型的建模方法，本研究針對這兩個項目縣選擇13個代表性水稻氮磷鉀田間肥效試驗結(jié)果作為研究案例。試驗采用“3414”設(shè)計方案，即：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”水平為試驗前當(dāng)?shù)氐租浲扑]施肥量，“0”水平表示不施肥，“1”水平和“3”水平的施肥量分別為“2”水平的50%和150%。供試水稻品種采用當(dāng)?shù)卮竺娣e種植的良種，試驗設(shè)計方案和田間管理措施與文獻(xiàn)[12]相同。供試土壤主要理化性狀采用常規(guī)方法[13]測定?；A(chǔ)土壤的主要理化性狀和處理 6) 的施肥量以及各處理的試驗產(chǎn)量結(jié)果，分別見表1和表2。

為客觀地評價三元非結(jié)構(gòu)肥效模型的科學(xué)性和實用性，截止2017年底，作者收集到近10年來福建早稻、晚稻和中稻的氮磷鉀“3414”設(shè)計試驗資料668個 (不包括表1和表2的相關(guān)試驗點)。這些試驗資料主要來自福建省水稻主產(chǎn)區(qū)完成的肥效試驗結(jié)果 (表3)。田間試驗設(shè)計方案與文獻(xiàn)[12]相同或相似，不再贅述。

表1 早稻代表性試驗點供試土壤理化性狀及其處理氮磷鉀施肥量Table 1 Physical and chemical properties of soils on the representative sites and fertilization rate

1.2 三元非結(jié)構(gòu)肥效模型的構(gòu)建

針對二次多項式肥效模型存在的模型設(shè)定偏誤和多重共線性等問題[3,14]，章明清等[11]根據(jù)水稻氮磷鉀單因素肥效試驗結(jié)果，提出了一元非結(jié)構(gòu)肥效模型：

式中：Y表示稻谷產(chǎn)量，X表示施肥量，s0為土壤供肥當(dāng)量，計量單位均為kg/hm2。c為施肥增產(chǎn)效應(yīng)系數(shù)；A表示施肥量X = 0時土壤肥力與稻谷產(chǎn)量之間的轉(zhuǎn)換系數(shù)，綜合反映了試驗地的土壤生產(chǎn)力。在(1) 式模型中，當(dāng)施肥量和土壤供肥當(dāng)量都等于零時，作物產(chǎn)量必等于零。因此，根據(jù)植物營養(yǎng)元素功能不可相互替代的原理，三元非結(jié)構(gòu)肥效模型可由如下最簡化的數(shù)學(xué)形式來描述：

式中：A=AN×AP×AK。數(shù)學(xué)理論分析表明，(2) 式在一定施肥量范圍內(nèi)模型存在一個稻谷產(chǎn)量峰值，該峰值對應(yīng)的施肥量即為最高產(chǎn)量施肥量。因此，根據(jù)微積分原理，令 (2) 式的水稻產(chǎn)量Y分別對N、P、K施肥量的導(dǎo)數(shù)等于零，得到最高產(chǎn)量施肥量計算式：

令 (2) 式的水稻產(chǎn)量Y分別對N、P、K施肥量的導(dǎo)數(shù)等于稻谷和肥料價格倒數(shù)比，得到經(jīng)濟(jì)產(chǎn)量施肥量計算式：

表2 早稻代表性試驗點14個施肥處理的稻谷產(chǎn)量 (kg/hm2)Table 2 Yields of early rice in the 14 treatments in the representative experimental sites

表3 福建省各地區(qū)水稻氮磷鉀田間肥效試驗數(shù)(n)Table 3 Number of NPK fertilization experiments with “3414” design collected from cities in Fujian Province

1.3 三元非結(jié)構(gòu)肥效模型的參數(shù)估計與統(tǒng)計檢驗

(2) 式模型是非線性模型，而且不能直接進(jìn)行線性化處理，模型參數(shù)估計需采用非線性最小二乘法[15]。假設(shè)非線性肥效模型為Y =f(X,a)，為求得參數(shù)a的估計值，可求解最小二乘問題：

(2) 式模型的回歸顯著性檢驗與三元二次多項式肥效模型相似，但其回歸自由度為6。本研究在計算機上的具體實現(xiàn)則使用MATLAB軟件的nlinfit功能函數(shù)進(jìn)行非線性肥效模型參數(shù)估計和統(tǒng)計檢驗，對三元二次多項式肥效模型回歸分析則使用regress功能函數(shù)。文中圖形采用MATLAB語言編程繪制，具體計算的數(shù)學(xué)原理和有關(guān)功能函數(shù)的使用方法可參閱相關(guān)專著[15-16]。

2 結(jié)果與分析

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

常用的三元二次多項式肥效模型可用如下數(shù)學(xué)形式來表達(dá)：

式中：Y表示模型擬合產(chǎn)量；N、P、K分別表示N、P2O5、K2O施肥量 (kg/hm2)；b0至b9表示模型的肥效系數(shù)。根據(jù)表1中各試驗點的氮磷鉀施肥量及表2對應(yīng)試驗點各處理的稻谷產(chǎn)量，采用普通最小二乘法對 (6) 式進(jìn)行回歸建模 (表4)。結(jié)果表明，除1號和2號試驗點外，其它各試驗點建立的肥效模型均達(dá)到統(tǒng)計顯著水平。

進(jìn)一步對表4的三元二次多項式肥效模型進(jìn)行典型性判別[10]，表明3號和4號試驗點構(gòu)建的肥效模型，雖然模型系數(shù)的代數(shù)符號滿足要求，但模型不存在全局極大值，屬于無最高產(chǎn)量點的非典型式；5號試驗點構(gòu)建的肥效模型，雖然模型系數(shù)的代數(shù)符號合理且模型存在最高產(chǎn)量點，但屬于推薦施肥量外推的非典型式。6號至13號試驗點建立的三元二次多項式肥效模型滿足作物施肥效應(yīng)的一般規(guī)律，均屬于典型式。

2.2 三元非結(jié)構(gòu)肥效模型 (TNFM) 的擬合效果

對表1和表2的同一批試驗數(shù)據(jù)，采用 (2) 式進(jìn)行非線性回歸建模 (表5)。統(tǒng)計檢驗顯示，1號和2號試驗點的顯著水平概率值P由 (6) 式模型未達(dá)顯著水平的0.058和0.087提高到顯著水平的0.017和0.010，回歸模型通過了顯著性檢驗；(6) 式模型能通過顯著性檢驗的其它11個試驗點，(2) 模型擬合結(jié)果同樣能通過顯著性檢驗，而且，除了10號試驗點外，其它12個試驗點的顯著性概率值P均明顯小于表4的 (6) 式模型相應(yīng)指標(biāo)。

表5的三元非結(jié)構(gòu)肥效模型典型性判別[10]結(jié)果表明，由 (6) 式模型回歸建模未能通過顯著性檢驗的1號和2號試驗點以及歸屬于非典型式的第3、4、5號試驗點，均轉(zhuǎn)化成了典型三元肥效模型，模型系數(shù)的代數(shù)符號、模型最高產(chǎn)量點和推薦施肥量等方面均滿足了作物施肥效應(yīng)的一般規(guī)律；(6) 式模型能得到典型式的6號至13號試驗點，(2) 式模型的建模結(jié)果同樣能得到典型式。

因此，統(tǒng)計顯著性檢驗指標(biāo)F值、擬合優(yōu)度R2和顯著水平概率值P以及典型式出現(xiàn)比例等方面都表明，(2) 式模型比 (6) 式模型具有更高的擬合精度和更寬的適用范圍。

2.3 三元非結(jié)構(gòu)肥效模型的推薦施肥量

根據(jù)表5中各試驗點的三元非結(jié)構(gòu)肥效模型的參數(shù)估計結(jié)果以及 (3) 式和 (4) 式，計算各試驗點的氮、磷、鉀最高施肥量，并以N 4.3元/kg、P2O55.0元/kg、K2O 5.0元/kg和稻谷2.0元/kg的平均市場價格計算經(jīng)濟(jì)施肥量。結(jié)果表明，1號至5號試驗點的推薦施肥量均落在試驗設(shè)計的施肥量范圍內(nèi)，相應(yīng)的模型預(yù)測產(chǎn)量也都落在試驗各處理的產(chǎn)量范圍之內(nèi)，無異常情況出現(xiàn)。

在根系養(yǎng)分吸收機理模型研究中，為評價機理模型養(yǎng)分吸收量預(yù)測值的可靠性，一般都采用在預(yù)測值與實測值之間進(jìn)行一元線性回歸分析[17]。本文采用相同方法考察（2）式模型在推薦施肥量上的可靠性，即：針對6號至13號試驗點的（6）式模型參數(shù)估計結(jié)果，采用邊際產(chǎn)量導(dǎo)數(shù)法計算相應(yīng)試驗點的氮、磷、鉀的最高施肥量和經(jīng)濟(jì)施肥量，并與（2）式模型對應(yīng)推薦施肥量繪制成圖1。結(jié)果顯示，兩種肥效模型的氮、磷、鉀的最高施肥之間以及經(jīng)濟(jì)施肥量之間都存在顯著水平的線性正相關(guān)，說明在三元典型肥效模型前提下，新模型的推薦施肥量對二次多項式模式推薦施肥量具有繼承性和可靠性。同時，線性回歸模型的一次項系數(shù)分別為0.876和0.940，表明當(dāng)三元二次多項式肥效模型推薦施肥量每增加1kg時，三元非結(jié)構(gòu)肥效模型的最高施肥量和經(jīng)濟(jì)施肥量平均只增加了0.876kg和0.940kg，新模型較好地克服了二次多項式肥效模型推薦施肥量偏高的不足[2,18]。

圖1 三元非結(jié)構(gòu)肥效模型和三元二次多項式肥效模型推薦施肥量的相關(guān)性Fig. 1 Agreement of the recommended fertilizer rate between ternary non-structural and ternary quadratic polynomial fertilizer response models (TNFM and TPFM)

2.4 三元非結(jié)構(gòu)肥效模型擬合效果評價

分別采用 (2) 式和 (6) 式對收集的 668個“3414”田間試驗建三元肥效模型 (表6)。統(tǒng)計結(jié)果表明，目前常用的三元二次多項式肥效模型的典型式平均比例僅為19.5%，而三元非結(jié)構(gòu)肥效模型的平均比例提高到39.1%，是前者的2.0倍，顯著提升了建模成功率。

進(jìn)一步分析還表明，受多重共線性危害的影響[3,19]，二次多項式肥效模型一次項或二次項的系數(shù)代數(shù)符號不合理的非典型式平均比例達(dá)到32.3%，而消除了多重共線性危害的三元非結(jié)構(gòu)肥效模型的N0、P0、K0或c1、c2、c3等系數(shù)符號不合理的模型比例則平均為6.9%，大幅度降低了此類非典型式出現(xiàn)的機率。二次多項式肥效模型的一次項和二次項系數(shù)代數(shù)符號合理但模型無最高產(chǎn)量點的非典型式平均比例為14.4%；由于非結(jié)構(gòu)肥效模型所具有的數(shù)學(xué)結(jié)構(gòu)特點，模型系數(shù)代數(shù)符號合理則必有最高產(chǎn)量點。二次多項式模型系數(shù)代數(shù)符號合理且存在最高產(chǎn)量點但推薦施肥量外推的非典型式平均比例為4.0%，但非結(jié)構(gòu)肥效模型的此類非典型式比例提高到30.7%，顯示兩種肥效模型在無最高產(chǎn)量點和推薦施肥量外推的非典型式類型出現(xiàn)比例方面有明顯的差別。

3 討論

3.1 三元二次多項式肥效模型的設(shè)定偏誤

在肥料效應(yīng)函數(shù)中，單位養(yǎng)分增產(chǎn)量與施肥量的函數(shù)關(guān)系的不同假設(shè)，就會得到數(shù)學(xué)形式各異和適用性不同的肥效模型[20]。假設(shè)增施單位量肥料的增產(chǎn)量和該養(yǎng)分最高產(chǎn)量施肥量與現(xiàn)有施肥量之差成正比，由此推導(dǎo)可得到一元二次多項式肥效模型。該模型被國內(nèi)外的大量肥效試驗尤其是氮肥試驗結(jié)果所證實[20-23]。Cerrato等[24]、楊靖一等[25]和陳新平等[18]對蔬菜、冬小麥和夏玉米的多項式、線性或二次函數(shù)+平臺等氮肥效應(yīng)模型進(jìn)行了比較研究，表明二次型模型具有較強的通用性。毛達(dá)如等[26]對2次、1.5次、0.75次和0.5次的氮磷二元多項式肥效模型比較也顯示，在灌溉地上二次多項式能較好地反映冬小麥的肥效規(guī)律。可惜的是，668個水稻氮磷鉀田間肥效試驗表明，三元二次多項式肥效模型典型式的平均比例僅為19.5% (表6)，過低的建模成功率令人不得不懷疑二次多項式肥效模型設(shè)定本身的合理性。

對二次多項式的分析表明，該類模型是假設(shè)單位養(yǎng)分增產(chǎn)量與施肥量之間的函數(shù)關(guān)系為線性關(guān)系，結(jié)果導(dǎo)致最高施肥量之前和最高施肥量之后的施肥效應(yīng)是對稱關(guān)系[3]。這種模型設(shè)定忽略了高產(chǎn)耐肥新品種在過量施肥時因作物耐肥特性，使產(chǎn)量降低幅度得到較大程度緩解的當(dāng)前生產(chǎn)實際。同時也忽略了土壤對養(yǎng)分的緩沖能力從而減輕了過量施肥對作物產(chǎn)量負(fù)效應(yīng)的作用。此外，二次多項式模型的回歸變量間存在強烈的多重共線性[3,14]，制約了普通最小二乘法回歸建模的有效性和統(tǒng)計檢驗的可靠性。上述兩個缺陷是導(dǎo)致當(dāng)前常用的三元二次多項式肥效模型的典型式比例明顯偏低的重要原因。

針對多項式統(tǒng)計模型的多重共線性問題，數(shù)理統(tǒng)計學(xué)家已經(jīng)提出了諸如嶺回歸、主成分回歸、偏最小二乘回歸等有偏估計方法[19,27]，來消除或削弱這種多重共線性的危害。例如，采用主成分回歸技術(shù)對福建171個早稻氮磷鉀“3414”試驗結(jié)果建立三元二次多項式肥效模型，其典型式的比例從普通最小二乘法的27.5%提高到43.3%[28]。但是，有偏估計只能消除或緩解多重共線性危害，不能有效地解決肥效模型設(shè)定的偏誤問題。

表6 三元二次多項式肥效模型和三元非結(jié)構(gòu)肥效模型擬合效果比較Table 6 Fitting effect of ternary non-structural fertilizer response model compared with that of ternary quadratic polynomial fertilizer response model

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

在肥效模型研究中，常用的邊際產(chǎn)量導(dǎo)數(shù)法計算推薦施肥量的方法僅適用于典型肥效模型，對非典型肥效模型的計算結(jié)果則會產(chǎn)生誤導(dǎo)。因此，三元肥效模型典型性判別在推薦施肥中具有重要作用?？偨Y(jié)福建省668個水稻氮磷鉀田間肥效試驗結(jié)果以及先前的相關(guān)田間肥效試驗資料[10]，三元二次多項式肥效模型的非典型式有三種類型。在通過統(tǒng)計顯著性檢驗前提下，從直觀上看，如果肥效模型的一次項系數(shù)至少有一個為負(fù)數(shù)或二次項系數(shù)至少有一個為正數(shù) (不考慮交互項系數(shù)符號)，因其不滿足作物施肥效應(yīng)的一般規(guī)律，此類肥效模型就屬于模型系數(shù)不合理的非典型式，平均占到試驗點總數(shù)的32.3% (表6)；第二種類型是模型一次項或二次項系數(shù)符號合理，但肥效模型不存在全局最高產(chǎn)量點，此類模型屬于無最高產(chǎn)量點的非典型式，其平均比例占試驗點總數(shù)的14.4% (表6)；第三種類型是模型一次項或二次項系數(shù)符號合理，而且肥效模型存在最高產(chǎn)量點，但推薦施肥量屬于遠(yuǎn)外推結(jié)果，此類模型屬于推薦施肥量外推的非典型式，其比例平均占試驗點總數(shù)的4.0% (表6)。

與之對應(yīng)，在表6的建模結(jié)果中，三元非結(jié)構(gòu)肥效模型也會存在不同類型的非典型式。在通過統(tǒng)計顯著性檢驗前提下，當(dāng)模型參數(shù)N0、P0、K0、c1、c2、c3中有任何一個或一個以上的參數(shù)為負(fù)數(shù)時，因其違背了施肥效應(yīng)的一般規(guī)律，該類肥效模型就屬于參數(shù)符號不合理的非典型式，平均占試驗點總數(shù)的6.9%，大幅度低于三元二次多項式肥效模型相同類型的比例。當(dāng)模型參數(shù)N0、P0、K0、c1、c2、c3均為正數(shù)時，非結(jié)構(gòu)肥效模型因數(shù)學(xué)結(jié)構(gòu)特點必有最高產(chǎn)量點，因而不存在無最高產(chǎn)量點的非典型式類型。然而，模型參數(shù)代數(shù)符號合理但推薦施肥量屬于遠(yuǎn)外推的非典型式類型平均達(dá)到試驗點總數(shù)的30.7%，遠(yuǎn)高于三元二次多項式肥效模型相同類型的比例?？梢韵胂?，相信為降低遠(yuǎn)外推模型比例，非結(jié)構(gòu)肥效模型對試驗施肥量設(shè)計具有更高的要求。幸運的是，這一要求在試驗設(shè)計中容易做到。

模型擬合過程還表明，“3414”試驗采用多點分散不設(shè)重復(fù)的試驗方法，無論是多項式肥效模型還是非結(jié)構(gòu)肥效模型，都有相當(dāng)大比例的試驗點未能通過顯著性檢驗，其中，二次多項式肥效模型的平均比例達(dá)到30.1%，非結(jié)構(gòu)肥效模型則為23.4%。另外，未能通過統(tǒng)計顯著性檢驗或者得到典型肥效模型的相關(guān)試驗點經(jīng)常會相對集中出現(xiàn)，說明建模成功與否還和田間試驗質(zhì)量密切相關(guān)。

3.3 三元非結(jié)構(gòu)肥效模型的適用性

研究表明，一元二次多項式肥效模型是一元非結(jié)構(gòu)肥效模型的簡化式和特例[11]。針對三元非結(jié)構(gòu)肥效模型的指數(shù)項，根據(jù)高等數(shù)學(xué)的泰勒展開式，即：其中，由于三元非結(jié)構(gòu)肥效模型的參數(shù)c1、c2、c3均在10-3量級，若只取展開式的前兩項，則 (2) 式模型可以轉(zhuǎn)化為Y =A(N0+BN-c1N2) (P0+CP -c2P2) (K0+DK -c3K2) , 其中，B=1 -N0c1，C= 1 -P0c2，D= 1 -K0c3。對該式進(jìn)行代數(shù)式展開，并忽略c1、c2、c3的兩兩乘積項以及c1c2c3的乘積項和NPK三因子交互項，則 (2) 式可以進(jìn)一步轉(zhuǎn)化為Y =A(N0P0K0+BP0K0N +CN0K0P +DN0P0K -c1P0K0N2-c2N0K0P2-c3N0P0K2+BCK0NP +BDP0NK +CDN0PK)，結(jié)果顯示與 (6) 式具有相同的數(shù)學(xué)形式。由此可見，三元二次多項式肥效模型是三元非結(jié)構(gòu)肥效模型的簡化式和特例。當(dāng)某些試驗結(jié)果確實使上述忽略項對作物產(chǎn)量影響足夠小時，(2) 式和 (6) 式模型都能得到很好的擬合效果；反之，三元二次多項式模型可能會因過分簡化導(dǎo)致擬合效果較差，而三元非結(jié)構(gòu)肥效模型因未進(jìn)行這種簡化則可能較好地擬合相關(guān)試驗結(jié)果。

三元非結(jié)構(gòu)肥效模型假設(shè)單位養(yǎng)分增產(chǎn)量與施肥量之間的函數(shù)關(guān)系為非線性關(guān)系，較好地克服了二次多項式肥效模型的設(shè)定偏誤；模型本身不能直接進(jìn)行線性化轉(zhuǎn)換，較好地克服了多重共線性問題。在668個水稻氮磷鉀田間肥效試驗資料 (表6)中，三元非結(jié)構(gòu)肥效模型的典型式平均比例達(dá)到39.1%，是三元二次多項式肥效模型典型式比例的2.0倍。圖1分析表明，新模型推薦的最高施肥量和經(jīng)濟(jì)施肥量與三元二次多項式肥效模型的相應(yīng)推薦施肥量具有顯著水平的線性正相關(guān)，同時較好地克服了推薦施肥量偏高的問題。因此，新模型具有較寬的適用范圍。

4 結(jié)論

三元非結(jié)構(gòu)肥效模型較好地克服了三元二次多項式肥效模型的設(shè)定偏誤和多重共線性危害，在水稻氮磷鉀施肥效應(yīng)建模中，具有更高的擬合精度和更寬的適用范圍。