基于混合效應(yīng)的興安落葉松單木斷面積生長模型研究

覃鑫浩

(國家林業(yè)和草原局 調(diào)查規(guī)劃設(shè)計(jì)院，北京 100714)

準(zhǔn)確描述林分或林木的生長規(guī)律，從而更好地為森林的經(jīng)營與利用提供參考依據(jù)。描述生長規(guī)律主要包括前人積累的定性經(jīng)驗(yàn)以及林分生長與收獲模型[1]。根據(jù)模擬目的和對象的不同將森林生長和收獲模型劃分為林分-徑階-單木3個(gè)水平。單木生長模型的模擬對象是單株樹木，將林木大小，競爭以及立地條件表述為單木生長量的數(shù)學(xué)表達(dá)式，單木模型能夠描述復(fù)雜的林分關(guān)系，不僅能適用于人工純林，其對異齡混交林和天然林也具有較強(qiáng)的適應(yīng)性[2-4]。此外，與徑階模型和全林分模型相比，單木模型針對不同經(jīng)營措施模擬得到的林分生長動(dòng)態(tài)更準(zhǔn)確，為森林質(zhì)量的精準(zhǔn)提升提供理論依據(jù)[5]。國內(nèi)學(xué)者對不同樹種的單木生長模型做了大量研究，邵國凡以紅松人工林為研究對象，構(gòu)建了紅松(Pinuskoraiensis)的單木生長模型，為單木模型在我國的普及和應(yīng)用奠定基礎(chǔ)[6]；孟憲宇等[7]采用生長量修正法構(gòu)建了華北落葉松(Larixprincipis-rupprechtii)人工林與距離無關(guān)的單木生長模型；呂理興[8]以栓皮櫟(Quercusvariabilis)天然林為研究對象，選用簡單競爭指數(shù)并建立以個(gè)體樹木生長為基礎(chǔ)的、與距離有關(guān)的單木生長模型。

傳統(tǒng)建模方法為最小二乘法，其假設(shè)數(shù)據(jù)獨(dú)立、等方差、正態(tài)分布等特點(diǎn)，然而林分生長和收獲模擬數(shù)據(jù)為連續(xù)多次重復(fù)觀測數(shù)據(jù)，數(shù)據(jù)之間普遍存在異方差和自相關(guān)，因而不滿足普通回歸分析中的獨(dú)立性假設(shè)[9-10]?；旌闲?yīng)模型能夠降低分組數(shù)據(jù)的異方差性和自相關(guān)性，從而提高預(yù)測精度并解釋隨機(jī)誤差的來源[11-12]。雷相東等[13]以天然落葉松云冷杉林為研究對象，研究發(fā)現(xiàn)考慮層次結(jié)構(gòu)的混合效應(yīng)模型的擬合效果優(yōu)于傳統(tǒng)模型，誤差、均方根誤差及其相對值均顯著地減少；杜志等[14]考慮樣地效應(yīng)的馬尾松(Pinusmassoniana)天然次生林單木胸高斷面積生長模型優(yōu)于傳統(tǒng)模型；龔召松等[15]考慮樣地效應(yīng)的楠木(Phoebezhennan)次生林林分?jǐn)嗝娣e與蓄積相容性混合效應(yīng)生長收獲預(yù)估模型在精度、平均誤差和均方根誤差都優(yōu)于傳統(tǒng)模型。

興安落葉松(Larixgmelinii)是我國最為典型的東部高寒地區(qū)的優(yōu)勢樹種，喜光性強(qiáng)，對水分要求較高，在內(nèi)蒙古大興安嶺地區(qū)分布廣泛，是該地區(qū)的地帶性植被[16]，同時(shí)他也是我國經(jīng)濟(jì)建設(shè)主要用材樹種[17-18]。以興安落葉松為主要建群樹種的森林對內(nèi)蒙古大興安嶺林區(qū)的水源涵養(yǎng)、空氣凈化等功能發(fā)揮著重要作用?？茖W(xué)地經(jīng)營興安落葉松林才能更好地發(fā)揮大興安嶺林區(qū)的生態(tài)功能和價(jià)值，在經(jīng)營目標(biāo)轉(zhuǎn)變和森林功能多元化影響下，及時(shí)準(zhǔn)確地掌握森林資源現(xiàn)狀和動(dòng)態(tài)，開展以森林提升為目標(biāo)的興安落葉松林密度調(diào)控技術(shù)，提高內(nèi)蒙古大興安嶺林區(qū)森林資源的現(xiàn)代化經(jīng)營管理水平。

1 研究區(qū)概況

內(nèi)蒙古大興安嶺林區(qū)為大興安嶺的西半部分，南北長696 km，東西寬384 km，地理坐標(biāo)119°36′20″-125°20′50″E，46°08′40″-53°20′00″E，海拔425～1 760 m[19]。地處歐亞大陸中高緯度地區(qū)，氣候類型為寒溫帶的大陸性季風(fēng)氣候，春季干旱且多風(fēng)，夏季炎熱多雨，年平均氣溫-2～4 ℃，年均降水量450 mm。土壤類型主要包括漂灰土、黑土、暗棕壤、灰灰黑土和黑鈣土[20]。主要喬木樹種是興安落葉松、白樺(Betulaplatyphylla)、山楊(Populusdavidiana)和櫟類等，植被類型主體屬于歐亞針葉林植物區(qū)，以中、低山森林植被占優(yōu)勢，共有維管束植物1 339種，其中種子植物1 297種，蕨類植物42種。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)

采用第六次至第九次內(nèi)蒙古大興安嶺重點(diǎn)國有林區(qū)森林資源連續(xù)清查數(shù)據(jù)，以興安落葉松天然林為研究對象，以單水平線性混合效應(yīng)模型的方法，構(gòu)建單木胸高斷面積生長的混合效應(yīng)模型。樣地面積0.067 hm2。外業(yè)調(diào)查因子包括林分調(diào)查因子(樣地的地理坐標(biāo)(經(jīng)緯度)、坡向、坡位、坡地、地貌、起源、齡組、優(yōu)勢樹種、平均年齡、平均胸徑、林種和起源等)和樣木調(diào)查因子(樹種名稱、林木的胸徑)。樣地分布如圖1所示。

圖1 興安落葉松天然林樣地分布Fig.1 Distribution map of Larix gmelinii forest sample plots

2.2 研究方法

2.2.1 變量選擇及基礎(chǔ)模型的構(gòu)建 林木生長的主要影響因子有林木的大小、林木所處的立地條件以及林木之間的競爭等，根據(jù)前人研究經(jīng)驗(yàn)，本研究將林木生長描述為三者的線性函數(shù)[2-3]。采用林木大小、競爭及立地條件來描述單木胸高斷面積的生長。其中林木大小為林木的胸徑以及胸徑的變形形式；林木競爭是林分生長必不可少的因素，對林木的生長趨勢有重要的影響，本研究選擇與距離無關(guān)的競爭指數(shù)參與模型構(gòu)建，與距離無關(guān)的競爭指數(shù)主要包含以下4個(gè)變量，即每公頃株數(shù)(N,trees·hm-2),每公頃斷面積(BA,m2·hm-2)，相對密度指數(shù)RD(對象木胸徑與林分平均胸徑之比：DBH/QMD)以及比對象木胸徑大的斷面積之和(BAL,m2·hm-2)；立地條件是影響林分生長的重要因子，本研究主要考慮海拔(EL)，坡度(SL)以及坡向(ASP)對于單木生長的影響。其中，采用Stage提出的公式對ASP和SL進(jìn)行組合[1,13,21]，即SLcos=tan(SL)×cos(ASP)，SLsin=tan(SL)×sin(ASP)。

因變量的形式會(huì)對模型的表現(xiàn)產(chǎn)生很大的影響，在模型模擬前，分別采用不同的因變量對模型進(jìn)行預(yù)測，根據(jù)殘差圖和評價(jià)指標(biāo)，當(dāng)自變量為ln(DBH22-DBH2)時(shí)，模型的模擬效果最好，因此選擇ln(DBH22-DBH2)作為模型的因變量。

利用多元逐步回歸方法，建立傳統(tǒng)的興安落葉松基礎(chǔ)模型，根據(jù)樹木的生長規(guī)律，選擇對興安落葉松生長影響性強(qiáng)的變量進(jìn)入模型，并對模擬的結(jié)果采用多重共線性檢驗(yàn)，最終選擇對興安落葉松生長影響顯著，即回歸系數(shù)顯著(P<0.05)且方差膨脹因子(VIF)<5的變量進(jìn)入模型。

篩選出興安落葉松樣地626塊，樣木株數(shù)60 279株。其中501塊樣地?cái)?shù)據(jù)，興安落葉松47 524株，用于構(gòu)建興安落葉松的胸高斷面積預(yù)估模型，采用剩余的125塊樣地，興安落葉松12 755株對構(gòu)建的模型進(jìn)行檢驗(yàn)。建模數(shù)據(jù)統(tǒng)計(jì)如表1所示。

表1 興安落葉松建模數(shù)據(jù)變量統(tǒng)計(jì)Table 1 Statistical table of L.gmelinii modeling data variables

由表1可以看出，興安落葉松建模數(shù)據(jù)期初最小胸徑為5 cm，最大胸徑為74.1 cm；期末最小胸徑為5 cm，最大胸徑為79 cm；胸徑的增量從0.1～5.0 cm；胸高斷面積之和最大為32.08，最小為0.16，BAL最大值為2.14 m2，最小為0，林分每公頃株數(shù)最大為3 776，最小為15；坡度正切值變化為0.0～0.7，樣地海拔分布266～1 650 m。

2.2.2 混合效應(yīng)模型的構(gòu)建 采用單水平線性混合效應(yīng)模型的方法構(gòu)建內(nèi)蒙古大興安嶺重點(diǎn)國有林區(qū)興安落葉松的單木胸高斷面積生長預(yù)估模型，模型的一般形式如下[22]

(1)

式中:yi為第i個(gè)樣地的單木斷面積生長量；Xi為固定效應(yīng)即與林木大小、競爭和立地條件相關(guān)的設(shè)計(jì)矩陣；β為由固定效應(yīng)參數(shù)構(gòu)成的向量；Zi為隨機(jī)效應(yīng)設(shè)計(jì)矩陣；bi為模型隨機(jī)參數(shù)構(gòu)成的向量；ei為誤差向量；n為樣地個(gè)數(shù)；D為隨機(jī)效應(yīng)的方差-協(xié)方差矩陣；σ2為模型方差；Ri為誤差方差-協(xié)方差結(jié)構(gòu)。

混合效應(yīng)模型既有固定效應(yīng)參數(shù)又有隨機(jī)效應(yīng)參數(shù)，其中隨機(jī)效應(yīng)包括隨機(jī)截距模型、隨機(jī)截距和斜率模型。本研究將基礎(chǔ)模型的參數(shù)均作為隨機(jī)效應(yīng)參數(shù)，需計(jì)算隨機(jī)截距和斜率模型，根據(jù)計(jì)算結(jié)果對不同的隨機(jī)效應(yīng)參數(shù)進(jìn)行組合，并對不同組合進(jìn)行模擬，最后篩選出最優(yōu)的參數(shù)組合。

在混合效應(yīng)的參數(shù)估算過程中，為了確定合適的方差協(xié)方差結(jié)構(gòu)，對復(fù)合對稱(CS)、一階自回歸[AR(1)]與一階自回歸移動(dòng)平均結(jié)構(gòu)[ARMA(1,1)]進(jìn)行檢驗(yàn)。采用指數(shù)函數(shù)(Exponent)、冪函數(shù)(Power)和常數(shù)加冪函數(shù)(ConstPower)對方差異質(zhì)性進(jìn)行描述。

2.2.3 模型的選擇與檢驗(yàn)

1)模型的篩選依據(jù)AIC(AIC)、BIC(BIC)和-2LL(-2Loglik,-2Loglik)越小越好原則，同時(shí)要對篩選出的含有不同隨機(jī)參數(shù)個(gè)數(shù)的混合模型進(jìn)行似然比檢驗(yàn)(LRT)避免過多參數(shù)化問題的產(chǎn)生。

AIC=-2Loglik+2d

(2)

BIC=-2Loglik+dln(n)

(3)

LRT=2[log(L2)-log(L1)]

(4)

2)計(jì)算模型的隨機(jī)效應(yīng)參數(shù)值。

(5)

式中:D為隨機(jī)效應(yīng)方差-協(xié)方差矩陣；ZK為設(shè)計(jì)矩陣；RK為誤差方差-協(xié)方差矩陣；ek為誤差項(xiàng)。

3)采用平均絕對誤差(Bias)、均方根誤差(RMSE)和決定系數(shù)(R2)對模型的擬合效果進(jìn)行評價(jià)。

2.2.4 統(tǒng)計(jì)處理 樣地?cái)?shù)據(jù)的篩選與處理在Excel軟件中完成，基礎(chǔ)模型和混合效應(yīng)模型的構(gòu)建以及模型的檢驗(yàn)在R語言軟件中完成。

3 結(jié)果與分析

3.1 興安落葉松基礎(chǔ)模型結(jié)果

基于內(nèi)蒙古大興安嶺重點(diǎn)國有林區(qū)第六次至第九次4期森林資源連續(xù)清查數(shù)據(jù)，篩選出興安落葉松樣木株數(shù)47524株用于基礎(chǔ)模型的構(gòu)建，得到最終的預(yù)估模型

ln(DBH22-DBH2)=β1+β21/DBH+β3NT+β4BAL+β5SLcos+eij

(6)

式中:DBH表示期初胸徑，DBH2表示期末胸徑，β1、β2、β3、β4和β5表示固定效應(yīng)參數(shù)，NT表示每公頃株數(shù)，BAL表示樣地內(nèi)比對象木大的所有胸高斷面積之和，即SLcos=SL×cos(ASP)，eij表示模型誤差。

模型模擬結(jié)果如表2所示。從表2可以看出，1/DBH、NT、BAL和SLcos對興安落葉松單木斷面積的生長有顯著影響。興安落葉松斷面積生長模型的決定系數(shù)R2為0.398。根據(jù)構(gòu)建的模型使用R(ggplot2)繪制興安落葉松的殘差圖(圖2)，從殘差圖發(fā)現(xiàn)模型存在異方差性。

表2 傳統(tǒng)方法模擬興安落葉松單木胸高斷面積預(yù)估模型的結(jié)果Table 2 The results of L.gmelinii basal area growth model based on traditional regression

圖2 基礎(chǔ)模型的殘差圖和QQ圖Fig.2 Residual plot and quantile-quantile (QQ) plot of the basic model

3.2 混合效應(yīng)模型

在基礎(chǔ)模型的基礎(chǔ)上，引入樣地效應(yīng)，通過對常數(shù)項(xiàng)(int)、1/DBH、NT、BAL和SLcos 5個(gè)變量不同隨機(jī)參數(shù)組合進(jìn)行擬合，構(gòu)建興安落葉松斷面積生長混合效應(yīng)模型，擬合收斂的情況共有31種(表3)。通過比較擬合的AIC、BIC以及-2LL發(fā)現(xiàn)，相比于傳統(tǒng)模型，基于樣地效應(yīng)的混合效應(yīng)模型的擬合效果更好。其中，只含有1個(gè)隨機(jī)效應(yīng)參數(shù)的擬合方程，常數(shù)項(xiàng)(模擬2)的擬合效果最好；當(dāng)擬合方程含有2個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，常數(shù)項(xiàng)和1/DBH的組合(模擬7)擬合效果最好；當(dāng)包含3個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，常數(shù)項(xiàng)、1/DBH和BAL的組合(模擬18)擬合效果最好；當(dāng)含有4個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，常數(shù)項(xiàng)、1/DBH、BAL和SLcos的組合(模擬29)擬合效果最好；當(dāng)含有常數(shù)項(xiàng)、1/DBH、NT、BAL和SLcos共5個(gè)隨機(jī)效應(yīng)參數(shù)時(shí)，得到模擬32。當(dāng)擬合方程包含常數(shù)項(xiàng)、1/DBH、NT、BAL和SLcos(模擬32)的擬合效果優(yōu)于其他任意一個(gè)組合的擬合效果，模擬精度最高。

對于包含不同參數(shù)的5種最優(yōu)模型，似然比檢驗(yàn)結(jié)果(表3)發(fā)現(xiàn)，模擬32>模擬29>模擬18>模擬7>模擬2，隨著隨機(jī)效應(yīng)參數(shù)個(gè)數(shù)的增加，模型的模擬效果越好。最終選擇模擬32(式7)作為最佳參數(shù)組合的混合效應(yīng)模型。

表3 胸高斷面積生長量線性混合模型收斂結(jié)果Table 3 The lack-of-fit statistics of the linear mixed effects individual-tree basal area growth model

ln(DBH22-DBH2)=(β1+b1)+(β2+b2)1/DBH+(β3+b3)NT+(β4+b4)BAL+(β5+b5)SLcos+eij

(7)

式中:β1～β5是固定效應(yīng)參數(shù)；b1～b4是隨機(jī)效應(yīng)參數(shù)。

由于數(shù)據(jù)存在異方差性，本研究用指數(shù)函數(shù)、冪函數(shù)和常數(shù)加冪函數(shù)作為異方差結(jié)構(gòu)消除數(shù)據(jù)的異方差性。表4顯示，混合效應(yīng)模型的3個(gè)異方差結(jié)構(gòu)均顯著改善了AIC值、BIC值和Loglik值，提高了模型的預(yù)估能力。從表4計(jì)算的LRT結(jié)果發(fā)現(xiàn)，指數(shù)函數(shù)在消除數(shù)據(jù)的異方差性中表現(xiàn)最佳，作為該模型最終的異方差函數(shù)。此外，由于數(shù)據(jù)為多期重復(fù)測量數(shù)據(jù)，因此存在自相關(guān)性。采用3種常用的自相關(guān)結(jié)構(gòu)對數(shù)據(jù)的相關(guān)性進(jìn)行表達(dá)，結(jié)果表明除ARMA(1,1)未收斂外，其余2種自相關(guān)結(jié)構(gòu)均改善了模型的表型，其中AR(1)表現(xiàn)最好，確定為模型最優(yōu)自相關(guān)結(jié)構(gòu)。

表4 考慮異方差函數(shù)和自相關(guān)結(jié)構(gòu)后混合效應(yīng)模型結(jié)果比較Table 4 The lack-of-fit statistics of the linear mixed effects individual-tree basal area growth model using different error variance-covariance structures

確定了參數(shù)效應(yīng)和誤差方差-協(xié)方差結(jié)構(gòu)后，基于限制性極大似然法(REML)，得到最終的興安落葉松單木胸高斷面積生長混合效應(yīng)模型

ln(DBH22-DBH2)=(4.477 5+b1)+(-11.440 8+b2)1/DBH+(-0.000 2+b3)NT+(0.749 2+b4)BAL+(0.066 6+b5)SLcos+eij

(8)

圖3為混合效應(yīng)模型的殘差圖和QQ圖。通過對比發(fā)現(xiàn)，模型殘差圖在引入隨機(jī)效應(yīng)之后得到了改善。

圖3 混合效應(yīng)模型的殘差圖與QQ圖Fig.3 Residual plot and QQ plot of our final mixed effects model

3.3 模型檢驗(yàn)

采用剩余的125塊樣地、12 755株興安落葉松對構(gòu)建的混合模型與傳統(tǒng)模型進(jìn)行檢驗(yàn)，比較模型的均方根誤差、平均絕對誤差和決定系數(shù)?；A(chǔ)模型的檢驗(yàn)只需要考慮固定效應(yīng)，將自變量參數(shù)代入到模型中，與實(shí)測值進(jìn)行比較即可；考慮到混合效應(yīng)模型的固定效應(yīng)和樣地的隨機(jī)效應(yīng)，本研究將獲得的固定效應(yīng)參數(shù)與隨機(jī)效應(yīng)參數(shù)求和，帶入模型得到擬合值并與實(shí)測值進(jìn)行比較，判斷模型的模擬效果。結(jié)果如表5和圖4所示，傳統(tǒng)最小二乘方法的3個(gè)統(tǒng)計(jì)值計(jì)算結(jié)果分別為平均絕對誤差0.671 9，均方根誤差0.820 7，模型的決定系數(shù)R2為0.353 7?？紤]樣地效應(yīng)的混合效應(yīng)模型的平均絕對誤差為0.588 9，均方根誤差為0.730 0，模型的決定系數(shù)R2為0.488 6。結(jié)果顯示考慮樣地效應(yīng)的混合效應(yīng)模型的平均絕對殘差與均方根誤差都比基礎(chǔ)模型的低，模型的決定系數(shù)R2由0.353 7提高到0.488 6。模型檢驗(yàn)結(jié)果證明考慮樣地效應(yīng)的混合效應(yīng)模型的模擬效果優(yōu)于一般線性模型。

圖4 興安落葉松胸高斷面積基礎(chǔ)模型與混合效應(yīng)模型觀測值與預(yù)測值Fig.4 The observed and predicted values of the basic model and mixed effects model for L.gmelinii

表5 傳統(tǒng)模型與混合效應(yīng)模型驗(yàn)證擬合統(tǒng)計(jì)值比較Table 5 The lack-of-fit statistics of model validation using the basic model and mixed effects model

4 結(jié)論與討論

以內(nèi)蒙古大興安嶺地區(qū)興安落葉松天然林為研究對象，構(gòu)建了興安落葉松天然林中興安落葉松的胸高斷面積生長模型。林木的大小即胸徑的倒數(shù)(1/DBH)、每公頃株數(shù)(NT)、比對象木大的胸高斷面積之和(BAL)、坡度和坡向的組合(SLcos)4個(gè)參數(shù)明顯影響了興安落葉松天然林的生長效果，其中林分生長與胸徑的倒數(shù)呈負(fù)相關(guān)，即林木胸徑隨個(gè)體胸徑倒數(shù)的減小而增大，因此提升了林木獲取資源的能力，促使林木生長量增大。杜志等[14]以湖南地區(qū)馬尾松為研究對象，王建軍等[1]以福建省三明市杉木為研究對象，Zhao等[23]以硬闊樹種為研究對象，R.Cannell等[24]以云杉(Piceaasperata)和松樹為研究對象構(gòu)建各個(gè)樹種的生長模型發(fā)現(xiàn)了相同的規(guī)律。林分密度(NT)表示單位面積林木株數(shù)，它能間接反映單位面積林木間的競爭，胸高斷面積生長量與每公頃株數(shù)成反比，單位面積林木株數(shù)越多，獲取光和水分的幾率就越小，樹木間的競爭壓力越大，胸高斷面積生長量越小。BAL是比對象木大的所有林木的胸高斷面積之和，BAL值越大則表明林分中比對象木大的樹木越多，對象木受到的競爭壓力越大，因此與胸高斷面積生長量成反比。坡度和坡向的組合(SLcos)對林分的生長表現(xiàn)顯著。由于森林監(jiān)測數(shù)據(jù)具有連續(xù)重讀觀測的特點(diǎn)，與傳統(tǒng)方法會(huì)得到有偏估計(jì)，本研究通過對不同的隨機(jī)效應(yīng)參數(shù)進(jìn)行組合，篩選最優(yōu)的參數(shù)組合和誤差方差-協(xié)方差結(jié)構(gòu)提高了模型預(yù)測精度。與傳統(tǒng)模型相比混合效應(yīng)的決定系數(shù)從0.353 7提高到0.488 6，誤差、均方根誤差及其相對值均顯著減少，建立的模型有一定的實(shí)用性。本研究以興安落葉松天然林為研究對象，探究主要樹種胸高斷面積生長模型，發(fā)現(xiàn)林木大小、林地立地條件以及林木之間的競爭顯著影響著林木的生長，這個(gè)與很多專家的研究結(jié)果一致，并且上述因變量在林業(yè)調(diào)查工作中易于獲取。因此，本研究的興安落葉松胸高斷面積生長模型具有一定的生物學(xué)意義，統(tǒng)計(jì)結(jié)果可靠。可為以后興安落葉松林的采伐、經(jīng)營等提供基礎(chǔ)。

由于建模的數(shù)量較大，數(shù)據(jù)跨距的南北范圍較廣，溫度、降水、土壤類型等環(huán)境因子可能對樹木生長存在潛在的影響。因此，在未來的研究中有必要考慮氣候因子，以更精準(zhǔn)、更全面地描述樹木的生長。此外，構(gòu)建模型的競爭指數(shù)時(shí)只考慮了林木大小的相關(guān)參數(shù)沒有把樹種組成進(jìn)行討論在以后研究中可以將樹種組成加進(jìn)來。另外，混合效應(yīng)模型沒有明確的區(qū)域劃分，在以后研究中要考慮區(qū)域水平、樣地水平以及2水平的混合效應(yīng)模型，比較不同混合效應(yīng)的模型，找出最優(yōu)的混合效應(yīng)模型。