孔慶波，章明清，李娟，許文江，章贊德，姚建族

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

0 引言

【研究意義】肥料效應(yīng)函數(shù)法是中國測土配方施肥技術(shù)的一個重要分支體系。但根據(jù)田間肥效試驗(yàn)結(jié)果，應(yīng)用普通最小二乘法（OLS）建立肥效模型的成功率明顯偏低，嚴(yán)重制約了該法的計量精確性和實(shí)用價值。因此，探討如何提高田間肥效試驗(yàn)建模成功率具有意義。【前人研究進(jìn)展】在常用肥效模型中，二次多項(xiàng)式肥效方程能較好地反映作物產(chǎn)量和施肥量之間的數(shù)量關(guān)系，模型計算和參數(shù)估計簡便易行，是研究和應(yīng)用最多的模型種類[1-3]。但在施肥實(shí)踐中，OLS法回歸分析建立的一元肥效模型典型式僅占 60%左右，二元肥效模型典型式則只有40.2%[4-5]，而三元肥效模型的典型式更低至 23.6%[6]。針對肥效模型在應(yīng)用中存在的建模成功率偏低的問題，國內(nèi)外學(xué)者就肥效模型選擇與其適用性[7-9]、試驗(yàn)設(shè)計與參數(shù)估計方法[10-12]、類特征肥效模型構(gòu)建[13-15]及非典型式推薦施肥優(yōu)化方法[12]等方面進(jìn)行了深入研究，提出了許多改進(jìn)措施，但至今未能得到較好地解決。當(dāng)前廣泛應(yīng)用的“3414”設(shè)計方案的肥效試驗(yàn)結(jié)果既使采用蒙特卡洛（Monter Carlo）建模法，二元和三元模型典型式比率也僅有 56.7%和 37.3%[12]。這種成敗幾乎各半的建模結(jié)果，困擾著計量施肥研究和應(yīng)用。雖然任何田間試驗(yàn)都可能存在失誤從而導(dǎo)致失敗，但是，在不同肥料之間、肥料與其他產(chǎn)量限制因子之間存在著復(fù)雜的聯(lián)應(yīng)關(guān)系，我們沒有理由認(rèn)為其余的過半數(shù)試驗(yàn)均是失誤。【本研究切入點(diǎn)】隨著回歸分析理論的發(fā)展和完善，人們發(fā)現(xiàn)，當(dāng)回歸方程隨機(jī)誤差項(xiàng)的方差存在顯著水平差異時，普通最小二乘法回歸建模就不再適用[16]。此時其參數(shù)估計值就不再具有最小方差的優(yōu)良性，模型參數(shù)顯著性的t檢驗(yàn)、回歸方程顯著性的F檢驗(yàn)和擬合優(yōu)度R2檢驗(yàn)等都將失去可靠性，同時擬合參數(shù)具有較大的方差，回歸方程預(yù)測能力也將失效，結(jié)果導(dǎo)致肥效模型失去實(shí)用價值?！緮M解決的關(guān)鍵問題】本文利用近年來“3414”試驗(yàn)設(shè)計方案完成的早稻氮磷鉀田間肥效試驗(yàn)結(jié)果，探討三元二次多項(xiàng)式肥效模型的異方差診斷方法和可行廣義最小二乘回歸建模法（FGLS）在消除或緩解異方差危害方面的應(yīng)用效果，為提高三元肥效模型建?？煽啃院统晒β侍峁┮环N新方法。

1 材料與方法

1.1 早稻氮磷鉀肥效試驗(yàn)資料的收集整理

近年來，在水稻測土配方施肥工作中，福建在福州市、寧德市、南平市、三明市、龍巖市、莆田市、漳州的龍海市以及平和縣、泉州的南安市等水稻主產(chǎn)區(qū)的主要項(xiàng)目縣先后完成了171個氮磷鉀田間肥效試驗(yàn)。這些試驗(yàn)均采用“3414”設(shè)計方案[17]，即（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”水平為試驗(yàn)前的當(dāng)?shù)赝扑]施肥量，“0”水平表示不施肥，“1”水平和“3”水平的施肥量分別為“2”水平的50%和150%。本研究選擇8個代表性試驗(yàn)點(diǎn)的田間試驗(yàn)結(jié)果來探討三元二次多項(xiàng)式肥效模型的異方差診斷方法及其應(yīng)對技術(shù)。這些試驗(yàn)點(diǎn)基礎(chǔ)土樣的主要理化性狀采用常規(guī)方法[18]測定，處理（6）的施肥量采用當(dāng)?shù)赝扑]施肥量（表 1），試驗(yàn)設(shè)計方案和田間管理與文獻(xiàn)[17]相同。8個試驗(yàn)點(diǎn)的產(chǎn)量結(jié)果見表2。

1.2 多元線性回歸模型的異方差及其診斷

經(jīng)典線性回歸分析針對模型隨機(jī)誤差項(xiàng)u?的一個重要假設(shè)是u?在不同觀測點(diǎn)具有相同方差，即誤差項(xiàng)u?的方差是一個等于σ2的常數(shù)[16]。但是，這種等方差假設(shè)在實(shí)際問題中往往難以得到滿足。當(dāng)回歸模型存在嚴(yán)重異方差時，普通最小二乘法就喪失了實(shí)用價值[19-20]。異方差的定性診斷方法是模型殘差分析[16]。當(dāng)回歸模型的誤差項(xiàng)具有同方差時，殘差圖上n個點(diǎn)的散布是隨機(jī)的，無任何規(guī)律。如果存在異方差，殘差圖上點(diǎn)的分布就會呈現(xiàn)出一定的變化趨勢，如殘差值增大或減小等。

因線性回歸模型的異方差有多種不同表現(xiàn)形式[19]，其定量診斷方法也有多種不同方法。本研究選擇常用的懷特（White）檢驗(yàn)、帕克（Park）檢驗(yàn)和冠因克-巴塞特檢驗(yàn)（Koenker-Basett，KB檢驗(yàn)）等，具體使用方法可參考相關(guān)專著[19]。

表1 早稻8個代表性試驗(yàn)點(diǎn)的供試土壤主要理化性狀及其處理（6）氮磷鉀施肥量Table 1 Main physical and chemical properties of soils in 8 representative experience sites and fertilization rate in 6th treatment in early rice

表2 早稻8個代表性試驗(yàn)點(diǎn)各施肥處理的稻谷產(chǎn)量Table 2 Rice grain yield of each fertilizer treatment in 8 representative experiments in early rice

1.3 可行廣義最小二乘回歸建模法

目前常用的三元二次多項(xiàng)式肥效模型可用矩陣形式表示為：

式中，Y?表示擬合產(chǎn)量；X表示肥效試驗(yàn)設(shè)計矩陣，β表示肥效模型參數(shù)向量，u?表示模型擬合殘差向量。用普通最小二乘法（OLS）進(jìn)行回歸分析時，模型參數(shù)的計算式為：式中，Y表示各施肥處理的試驗(yàn)產(chǎn)量向量。針對異方差問題，數(shù)理統(tǒng)計學(xué)家已經(jīng)提出了消除或緩解線性回歸模型異方差的許多方法[19,21]，但應(yīng)用較廣泛的方法是可行廣義最小二乘法（FGLS）或稱之為可行加權(quán)最小二乘法（FWLS）。這種方法是對原回歸模型（1）式加權(quán)，使之變成一個新的不存在異方差性的模型，然后采用OLS估計其參數(shù)[19]。因此，F(xiàn)GLS就是對加了權(quán)重的殘差平方和實(shí)施OLS，即：

式中，Wi為給定的第i個處理產(chǎn)量的權(quán)數(shù)。如何獲得各處理產(chǎn)量的 Wi？一種可行的方法是對原模型進(jìn)行OLS估計，得到隨機(jī)誤差項(xiàng)的近似估計量u?，以此構(gòu)建權(quán)數(shù)估計量。即：，進(jìn)而構(gòu)造權(quán)矩陣W，然后用OLS估計新模型。記參數(shù)估計量為，即FGLS的參數(shù)估計量計算式為：

這就是原模型（1）式的可行廣義最小二乘估計，它是無偏和有效的估計量[16]。因?yàn)?*β是β的最佳線性無偏估計量，因此，F(xiàn)GLS的總離差平方和SST*、回歸平方和 SSR*以及殘差平方和 SSE*的計算式分別為：

進(jìn)而可以計算得到FGLS回歸模型的F值和擬合優(yōu)度R2等統(tǒng)計檢驗(yàn)參數(shù)。上述的具體數(shù)學(xué)原理和計算方法可參閱相關(guān)專著[16,21-22]。

1.4 肥效模型典型性判別方法

假設(shè)三元二次多項(xiàng)式肥效模型的數(shù)學(xué)形式如下：

式中，N、P、K 分別表示 N、P2O5和 K2O 施肥量（kg·hm-2）；ε為擬合殘差（kg·hm-2）；Y?為擬合產(chǎn)量（kg·hm-2）。根據(jù)無約束最優(yōu)化方法[6]，若函數(shù)f (N,P, K)在點(diǎn)X處的一階梯度向量g(X) （其中，X=（N,P, K））等于零向量；同時，若肥效函數(shù)的Hesse矩陣G(X )為負(fù)定，則此函數(shù)有極大值。對式（8）的三元二次多項(xiàng)式肥效模型而言，Hesse矩陣G(X)的各階主子式假設(shè)為G1、G2和G3，則它們的行列式值分別為：G1=2b4；G2=4b4b5-b72；G3=2(4b4b5b6+b7b8b9-b4b92-b5b82-b6b72)。因此，各類三元肥效模型的典型性判別方法如下：①若g(X* )=0，且G1＜0，G2＞0，G3＜0，則矩陣G(X )為負(fù)定，肥效模型具有全局最高產(chǎn)量點(diǎn)。同時，若g(X*)=0所代表的點(diǎn)X*對應(yīng)的N、P、K數(shù)值均落在試驗(yàn)設(shè)計的施肥量范圍內(nèi)，則該肥效模型屬于典型式。②若G(X*)為負(fù)定，但點(diǎn)X*對應(yīng)的N、P、K數(shù)值有一個或一個以上超過試驗(yàn)設(shè)計施肥量，則該肥效模型屬于推薦施肥量外推的非典型式。③若g(X*)=0，而且 G1＞0，G2＞0，G3＞0，則 Hesse矩陣G (X )為正定，此函數(shù)有全局最低產(chǎn)量點(diǎn)，肥效型屬于存在最低產(chǎn)量點(diǎn)的非典型式。④若g(X*)=0，G1、G2和 G3不滿足正定和負(fù)定的條件而且不等于零，則Hesse矩陣為不定，點(diǎn)X*為駐點(diǎn)，該肥效模型屬于無最高產(chǎn)量點(diǎn)的非典型式。

1.5 肥效模型異方差檢驗(yàn)和參數(shù)估計的計算機(jī)實(shí)現(xiàn)

三元二次多項(xiàng)式肥效模型的異方差檢驗(yàn)、OLS和FGLS回歸建模均涉及到繁雜的數(shù)學(xué)計算，本研究采用MATLAB軟件來完成相關(guān)工作，其中，OLS法回歸診斷調(diào)用 regstats功能函數(shù)，OLS回歸分析調(diào)用regress功能函數(shù)，F(xiàn)GLS法則調(diào)用fgls功能函數(shù)。根據(jù)上述統(tǒng)計學(xué)原理，在 MATLAB軟件平臺上編寫計算程序完成異方差檢驗(yàn)及其繪圖。

2 結(jié)果

2.1 OLS法回歸分析及其肥效模型典型性判別

根據(jù)表1的8個試驗(yàn)點(diǎn)施肥量和表2相應(yīng)的14個處理試驗(yàn)產(chǎn)量結(jié)果，應(yīng)用OLS進(jìn)行回歸分析（表3）。結(jié)果表明，1號和2號的兩個試驗(yàn)點(diǎn)未達(dá)到統(tǒng)計顯著水平，試驗(yàn)資料不可用；3—8號的6個試驗(yàn)點(diǎn)，回歸模型均達(dá)到統(tǒng)計顯著水平。因?yàn)橥ㄟ^顯著性檢驗(yàn)的三元二次多項(xiàng)式肥效模型有典型式和3種不同類別的非典型式[6]，因而，通過顯著性檢驗(yàn)的 6個回歸模型還需進(jìn)一步進(jìn)行典型性判別（表 3），以便確定推薦施肥方法。結(jié)果顯示，盡管3號試驗(yàn)點(diǎn)肥效模型的一次項(xiàng)參數(shù)符號為正，二次項(xiàng)參數(shù)符號為負(fù)，滿足了植物營養(yǎng)的一般規(guī)律，但是，該肥效模型不存在最高產(chǎn)量點(diǎn)，邊際產(chǎn)量導(dǎo)數(shù)法推薦施肥的結(jié)果不可靠。4號試驗(yàn)點(diǎn)有最高產(chǎn)量點(diǎn)，邊際產(chǎn)量導(dǎo)數(shù)法推薦施肥量沒有外推，但是，P的一次項(xiàng)參數(shù)為負(fù)數(shù)，屬于參數(shù)符號不合理的非典型式。5號試驗(yàn)點(diǎn)有最高產(chǎn)量點(diǎn)，N、P、K的一次項(xiàng)和二次項(xiàng)參數(shù)的符號滿足植物營養(yǎng)的一般規(guī)律，但鉀肥的最高施肥量超過試驗(yàn)設(shè)計施肥量，屬于推薦施肥結(jié)果遠(yuǎn)外推的非典型式。6—8號試驗(yàn)點(diǎn)的肥效模型存在最高產(chǎn)量點(diǎn)，模型一次項(xiàng)和二次項(xiàng)參數(shù)的符號正常，邊際產(chǎn)量導(dǎo)數(shù)法推薦施肥量在試驗(yàn)設(shè)計施肥量范圍內(nèi)，屬于典型肥效模型。

因此，采用目前常用的OLS回歸建模，1號和2號試驗(yàn)點(diǎn)屬于未達(dá)統(tǒng)計顯著水平的肥效模型，3—5號試驗(yàn)點(diǎn)建立的肥效模型屬于非典型式，只能棄之不用；只有 6—8號試驗(yàn)點(diǎn)的肥效模型滿足植物營養(yǎng)學(xué)的一般肥效規(guī)律，屬于三元典型肥效模型。

2.2 肥效模型的異方差及其統(tǒng)計檢驗(yàn)

OLS回歸建模要求模型隨機(jī)誤差項(xiàng)具有相同方差，否則，回歸建模結(jié)果就不可靠[16]，在應(yīng)用上易對肥效試驗(yàn)結(jié)果產(chǎn)生誤判。為此，需對表2的8個試驗(yàn)資料進(jìn)行回歸診斷。以各處理的試驗(yàn)產(chǎn)量為X軸，OLS擬合殘差為Y軸，繪制圖1的殘差散點(diǎn)圖。結(jié)果表明，無論是統(tǒng)計未達(dá)顯著水平的1號試驗(yàn)點(diǎn)肥效模型（圖1-a），還是3號試驗(yàn)點(diǎn)的非典型肥效模型（圖1-b）或者8號試驗(yàn)點(diǎn)的典型肥效模型（圖1-c），隨著水稻產(chǎn)量水平提高，其擬合殘差均呈現(xiàn)逐漸散開的漏斗狀分布趨勢，表明都不滿足隨機(jī)分布的假設(shè)。結(jié)果預(yù)示著14個處理產(chǎn)量的擬合殘差方差不相等，且隨著產(chǎn)量水平提高而增大，但1號試驗(yàn)點(diǎn)隨著稻谷產(chǎn)量水平增加，擬合殘差的增加絕對值最大，3號試驗(yàn)點(diǎn)則其次，8號試驗(yàn)點(diǎn)擬合殘差絕對值的變化幅度最小。

表3 三元二次多項(xiàng)式肥效模型的OLS回歸分析Table 3 OLS regression analysis of ternary quadratic polynomial’s fertilizer response model

圖1的擬合殘差圖給出了三元肥效模型異方差的定性分析結(jié)果，要回答異方差的嚴(yán)重程度，還需要做定量分析（表 4）。由于異方差在線性回歸模型中可能有不同的表現(xiàn)形式[19]，在不同的定量檢驗(yàn)方法中，只要有一種檢驗(yàn)結(jié)果達(dá)到統(tǒng)計顯著水平，就說明回歸模型存在該種形式的異方差。結(jié)果表明，擬合殘差絕對值變化幅度最小的8號試驗(yàn)點(diǎn)的OLS回歸模型異方差未達(dá)顯著水平，其余7個試驗(yàn)點(diǎn)建立的OLS回歸方程至少有一種異方差模型達(dá)到統(tǒng)計顯著水平，即存在嚴(yán)重的異方差。

2.3 可行廣義最小二乘法對肥效模型建模結(jié)果的影響

根據(jù)表4的異方差檢驗(yàn)結(jié)果，1—7號試驗(yàn)點(diǎn)肥效模型存在統(tǒng)計顯著水平的異方差，表3的OLS回歸建模結(jié)果就失去了可靠性。為此，采用可行廣義最小二乘法回歸建模（表5）。

圖1 三元肥效模型擬合殘差圖Fig. 1 Model fitting residual graph of ternary fertility efficiency

表4 肥效模型異方差的統(tǒng)計檢驗(yàn)Table 4 Statistical test of model’s heteroscedasticity

表5 可行廣義最小二乘法回歸建模結(jié)果Table 5 Modeling results using feasible generalized least squares regression

結(jié)果顯示，可行廣義最小二乘法明顯改善了肥效模型的建模效果。1號和2號試驗(yàn)點(diǎn)的F值和擬合優(yōu)度R2由OLS法的未達(dá)到顯著水平提高到極顯著水平，而且肥效模型轉(zhuǎn)化為典型式。其他6個試驗(yàn)點(diǎn)的肥效模型F值和擬合優(yōu)度R2同樣得到顯著的提升；3號、4號和5號試驗(yàn)點(diǎn)分別由OLS法的無最高產(chǎn)量點(diǎn)的、P項(xiàng)參數(shù)符號不合理、推薦施肥量外推的非典型式轉(zhuǎn)化為典型式。但是，OLS法屬于典型式的6號試驗(yàn)點(diǎn)，盡管模型存在異方差，若采用FGLS法則導(dǎo)致回歸模型由典型式變?yōu)榉堑湫褪?，結(jié)果反而變劣。對異方差未達(dá)顯著水平差異的8號試驗(yàn)點(diǎn)，采用FGLS回歸也能提高擬合優(yōu)度，模型同樣是典型式。

因此，在8個代表性試驗(yàn)點(diǎn)中，采用FGLS法回歸建模，顯著提高了肥效模型的擬合優(yōu)度和典型式的出現(xiàn)幾率，明顯地提高了建模成功率。

FGLS法為什么能提高建模成功率呢？比較1號、3號和8號試驗(yàn)點(diǎn)的兩種參數(shù)估計方法，F(xiàn)GLS回歸建模明顯地降低了肥效模型各個參數(shù)的方差（表 6）；表5和表3的均方誤差（MSE）結(jié)果也表明，F(xiàn)GLS回歸建模大幅度降低了誤差均方。結(jié)果說明FGLS大幅度提高了肥效模型的擬合精度和模型預(yù)測能力，從而提高了肥效模型的可靠性。

2.4 消除異方差后的肥效模型推薦施肥量

根據(jù)FGLS回歸建模得到的典型肥效模型，結(jié)合N、P2O5、K2O和稻谷價格分別為4.3、5.0、6.0和2.0元/kg的市場均價，應(yīng)用邊際產(chǎn)量導(dǎo)數(shù)法計算推薦施肥量（表7）。OLS回歸建模能得到典型肥效模型的7號和8號試驗(yàn)點(diǎn)，兩種建模方法所得肥效模型的推薦施肥量和預(yù)測產(chǎn)量水平十分相近；1—5號試驗(yàn)點(diǎn)，由OLS回歸建模未達(dá)顯著水平或是非典型肥效模型，用FGLS回歸建模轉(zhuǎn)化成了典型肥效模型，其推薦施肥量和預(yù)測產(chǎn)量水平也符合現(xiàn)有經(jīng)驗(yàn)認(rèn)識，尤其是，1號、3號和5號試驗(yàn)點(diǎn)的基礎(chǔ)土壤堿解氮含量分別為186、213、275 mg·kg-1（表 1），接近或超過早稻高產(chǎn)臨界指標(biāo)[17]，3個試驗(yàn)點(diǎn)的氮肥推薦用量明顯低于常規(guī)推薦水平。結(jié)果表明，F(xiàn)GLS回歸建模的推薦施肥量是可靠的。

表6 OLS與FGLS回歸的參數(shù)方差比較Table 6 Parameter variance compare between the regression methods of OLS and FGLS

表7 OLS肥效模型與FGLS肥效模型推薦施肥量的比較Table 7 Compare the recommended fertilization rate between the fertilizer response model of OLS and FGLS

2.5 肥效模型異方差出現(xiàn)幾率及其FGLS法效果評價

為了評價FGLS對三元肥效模型的建模效果，利用福建省早稻測土配方施肥項(xiàng)目在全省各地主要項(xiàng)目縣完成的171個“3414”設(shè)計方案的氮磷鉀田間肥效試驗(yàn)結(jié)果，逐個進(jìn)行回歸分析，統(tǒng)計結(jié)果見表8。FGLS大幅度提高了三元肥效模型的擬合優(yōu)度，每個試驗(yàn)點(diǎn)的R2和F值均達(dá)到統(tǒng)計顯著水平，而OLS則有16%試驗(yàn)點(diǎn)的肥效模型未達(dá)統(tǒng)計顯著水平。

模型典型性判別分析表明，在3種非典型肥效模型類型中，無論是OLS還是FGLS回歸建模，無最高產(chǎn)量點(diǎn)的非典型肥效模型所占比例高達(dá)40%左右，其次是參數(shù)符號不合理的非典型肥效模型所占試驗(yàn)點(diǎn)比例大致15%，而氮磷鉀推薦施肥量屬于外推的非典型肥效模型所占比例 6%左右，兩種方法的結(jié)果差異很?。籓LS回歸建模的典型肥效模型占試驗(yàn)點(diǎn)數(shù)的27.68%，而FGLS回歸建模的典型肥效模型比例達(dá)到34.99%，比OLS回歸建模提高了7.31個百分點(diǎn)。

表8 不同參數(shù)估計方法對肥效模型建模結(jié)果的影響Table 8 Effects of modeling results by using different methods of estimate parameters in fertilizer response model

3 討論

3.1 多項(xiàng)式肥效模型產(chǎn)生異方差的原因

在經(jīng)典回歸建模理論中，針對線性回歸方程隨機(jī)誤差項(xiàng)?u的一個重要假設(shè)是?u在不同觀測點(diǎn)具有同方差性。異方差將導(dǎo)致嚴(yán)重后果[16]，包括：（1）盡管OLS參數(shù)估計的無偏性仍然成立，并且基于此的產(chǎn)量預(yù)測也是無偏的，但參數(shù)估計值的方差不再是最小的，即最小二乘法失去了有效性，對大樣本也是如此。（2）由于異方差的影響，造成難以正確估計模型參數(shù)的標(biāo)準(zhǔn)誤，且無法辨別是正的還是負(fù)的偏差，導(dǎo)致回歸模型的t 統(tǒng)計量和F統(tǒng)計量失去了可靠性。（3）由于參數(shù)估計量不是有效的，從而對產(chǎn)量的預(yù)測也將不是有效的。二次多項(xiàng)式肥效模型是建立在經(jīng)典正態(tài)線性回歸模型理論基礎(chǔ)上的，要正確建立肥效模型就必須遵守回歸分析理論的基本假設(shè)[19]。長期以來，我們在建立作物肥效模型時，大都忽略了回歸分析理論有關(guān)同方差性等條件的基本假設(shè)，這是目前建模成功率低[4-6]的重要原因之一。

總結(jié)本研究結(jié)果并結(jié)合專業(yè)特點(diǎn)，多項(xiàng)式肥效模型產(chǎn)生異方差的主要原因有：（1）隨著試驗(yàn)設(shè)計和田間試驗(yàn)管理水平的提高，肥效模型的誤差方差σ2可能會減少。反之，若田間試驗(yàn)操作或管理出現(xiàn)失誤，可能導(dǎo)致部分處理出現(xiàn)異常產(chǎn)量結(jié)果，使模型誤差方差σ2出現(xiàn)非常數(shù)；（2）由于田間土壤肥力水平存在空間異質(zhì)性，可能導(dǎo)致部分試驗(yàn)處理出現(xiàn)異常產(chǎn)量值，導(dǎo)致出現(xiàn)異方差；（3）多項(xiàng)式肥效模型屬于經(jīng)驗(yàn)回歸模型，對某些試驗(yàn)點(diǎn)來說可能導(dǎo)致模型設(shè)定不正確；（4）在多數(shù)情況下，當(dāng)土壤養(yǎng)分肥力較低時，作物產(chǎn)量處于低產(chǎn)水平，此時通過施肥可迅速提高產(chǎn)量，而因施肥造成減產(chǎn)或平產(chǎn)的可能空間很小，即產(chǎn)量方差較?。坏诟弋a(chǎn)階段，通過施肥可能繼續(xù)增產(chǎn)，也可能不再增產(chǎn)，甚至可能因過量施肥造成減產(chǎn)，即產(chǎn)量變化的可能空間較大，此時方差必然較大。從這個意義上講，依據(jù)不同施肥量和產(chǎn)量建立的經(jīng)驗(yàn)肥效模型必然帶有異方差性質(zhì)。福建的171個早稻氮磷鉀肥效試驗(yàn)結(jié)果表明，這種異方差性達(dá)到統(tǒng)計顯著水平的試驗(yàn)點(diǎn)占到試驗(yàn)總數(shù)的25%左右（表8）。因此，通過控制試驗(yàn)田土壤空間異質(zhì)性、提高試驗(yàn)管理水平，降低試驗(yàn)產(chǎn)量異常值的出現(xiàn)幾率，是消減異方差和提高建模成功率的重要途徑。

3.2 多項(xiàng)式肥效模型異方差的檢驗(yàn)方法

針對異方差的嚴(yán)重危害性，在現(xiàn)代回歸分析理論中已經(jīng)研究提出了多種異方差定量診斷方法[19]。本研究選擇懷特（White）檢驗(yàn)、冠因克-巴塞特（KB）檢驗(yàn)、帕克（Park）檢驗(yàn)等3種常用的異方差定量檢驗(yàn)方法，探討三元肥效模型異方差性質(zhì)。結(jié)果表明，不同試驗(yàn)點(diǎn)建立的三元二次多項(xiàng)式肥效模型，誤差項(xiàng)的異方差有不同的表現(xiàn)形式。由于目前還難以明確在什么試驗(yàn)條件下會出現(xiàn)哪種形式的異方差，因而，在實(shí)用上，通常同時選用幾種不同異方差檢驗(yàn)方法進(jìn)行檢驗(yàn)。只要有一種檢驗(yàn)方法達(dá)到統(tǒng)計顯著水平，就說明肥效模型存在該種形式的異方差。

當(dāng)前，構(gòu)建異方差定量檢驗(yàn)方法的共同思路是，檢驗(yàn)隨機(jī)誤差項(xiàng)的方差與回歸模型自變量或因變量觀測值之間的相關(guān)性及其數(shù)學(xué)形式[19]，各種檢驗(yàn)方法正是在這個思路下發(fā)展起來的。本文選用3種常用的檢驗(yàn)方法，它們各有自己的特色和優(yōu)點(diǎn)。其中，懷特檢驗(yàn)不僅能夠檢驗(yàn)異方差的存在，同時在多變量的情況下，還能夠判斷出是哪一個變量引起的異方差。帕克檢驗(yàn)的優(yōu)點(diǎn)是不但能確定有無異方差性，而且還能給出異方差的具體函數(shù)形式。冠因克-巴塞特檢驗(yàn)的一個優(yōu)點(diǎn)是即使回歸模型中的誤差項(xiàng)不是正態(tài)分布的，它仍能適用。因此，這三種方法都是檢驗(yàn)肥效模型異方差的有效方法。

3.3 多項(xiàng)式肥效模型異方差的修正方法

用普通最小二乘法估計肥效模型參數(shù)時，對參數(shù)的估計是以產(chǎn)量擬合殘差平方和最小為條件的。在最小化過程中，對每個施肥處理的產(chǎn)量擬合殘差平方所提供的信息的重要程度是同等看待的，它們在決定肥效模型參數(shù)的過程中所起的作用是相同的，或者說取了相同的權(quán)數(shù)，即權(quán)數(shù)都為 1。這在同方差條件下，因不同施肥處理的作物產(chǎn)量水平在偏離均值程度上是相同的，這樣做是合理的。但是，在異方差條件下，由于各施肥處理的產(chǎn)量水平在偏離均值程度上相差很大，這樣做就未必合理。

本研究針對三元二次多項(xiàng)式肥效模型異方差問題，采用FGLS回歸建模。該法是對原模型加權(quán)，即：對較小 OLS擬合殘差平方的施肥處理產(chǎn)量賦予較大的權(quán)數(shù)，對較大OLS擬合殘差平方的施肥處理產(chǎn)量賦予較小的權(quán)數(shù)，使回歸模型轉(zhuǎn)化成不存在異方差性的新模型，然后再次采用OLS估計其參數(shù)。結(jié)果表明，與OLS回歸建模方法相比，F(xiàn)GLS方法能將絕大部分試驗(yàn)點(diǎn)的試驗(yàn)結(jié)果由 OLS方法的未達(dá)統(tǒng)計顯著水平轉(zhuǎn)化為達(dá)到統(tǒng)計顯著水平，或者將部分試驗(yàn)點(diǎn)的非典型式模型轉(zhuǎn)化為典型肥效模型（表3和表5），明顯提高了建模成功率。

在對付線性回歸模型異方差上，除了FGLS外，“OLS+穩(wěn)健標(biāo)準(zhǔn)誤”也是常用方法[22]。該法是在異方差情況下，仍然采用OLS，但對OLS估計量的標(biāo)準(zhǔn)差進(jìn)行修正。與OLS估計比較，肥效模型參數(shù)估計值沒有變化，但是參數(shù)估計值的方差和標(biāo)準(zhǔn)差變化明顯。經(jīng)過這種修正的模型參數(shù)標(biāo)準(zhǔn)差，模型參數(shù)的顯著性檢驗(yàn)和模型產(chǎn)量預(yù)測等方面就變?yōu)橛行У?。那么，究竟是使用“OLS+穩(wěn)健標(biāo)準(zhǔn)誤”還是FGLS？從肥效模型角度看，我們總是希望將OLS建模出現(xiàn)的非典型模型盡可能多地轉(zhuǎn)化為典型肥效模型。由于“OLS+穩(wěn)健標(biāo)準(zhǔn)誤”方法的模型參數(shù)沒有變化，顯然不能實(shí)現(xiàn)這一目標(biāo)，而FGLS方法則能把部分非典型肥效模型轉(zhuǎn)化為典型式。從這個意義上講，F(xiàn)GLS方法具有優(yōu)勢。同時，表7的建模結(jié)果使我們注意到，OLS回歸建模是典型肥效模型的試驗(yàn)資料，若采用FGLS回歸建模，可能出現(xiàn)非典型式，使結(jié)果變劣。因此，在實(shí)際工作中，我們也不應(yīng)對異方差反應(yīng)過度[19]。如果OLS方法能得到典型肥效模型，就無需考慮其他修正方法。

福建 171個早稻氮磷鉀肥效試驗(yàn)資料的總結(jié)結(jié)果表明，F(xiàn)GLS回歸建模的三元典型肥效模型比例僅達(dá)到34.99%，只比OLS回歸建模提高7.31個百分點(diǎn)。究其原因，多項(xiàng)式肥效模型不僅存在異方差，而且還存在嚴(yán)重的多重共線性[23]問題，制約了OLS法的有效性[19]，對FGLS回歸建模也不例外。只有較好地克服了這兩個建模問題，才可能較大幅度地提高建模成功率。如何同時克服多項(xiàng)式肥效模型的異方差和多重共線性問題，提出最優(yōu)建模策略，則有待于未來的進(jìn)一步研究。

4 結(jié)論

根據(jù)171個早稻“3414”設(shè)計的氮磷鉀田間肥效試驗(yàn)結(jié)果，現(xiàn)代回歸分析理論的可行廣義最小二乘法建模技術(shù)是消除三元二次多項(xiàng)式肥效模型異方差的有效方法。該法明顯減小模型各個參數(shù)的方差，改善了肥效模型的擬合精度和模型預(yù)測能力，進(jìn)而提高了田間肥效試驗(yàn)資料的建模成功率。

[1] VALKAMA E, UUSITALO R, TURTOLA E. Yield response models to phosphorus application: a research synthesis of Finnish field trials to optimize fertilizer P use of cereals. Nutrient Cycling in Agroecosystems,2011, 91(1): 1-15.

[2] GRWGORT M C, ZORITA M D, DARDANELLI J, Bongiovanni R G.Regional model for nitrogen fertilization of site-specific rainfed corn in haplustolls of the central Pampas, Argentina. Precision Agriculture,2011, 12(6): 831-849.

[3] 毛達(dá)如, 張承東. 推薦施肥技術(shù)中施肥模型與試驗(yàn)設(shè)計的研究. 土壤通報, 1991, 22(5): 216-218.MAO D R , ZHANG C D. Fertilization model and the experiment design of fertilization recommendation technology. Chinese Journa1 of Soil Science, 1991, 22(5): 216-218. (in Chinese)

[4] 王興仁. 二元二次肥料效應(yīng)曲面等產(chǎn)線圖在科學(xué)施肥中的位置（一）. 土壤通報, 1985, 16(1): 30-34.WANG X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization (part 1).Chinese Journa1 of Soil Science, 1985, 16(1): 30-34. (in Chinese)

[5] 王興仁. 二元二次肥料效應(yīng)曲面等產(chǎn)線圖在科學(xué)施肥中的位置（二）. 土壤通報, 1985, 16(2): 86-88.WANG X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization (part 2).Chinese Journa1 of Soil Science, 1985, 16(2): 86-88. (in Chinese)

[6] 章明清, 林仁塤, 林代炎, 姜永. 極值判別分析在三元肥效模型推薦施肥中的作用. 福建農(nóng)業(yè)學(xué)報, 1995, 10(2): 54-59.ZHANG M Q, LIN R X, LIN D Y, JIANG Y. Function of distinguish analysis on extreme value in recommendatory fertilization for three-fertilizer efficiency model. Fujian Journal of Agricultural Sciences, 1995, 10(2): 54-59. (in Chinese)

[7] CERRATO M E, BLACKMER A M. Comparison of models for describing corn yield response to nitrogen fertilizer. Agronomy Journal, 1990, 82(1): 138-143.

[8] 楊靖一, WADSWORTH G A, GREENWOOD D J. 三種蔬菜N肥效應(yīng)曲線的比較研究. 植物營養(yǎng)與肥料學(xué)報, 1995, 1(1): 71-78.YANG J Y, WADSWORTH G A, Greenwood D J. Comparative study on three kinds of vegetables response to N fertilizer curve. Plant Nutrition and Fertilizer Science, 1995, 1(1): 71-78. (in Chinese)

[9] 陳新平, 周金池, 王興仁, 張福鎖, 寶德俊. 小麥-玉米輪作制中氮肥效應(yīng)模型的選擇——經(jīng)濟(jì)和環(huán)境效益分析. 土壤學(xué)報, 2000,37(3): 346-354.CHEN X P, ZHOU J C, WANG X R, ZHANG F S, BAO D J.Economic and environmental evaluation on models for describing crop yield response to nitrogen fertilizers at winter-wheat and summer-corn rotation system. Acta Pedologica Sinica, 2000, 37(3):346-354. (in Chinese)

[10] 吳平, 陶勤南. 氮磷鉀肥效模型的研究及其應(yīng)用. 浙江農(nóng)業(yè)大學(xué)學(xué)報, 1989, 15(4): 383-388.WU P, TAO Q N. Study on N, P and K fertilizer response model and its application. Journal of Zhejiang Agricultural University, 1989,15(4): 383-388. (in Chinese)

[11] 章明清. 米氏肥效方程的灰色建模法. 土壤通報, 1997, 28(1):18-19.ZHANG M Q. Gray modeling method of Mitsherlich fertilizer response equations. Chinese Journa1 of Soil Science, 1997, 28(1):18-19. (in Chinese)

[12] 章明清, 徐志平, 姚寶全, 林瓊, 顏明娟. Monte Carlo 法在多元肥效模型參數(shù)估計和推薦施肥中的應(yīng)用. 植物營養(yǎng)與肥料學(xué)報, 2009,15(2): 366-373.ZHANG M Q, XU Z P, YAO B Q, LIN Q, YAN M J. Using Monte Carlo method for parameter estimation and fertilization recommendation of multivariate fertilizer response model. Plant Nutrition and Fertilizing Science, 2009, 15(2): 366-373. (in Chinese)

[13] COLWELL J D. The erivation of fertilizer recommendations for crops in non-uniform environment. Crop Quality and Economy, 1974:936-961.

[14] 王興仁, 陳倫壽, 毛達(dá)如, 張才東. 分類回歸綜合法及其在區(qū)域施肥決策中的應(yīng)用. 土壤通報, 1989, 20(1): 17-21.WANG X R, CHEN L S, MAO D R, ZHANG C D. Classification of regression synthesis and its application for regional fertilization decision-making. Chinese Journa1 of Soil Science, 1989, 20(1): 17-21.(in Chinese)

[15] 毛達(dá)如, 張承東. 多點(diǎn)肥料效應(yīng)函數(shù)的動態(tài)聚類方法. 北京農(nóng)業(yè)大學(xué)學(xué)報, 1991, 17(2): 49-54.MAO D R, ZHANG C D. Cluster analysis of quadratic response function on the fertilizer dispersed experiment. Acta Agriculturae Universitatis Pekinensis, 1991, 17(2): 49-54. (in Chinese)

[16] MONTGOMERY D C, PECK E A, VINING G G. 線性回歸分析導(dǎo)論. 第5版. 王辰勇, 譯. 北京: 機(jī)械工業(yè)出版社, 2016: 91-268.MONTGOMERY D C, PECK E A, VINING G G. Introduction to Linear Regression Analysis. 5th ed. WANG C Y, trans. Beijing: China Machine Press, 2016: 91-268. (in Chinese)

[17] 李娟, 章明清, 孔慶波, 姚寶全, 顏明娟. 福建早稻測土配方施肥指標(biāo)體系研究. 植物營養(yǎng)與肥料學(xué)報, 2010, 16(4): 938-946.LI J, ZHANG M Q, KONG Q B, YAO B Q, YAN M J. Soil testing and formula fertilization index for early rice in Fujian Province. Plant Nutrition and Fertilizer Science, 2010, 16(4): 938- 946. (in Chinese)

[18] 魯如坤. 土壤農(nóng)業(yè)化學(xué)分析方法. 北京: 中國農(nóng)業(yè)科技出版社,2000: 146-196.LU R K. Soil Agricultural Chemical Analysis Method. Beijing: China's Agricultural Science and Technology Press, 2000: 146-196. (in Chinese)

[19] GUJARATI D N, PORTER D C. 計量經(jīng)濟(jì)學(xué)基礎(chǔ) (上冊). 第5版.費(fèi)劍平, 譯. 北京: 中國人民大學(xué)出版社, 2011: 313-399.GUJARATI D N, PORTER D C. Econometrics Foundation. 5th ed.FEI J P, trans. Beijing: China Renmin University Press, 2011: 313-399.(in Chinese)

[20] 尹光霞. 多元線性回歸模型中的異方差性問題. 湖北大學(xué)學(xué)報(自然科學(xué)版), 2003, 25(2): 121-125.YIN G X. Problems on heteroscedasticity in multilinear regression models. Journal of Hubei University (Natural Science Edition), 2003,25(2): 121-125. (in Chinese)

[21] 何曉群，劉文卿. 應(yīng)用回歸分析. 3版. 北京: 中國人民大學(xué)出版社,2013: 158-170.HE X Q, LIU W Q. Applied Regression Analysis. 3rd ed. Beijing:China Renmin University Press, 2013: 158-170. (in Chinese)

[22] 陳強(qiáng). 高級計量經(jīng)濟(jì)學(xué)及Stata應(yīng)用. 2版. 北京: 高等教育出版社,2015: 87-99.CHEN Q. Advanced Econometric and Stata Applications.2nd ed.Beijing: Higher Education Press, 2015: 87-99. (in Chinese)

[23] 章明清, 李娟, 孔慶波, 嚴(yán)芳. 作物肥料效應(yīng)函數(shù)模型研究進(jìn)展與展望. 土壤學(xué)報, 2016, 53(6): 1143-1156.ZHANG M Q, LI J, KONG Q B, YAN F. Progress and prospect of the study on crop-response-to -fertilization function model. Acta Pedologica Sinica, 2016, 53(6): 1143-1156. (in Chinese)