高東啟，鄧華鋒*，蔣 益，陳虹兵，程志楚

(1.北京林業(yè)大學(xué) 林學(xué)院，北京 100083;2.國(guó)家林業(yè)局調(diào)查規(guī)劃設(shè)計(jì)院，北京 100714)

單木生長(zhǎng)模型與林分生長(zhǎng)模型各有特點(diǎn)，利用單木生長(zhǎng)模型預(yù)測(cè)林分變量時(shí)，可以判定各單株木的生長(zhǎng)狀況和生長(zhǎng)潛力，但是存在著誤差積累、復(fù)雜性等缺點(diǎn);通過(guò)林分生長(zhǎng)模型預(yù)測(cè)林分變量時(shí)，可以直接提供林分的收獲量，但卻無(wú)法反映單木水平的詳細(xì)信息［1-5］，且單木生長(zhǎng)模型與林分生長(zhǎng)模型預(yù)測(cè)的同一林分變量可能不一致［5-6］。因此，有必要對(duì)單木生長(zhǎng)模型與林分生長(zhǎng)模型進(jìn)行整合，使兩類模型能夠有機(jī)的耦合，緊密地聯(lián)系起來(lái)形成一個(gè)統(tǒng)一的整體，從而保證預(yù)測(cè)結(jié)果的一致。鄧成等［6-7］提出了林分生長(zhǎng)與收獲模型的整體化研究思路，但是沒(méi)有考慮撫育間伐等經(jīng)營(yíng)管理措施對(duì)林分生長(zhǎng)的影響。

而現(xiàn)實(shí)的林分由于受撫育間伐等經(jīng)營(yíng)管理措施的影響，與自然生長(zhǎng)的林分相比，其生長(zhǎng)規(guī)律會(huì)發(fā)生變化［8］。在林分生長(zhǎng)與收獲預(yù)估模型的研究中，撫育間伐作為獨(dú)立于林齡、立地質(zhì)量、密度之外的第4 個(gè)建模因子越來(lái)越受研究者的重視［9］，相關(guān)學(xué)者已經(jīng)研究了撫育間伐對(duì)落葉松、杉木、側(cè)柏、馬尾松等的影響［10-16］，而間伐措施對(duì)北京市油松(Pinus tabulaeformis)生長(zhǎng)變化規(guī)律的影響還有待進(jìn)一步研究。本文以油松林分?jǐn)嗝娣e、蓄積量生長(zhǎng)模型為例，在考慮林分間伐效應(yīng)的基礎(chǔ)上，研究單木生長(zhǎng)模型與林分生長(zhǎng)模型耦合的方法。

1 材料與方法

1.1 研究地概況

北京市位于北緯39°28'～41°05'，東經(jīng)115°25'～117°30'，地處華北平原北端，北以燕山山地與內(nèi)蒙古高原接壤，西以太行山與陜西高原毗連，東北與松遼平原相通，南與黃淮海平原連片;屬于暖溫帶半濕潤(rùn)大陸性季風(fēng)型特點(diǎn)，四季分明，夏季炎熱多雨，冬季寒冷干燥，夏季降水量占全年降水量的74%。

1.2 數(shù)據(jù)來(lái)源與處理

本研究所采用的數(shù)據(jù)來(lái)源于國(guó)家林業(yè)局調(diào)查規(guī)劃設(shè)計(jì)院，為國(guó)家一類清查數(shù)據(jù)，分別調(diào)查于1996、2001、2006 年，每個(gè)樣地的面積為0.066 7 hm2。樣地主要調(diào)查因子有:林木胸徑、林分年齡、林分平均高、林分蓄積、枯損蓄積、采伐蓄積、郁閉度、水平距、坡向、坡位、坡度、海拔高度、土層厚度等，樣地固定且每隔5 年復(fù)測(cè)1 次，調(diào)查期內(nèi)部分林分進(jìn)行了撫育間伐。

從收據(jù)的數(shù)據(jù)中共挑選出北京市油松人工林樣地133 塊，經(jīng)過(guò)間伐的樣地69 塊，未經(jīng)間伐的樣地64 塊。對(duì)133 塊樣地進(jìn)行隨機(jī)抽樣，共有80 塊樣地用于建模，剩余53 塊樣地用于檢驗(yàn)，統(tǒng)計(jì)樣地基本情況如表1 所示，主要有林分的平均年齡、平均胸徑(即平方平均胸徑)、平均樹(shù)高、優(yōu)勢(shì)木平均高、平均每公頃株數(shù)、每公頃斷面積和每公頃蓄積。

油松樣地資料中沒(méi)有林分優(yōu)勢(shì)高，表1 中的林分優(yōu)勢(shì)高是通過(guò)林分優(yōu)勢(shì)高與林分平均高的關(guān)系計(jì)算出來(lái)的，參考賀姍姍［17］的研究結(jié)果，其通過(guò)解析優(yōu)勢(shì)木的方法，得到油松人工林優(yōu)勢(shì)高與林分平均高的相關(guān)關(guān)系為:H=1.115 h+0.429 3(H 為油松林分優(yōu)勢(shì)高，h 為油松林分平均高，決定系數(shù)為0.952 1)。通過(guò)上式計(jì)算出油松各樣地的林分優(yōu)勢(shì)高，查詢北京油松人工林的地位指數(shù)表［18］，作為評(píng)價(jià)林分立地質(zhì)量的指標(biāo)。

表1 樣地基本情況Tab.1 The condition of plots

表2 統(tǒng)計(jì)了隨機(jī)抽樣的樣地、樣木分布情況，其中，在建模數(shù)據(jù)中，間伐林分有41 塊樣地、2 859株樣木，未間伐林分有39 塊樣地、2 357株樣木;在檢驗(yàn)數(shù)據(jù)中，間伐林分有28 塊樣地、1 683株樣木，未間伐林分有25 塊樣地、1 427株樣木。

表2 樣地、樣木分布情況Tab.2 The distribution of plots and trees

1.3 研究方法

相關(guān)研究表明:通過(guò)單木生長(zhǎng)模型推導(dǎo)林分生長(zhǎng)模型的方法可以使得兩類模型有機(jī)的耦合，該方法由鄧成等［6-7］在林分生長(zhǎng)與收獲模型的整體化研究中提出，其核心是不同類型的模型之間的相關(guān)推導(dǎo)和變換，先尋找各個(gè)水平模型相互之間的關(guān)系，然后再進(jìn)行變形和轉(zhuǎn)換。但是鄧成等沒(méi)有考慮林分經(jīng)營(yíng)措施對(duì)林分生長(zhǎng)的影響，下面采用單木生長(zhǎng)模型推導(dǎo)林分生長(zhǎng)模型的方法對(duì)兩類模型進(jìn)行耦合，并考慮間伐措施對(duì)油松林分?jǐn)嗝娣e、蓄積量生長(zhǎng)的影響，首先建立單木水平生長(zhǎng)模型。

1.3.1 單木水平生長(zhǎng)模型的研究 單木水平主要是建立油松單木斷面積和材積生長(zhǎng)模型，單木生長(zhǎng)模型是地位指數(shù)(SI)、單木年齡(t)、林分密度和單木競(jìng)爭(zhēng)指標(biāo)等的函數(shù)［6-7，19］。以每公頃株數(shù)(N)作為反映林分密度的指標(biāo)，單木競(jìng)爭(zhēng)指標(biāo)選用單木的相對(duì)直徑Rd，選用Richards 方程研究單木斷面積和材積生長(zhǎng)模型，其基本形式如下:

式中，g 為單木斷面積，v 為單木材積，SI 為地位指數(shù)，N 為每公頃株數(shù)，Rd為相對(duì)直徑，t 為單木年齡，a1～a6為待定參數(shù)。

考慮到間伐林分與未間伐林分生長(zhǎng)規(guī)律不同，因此通過(guò)在參數(shù)中引入啞變量的方法［20-21］，將油松間伐林分、未間伐林分分別用啞變量(即定性代碼)來(lái)表示，從而可以將兩種類型的林分合在一起建立統(tǒng)一的油松單木生長(zhǎng)模型。經(jīng)過(guò)初步試驗(yàn)，發(fā)現(xiàn)啞變量主要影響模型(1)、(2)式中的參數(shù)a1，即將啞變量引入?yún)?shù)a1中時(shí)，模型的模擬效果最好，因此確定含啞變量的單木生長(zhǎng)模型的形式如下:

其中，g 為單木斷面積，v 為單木材積，K1、K2(引入的啞變量)分別為間伐林分和未間伐林分的定性代碼，SI 為地位指數(shù)，N 為每公頃株數(shù)，Rd為相對(duì)直徑，t 為單木年齡，a0～a6為待定參數(shù)。

1.3.2 林分水平生長(zhǎng)模型的研究 為了使得單木水平模型與林分水平模型能夠有機(jī)地耦合，下面采用上面的單木斷面積、材積生長(zhǎng)模型(3)、(4)式分別推導(dǎo)林分?jǐn)嗝娣e、蓄積量生長(zhǎng)模型。

對(duì)于單木斷面積、材積生長(zhǎng)模型，當(dāng)令(3)、(4)式中的相對(duì)直徑Rd=1 時(shí)，即單木的胸徑等于林分的平均胸徑，此時(shí)就變成了林分平均標(biāo)準(zhǔn)木的斷面積、材積生長(zhǎng)模型。在模型(3)、(4)式中，由于Rd=1，因此參數(shù)a5可去除，保持其他參數(shù)的形式不變，得到含啞變量的林分平均標(biāo)準(zhǔn)木的斷面積、材積生長(zhǎng)模型為:

對(duì)于林分平均標(biāo)準(zhǔn)木的斷面積、材積生長(zhǎng)模型(5)、(6)式，再乘以林分株數(shù)n，就得到了含啞變量的林分?jǐn)嗝娣e和蓄積量生長(zhǎng)預(yù)估模型，具體過(guò)程如下:

對(duì)于(7)、(8)式，為了使簡(jiǎn)化模型形式，減少模型中自變量的個(gè)數(shù)，進(jìn)行如下的變換。林分每公頃株數(shù)N=n/0.066 7，則有林分株數(shù)n=0.066 7N，將其代入(7)、(8)式，得到最終的林分?jǐn)嗝娣e、蓄積量生長(zhǎng)預(yù)估模型，形式為:

其中，G 為林分的斷面積，V 為林分的蓄積量，K1、K2(引入的啞變量)分別為間伐林分和未間伐林分的定性代碼，SI 為地位指數(shù)，N 為每公頃株數(shù)，t 為林分年齡，a0～a6為待定參數(shù)。

1.3.3 參數(shù)估計(jì)與模型檢驗(yàn) 綜合應(yīng)用ForStat［22］、Excel、SPSS 進(jìn)行數(shù)據(jù)計(jì)算和參數(shù)估計(jì)，并對(duì)模型進(jìn)行t 檢驗(yàn)，同時(shí)計(jì)算平均絕對(duì)偏差(MAD)、均方根誤差(RMSE)、決定系數(shù)(R2)和預(yù)估精度(P%)等幾個(gè)指標(biāo)來(lái)檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力，它們的數(shù)學(xué)表達(dá)式如下:

式中:yi為實(shí)測(cè)值(林分平均直徑、林分平均高)為模型預(yù)估值為模型預(yù)估值的平均值，n 為樣本數(shù)，p 為模型參數(shù)個(gè)數(shù)，t0.05為置信水平a=0.05 時(shí)的t 分布值。

2 結(jié)果與分析

對(duì)單木水平生長(zhǎng)模型(3)、(4)式進(jìn)行參數(shù)估計(jì)，統(tǒng)計(jì)油松單木斷面積、材積生長(zhǎng)模型參數(shù)a1～a6的估計(jì)值和R2如下表3 所示，其中油松單木斷面積生長(zhǎng)模型的決定系數(shù)為0.876 8，單木材積生長(zhǎng)模型的決定系數(shù)為0.864 5，R2都較高，說(shuō)明模擬效果較好。

表3 單木水平生長(zhǎng)模型參數(shù)統(tǒng)計(jì)Tab.3 Parameter statistics of individual tree-level growth models

利用建模數(shù)據(jù)和檢驗(yàn)數(shù)據(jù)對(duì)單木水平生長(zhǎng)模型(3)、(4)式進(jìn)行檢驗(yàn)，統(tǒng)計(jì)平均絕對(duì)偏差(MAD)、均方根誤差(RMSE)、決定系數(shù)(R2)、預(yù)估精度(P%)等指標(biāo)如下表4 所示。單木斷面積生長(zhǎng)模型對(duì)建模數(shù)據(jù)的預(yù)測(cè)精度和R2分別為99.25%、0.876 8，對(duì)于檢驗(yàn)數(shù)據(jù)的預(yù)測(cè)精度和R2分別為98.64%、0.787 1，MAD 和RMSE 等誤差都較小，不超過(guò)0.005;單木材積生長(zhǎng)模型對(duì)建模數(shù)據(jù)的預(yù)測(cè)精度和R2分別為99.03%、0.864 5，對(duì)于檢驗(yàn)數(shù)據(jù)的預(yù)測(cè)精度和R2分別為98.41%、0.777 2，MAD 和RMSE 均小于0.021，誤差較小，說(shuō)明模型預(yù)測(cè)效果很好。進(jìn)行t 檢驗(yàn)，結(jié)果為Sig.＜0.05，說(shuō)明實(shí)測(cè)值與預(yù)測(cè)值在0.05 的顯著水平上無(wú)顯著差異，模型建立的比較合理。

表4 單木水平生長(zhǎng)模型檢驗(yàn)Tab.4 Inspection of individual tree-level growth models

進(jìn)一步檢驗(yàn)單木生長(zhǎng)模型(3)、(4)式對(duì)間伐林分和未間伐林分的預(yù)測(cè)效果，統(tǒng)計(jì)相關(guān)結(jié)果如下表5 所示。單木模型對(duì)間伐林分單木斷面積的預(yù)測(cè)精度和R2分別為98.77%、0.885 8，MAD 為0.002 1，RMSE 為0.003 1;對(duì)間伐林分單木材積的預(yù)測(cè)精度和R2分別為98.55%、0.881 1，MAD 為0.009 9，RMSE為0.014 7。單木模型對(duì)未間伐林分單木斷面積的預(yù)測(cè)精度和R2分別為97.22%、0.700 3，MAD 為0.003 6，RMSE 為0.005 6;對(duì)材積的預(yù)測(cè)精度和R2分別為96.74%、0.680 6，MAD 為0.016 9，RMSE 為0.026 3。單木模型對(duì)間伐林分和未間伐林分的預(yù)測(cè)精度都在96%以上，各項(xiàng)誤差均較小，說(shuō)明所建啞變量模型比較合理。t 檢驗(yàn)結(jié)果都為Sig.＜0.05，即實(shí)測(cè)值與預(yù)測(cè)值在0.05的顯著水平上無(wú)顯著差異。

表5 間伐林分和未間伐林分單木生長(zhǎng)模型預(yù)測(cè)結(jié)果檢驗(yàn)Tab.5 Inspectiing forecast results of tree growth models for thinned stands and unthinned stands

由于林分水平生長(zhǎng)模型(9)、(10)式是通過(guò)單木水平生長(zhǎng)模型(3)、(4)式推導(dǎo)而來(lái)的，因此將相關(guān)的參數(shù)a1～a6(表3)帶入(9)、(10)式得到林分水平生長(zhǎng)模型(11)、(12)式:

利用建模數(shù)據(jù)和檢驗(yàn)數(shù)據(jù)對(duì)林分水平生長(zhǎng)模型(15)式、(16)式進(jìn)行檢驗(yàn)，結(jié)果如下表6 所示。對(duì)于建模數(shù)據(jù)而言，林分?jǐn)嗝娣e生長(zhǎng)模型的預(yù)測(cè)精度和R2分別為94.17%、0.848 5，MAD 為0.124 8，RMSE為0.167 9;林分蓄積量生長(zhǎng)模型的預(yù)測(cè)精度和R2分別為92.60%、0.806 5，MAD 為0.586 8，RMSE 為0.810 0。對(duì)于檢驗(yàn)數(shù)據(jù)而言，林分?jǐn)嗝娣e生長(zhǎng)模型的預(yù)測(cè)精度和R2分別為88.59%、0.763 4，MAD 為0.165 0，RMSE 為0.273 2;林分蓄積量生長(zhǎng)模型的預(yù)測(cè)精度和R2分別為86.44%、0.711 9，MAD 為0.775 8，RMSE 為1.280 1。通過(guò)檢驗(yàn)，林分?jǐn)嗝娣e、蓄積量生長(zhǎng)模型的預(yù)估精度較高超過(guò)了86%，誤差較小，t 檢驗(yàn)結(jié)果為Sig.＜0.05，說(shuō)明實(shí)測(cè)值與預(yù)測(cè)值無(wú)顯著差異，模型比較合理，含啞變量的林分模型使得油松間伐林分與未間伐林分能夠相容從而合到一起。

表6 林分水平生長(zhǎng)模型檢驗(yàn)Tab.6 Inspection of stand-level growth models

進(jìn)一步檢驗(yàn)林分生長(zhǎng)模型(15)、(16)式對(duì)間伐林分和未間伐林分的預(yù)測(cè)效果，結(jié)果如下表7 所示。林分模型對(duì)間伐林分?jǐn)嗝娣e的預(yù)測(cè)精度和R2分別為90.68%、0.832 0，對(duì)間伐林分蓄積量的預(yù)測(cè)精度和R2分別為89.08%、0.799 2;林分模型對(duì)未間伐林分?jǐn)嗝娣e的預(yù)測(cè)精度和R2分別為87.01%、0.793 0，對(duì)未間伐林分蓄積量的預(yù)測(cè)精度和R2分別為85.90%、0.739 8，模型對(duì)間伐林分的預(yù)測(cè)效果要稍微好于未間伐林分。進(jìn)行t 檢驗(yàn)，結(jié)果為Sig.＜0.05，說(shuō)明實(shí)測(cè)值與預(yù)測(cè)值在0.05 的顯著水平上無(wú)顯著差異;同時(shí)林分模型對(duì)間伐林分和未間伐林分的預(yù)測(cè)精度都在85%以上，說(shuō)明所建啞變量模型比較合理，預(yù)測(cè)效果較好。

表7 間伐林分和未間伐林分生長(zhǎng)模型預(yù)測(cè)結(jié)果檢驗(yàn)Tab.7 Inspectiing forecast results of stands growth models for thinned stands and unthinned stands

3 結(jié)論與討論

本文在考慮林分間伐效應(yīng)的基礎(chǔ)上，引入啞變量將間伐林分與未間伐林分整合到一起建立統(tǒng)一的單木斷面積、材積生長(zhǎng)模型，進(jìn)一步對(duì)模型進(jìn)行變換推導(dǎo)出林分生長(zhǎng)模型，從而使得單木、林分兩個(gè)水平的模型能夠有機(jī)的耦合，形成一個(gè)聯(lián)系緊密的整體，使模型更具有生物學(xué)意義和可解釋性。

所建單木模型對(duì)油松單木斷面積、材積的預(yù)測(cè)精度在98%以上，對(duì)間伐林分、未間伐林分的預(yù)測(cè)精度都在96%以上;而推導(dǎo)出的林分模型對(duì)油松林分?jǐn)嗝娣e、蓄積量的預(yù)測(cè)精度分別為88.59%、86.44%，對(duì)間伐林分、未間伐林分的預(yù)測(cè)精度都在85%以上。模型對(duì)林分蓄積量的預(yù)測(cè)精度稍低于林分?jǐn)嗝娣e，可能的原因是林分蓄積量由林分中各單木的材積推算而來(lái)，而單木材積并非實(shí)測(cè)值，而是由北京市一元立木材積表所得，存在一定的誤差，所以預(yù)測(cè)精度稍低。

單木生長(zhǎng)模型推導(dǎo)林分生長(zhǎng)模型的方法，使得兩個(gè)水平的模型能夠整合在一起，并保證了兩個(gè)水平模型對(duì)林分?jǐn)嗝娣e、蓄積量預(yù)測(cè)精度的一致，該方法不僅適用于預(yù)測(cè)林分?jǐn)嗝娣e、蓄積量，也適用于預(yù)測(cè)其他林分因子。本文所用數(shù)據(jù)調(diào)查時(shí)間不同，沒(méi)有考慮時(shí)間效應(yīng)對(duì)建模的影響，下步可建立混合模型進(jìn)行更深入的研究。

本文所建模型對(duì)北京地區(qū)的油松人工林是適用的，可以用來(lái)模擬其單木或林分的生長(zhǎng)，從而為油松林分的經(jīng)營(yíng)管理決策提供參考依據(jù)。

［1］趙俊卉.長(zhǎng)白山云冷杉混交林生長(zhǎng)模型的研究［D］.北京:北京林業(yè)大學(xué)，2010.

［2］孟憲宇.測(cè)樹(shù)學(xué)［M］.3 版.北京:中國(guó)林業(yè)出版社，2008.

［3］Qin J，Cao Q V.Using disaggregation to link individual-tree and whole-stand growth models［J］.Canadian Journal of Forest Research，2006，36(4):953-960.

［4］張雄清，雷淵才，陳新美，等.組合預(yù)測(cè)法在林分?jǐn)嗝娣e生長(zhǎng)預(yù)估中的應(yīng)用［J］.北京林業(yè)大學(xué)學(xué)報(bào)，2010，32(4):6-11.

［5］胡煥香，佘濟(jì)云，李俊，等.湖北桂花林場(chǎng)檫木次生林單木生長(zhǎng)模型的研究［J］.中南林業(yè)科技大學(xué)學(xué)報(bào)，2013，33(4):61-65.

［6］鄧成.林分生長(zhǎng)和收獲模型整體化研究［D］.長(zhǎng)沙:中南林業(yè)科技大學(xué)，2008.

［7］鄧成，呂勇，雷淵才，等.林分生長(zhǎng)和收獲模型整體化思路探討［J］.林業(yè)調(diào)查規(guī)劃，2008，33(4):7-9.

［8］謝華.杉木人工林撫育間伐效應(yīng)研究［D］.合肥:安徽農(nóng)業(yè)大學(xué)，2007.

［9］Sun H G，Zhang J G，Duan A G，et al.A review of stand basal area growth models［J］.Forestry Studies in China，2007，9(1):85-94.

［10］畢永華.興安落葉松林不同間伐強(qiáng)度林分生長(zhǎng)模型的研究［D］.哈爾濱:東北林業(yè)大學(xué)，2001.

［11］雷相東，陸元昌，張會(huì)儒，等.撫育間伐對(duì)落葉松云冷杉混交林的影響［J］.林業(yè)科學(xué)，2005，41(4):78-85.

［12］李春明，杜紀(jì)山，張會(huì)儒.間伐林分的斷面積生長(zhǎng)模型研究［J］.林業(yè)資源管理，2004(3):52-55.

［13］李春明，杜紀(jì)山，張會(huì)儒.撫育間伐對(duì)人工落葉松斷面積和蓄積生長(zhǎng)的影響［J］.林業(yè)資源管理，2007(3):90-93.

［14］段劼，馬履一，賈黎明，等.撫育間伐對(duì)側(cè)柏人工林及林下植被生長(zhǎng)的影響［J］.生態(tài)學(xué)報(bào)，2010，30(6):1431-1441.

［15］明安剛，張治軍，諶紅輝，等.撫育間伐對(duì)馬尾松人工林生物量與碳貯量的影響［J］.林業(yè)科學(xué)，2013，49(10):1-6.

［16］趙西哲.湘中南地區(qū)主要用材林森林撫育綜合效益研究［D］.長(zhǎng)沙:中南林業(yè)科技大學(xué)，2007.

［17］賀珊珊.北京山區(qū)油松人工林林分結(jié)構(gòu)與生長(zhǎng)結(jié)構(gòu)研究［D］.北京:北京林業(yè)大學(xué)，2009.

［18］段劼，馬履一，賈黎明，等.北京低山地區(qū)油松人工林立地指數(shù)表的編制及應(yīng)用［J］.林業(yè)科學(xué)，2009，45(3):7-12.

［19］郎榮，許建初，Timm Tennigkeit，等.基于樣方數(shù)據(jù)的云南松林分生長(zhǎng)模型研究—以云南省保山市楊柳白族彝族鄉(xiāng)為例［J］.植物分類與資源學(xué)報(bào)，2011，33(3):357-363.

［20］李忠國(guó)，孫曉梅，陳東升，等.基于啞變量的日本落葉松生長(zhǎng)模型研究［J］.西北農(nóng)林科技大學(xué)學(xué)報(bào)，2011，39(8):69-74.

［21］高東啟，鄧華鋒，王海賓，等.基于啞變量的蒙古櫟林分生長(zhǎng)模型研究［J］.東北林業(yè)大學(xué)學(xué)報(bào)，2014，42(1):61-64.

［22］唐守正，郎奎建，李?？?統(tǒng)計(jì)與生物數(shù)學(xué)模型計(jì)算(ForStat 教程)［M］.北京:科學(xué)出版社，2009.