劉恩斌 姚鴻文 任澤茜 周國(guó)模 杜華強(qiáng)

(1.浙江農(nóng)林大學(xué)省部共建亞熱帶森林培育國(guó)家重點(diǎn)實(shí)驗(yàn)室 浙江省森林生態(tài)系統(tǒng)碳循環(huán)與固碳減排重點(diǎn)實(shí)驗(yàn)室 浙江農(nóng)林大學(xué)環(huán)境與資源學(xué)院 杭州 311300;2.倫敦大學(xué)學(xué)院地理系 倫敦 WC1E 6BT；3.浙江省森林資源監(jiān)測(cè)中心 杭州 310020)

毛竹(Phyllostachysedulis)是我國(guó)南方重要的森林資源之一，生長(zhǎng)快、成材早、產(chǎn)量高、伐期短、收益期長(zhǎng),不僅具有較高的經(jīng)濟(jì)價(jià)值和社會(huì)價(jià)值，而且其光合能力較強(qiáng)，全年都能維持良好的生物量和碳儲(chǔ)量，具有重要的生態(tài)功能(周國(guó)模等，2004；方晰等，2002；李意德等，1998；阮宏華等，1997)。毛竹生長(zhǎng)過程中，主要依靠地下莖(俗稱竹鞭)上的筍芽逐步發(fā)育成竹筍，再長(zhǎng)成新竹，新竹在第一年內(nèi)即可完成胸徑、高度生長(zhǎng)，之后雖然胸徑和竹高不再發(fā)生改變，但生物量會(huì)隨年齡增長(zhǎng)而變化，因此，胸徑和年齡是影響毛竹林生物量易測(cè)易辨的獨(dú)立變量。

毛竹林胸徑與年齡聯(lián)合分布較一元分布更能完整反映其資源結(jié)構(gòu)特征，精確預(yù)估各徑階、各齡階株數(shù)以及產(chǎn)量、枯損量等，為毛竹林經(jīng)營(yíng)決策提供科學(xué)依據(jù)；且結(jié)合毛竹胸徑年齡二元單竹生物量模型(周國(guó)模，2006)，還能對(duì)區(qū)域尺度毛竹林生物量進(jìn)行精確估算，為毛竹林碳庫(kù)精準(zhǔn)測(cè)算奠定基礎(chǔ)。到目前為止，二元廣義β、二元Johnson’s Sbb、二元Weibull等分布函數(shù)常被用來(lái)描述森林結(jié)構(gòu)特征因子的二元分布，如采用二元廣義β和 Johnson’s Sbb分布函數(shù)描述胸徑樹高二元聯(lián)合分布(Lietal.，2002；Gorgosovarelaetal.，2019；Cardil Forradellasetal.，2016；Wangetal.，2007；Mnness，2015)基于二元Weibull和二元β分布函數(shù)描述樹冠內(nèi)細(xì)枝葉面積分布的二維結(jié)構(gòu)以及毛竹胸徑年齡聯(lián)合分布(Temesgenetal.，2003；葛宏立等，2008；劉恩斌等，2010)等，森林結(jié)構(gòu)特征因子的概率分布對(duì)森林生物量估算及森林經(jīng)營(yíng)具有非常重要的作用。但已有森林結(jié)構(gòu)特征因子二元聯(lián)合分布研究方法存在如下缺陷：1)2個(gè)森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)類型必須相同，如二元β、二元Weibull和二元Sbb的邊緣分布函數(shù)必須是一元β、一元Weibull和一元Sbb分布函數(shù)，在實(shí)際應(yīng)用過程中，當(dāng)2個(gè)森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)類型不同時(shí)，常用二元分布函數(shù)不能描述2個(gè)因子的聯(lián)合分布，適用性不廣；2)不能由邊緣分布確定對(duì)應(yīng)的聯(lián)合分布，森林結(jié)構(gòu)特征因子二元聯(lián)合分布研究通常先建立聯(lián)合分布函數(shù)，后推導(dǎo)聯(lián)合分布對(duì)應(yīng)的邊緣分布(葛宏立等，2008；劉恩斌等，2010；Hafleyetal.，1976)，建立的聯(lián)合分布函數(shù)不是由對(duì)應(yīng)的邊緣分布函數(shù)確定的，事實(shí)上，邊緣分布函數(shù)的類型和理論較聯(lián)合分布函數(shù)更多、更完善，但理論和應(yīng)用價(jià)值遠(yuǎn)不及聯(lián)合分布函數(shù)，根據(jù)森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)建立其對(duì)應(yīng)的二元聯(lián)合分布函數(shù)具有重要意義；3)擬合參數(shù)需要森林結(jié)構(gòu)特征因子聯(lián)合分布(密度)的實(shí)測(cè)值，這在很多情況下難以滿足，如要研究2014年浙江省毛竹胸徑年齡的二元分布結(jié)構(gòu)，從已有數(shù)據(jù)資料能夠獲得每塊樣地各徑階、各齡級(jí)的毛竹分別有多少株，但無(wú)法獲得每塊樣地中每株毛竹的胸徑和年齡實(shí)測(cè)數(shù)據(jù)，即只能得到全省毛竹胸徑和年齡的邊緣分布值，無(wú)法獲得全省毛竹胸徑和年齡的聯(lián)合分布值；再如要研究某區(qū)域杉木(Cunninghamialanceolata)胸徑樹高二元聯(lián)合分布，從各種文獻(xiàn)中獲取滿足該區(qū)域抽樣精度的2套樣本，一套可得到胸徑分布值，另一套可得到樹高分布值，但由這2套樣本無(wú)法得到杉木胸徑樹高聯(lián)合分布值，也就無(wú)法使用常用二元分布函數(shù)構(gòu)建其對(duì)應(yīng)的聯(lián)合分布模型。

隨著科學(xué)研究的逐步推進(jìn)，在統(tǒng)計(jì)領(lǐng)域快速發(fā)展的Copula函數(shù)理論為解決上述相關(guān)非正態(tài)變量建模問題提供了新思路(Nelsen，1998)。Copula函數(shù)理論的核心思想是任意一個(gè)二元聯(lián)合分布函數(shù)均可拆分為1個(gè)包含相關(guān)性大小和相關(guān)結(jié)構(gòu)信息的Copula分布(密度)函數(shù)與2個(gè)對(duì)應(yīng)的邊緣分布函數(shù)，其在建模時(shí)的靈活性和適應(yīng)性具有無(wú)可比擬的優(yōu)勢(shì)，如2個(gè)森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)類型可以不同、可由邊緣分布函數(shù)確定其對(duì)應(yīng)的聯(lián)合分布(密度)函數(shù)、對(duì)實(shí)測(cè)值依賴性不強(qiáng)等，已在金融和水文等領(lǐng)域廣泛應(yīng)用(張堯庭，2002；謝華等，2012；Sklar，1959；Fangetal.，2020；Wangetal.，2016)。在林業(yè)領(lǐng)域也有相關(guān)研究，如采用三元Copula函數(shù)對(duì)單株材積進(jìn)行估算(Petrauskasetal.，2011；Wangetal.，2010)，采用二元Copula函數(shù)描述樹木胸徑樹高聯(lián)合分布與樹高曲線(Wangetal.，2008；Oganaetal.，2020；Gorgoso-Varelaetal.，2016；Rup?ysetal.，2012)，還用于河岸林下植被蓋度的估算(Eskelsonetal.，2011)以及林分結(jié)構(gòu)特征因子的模擬(Kershawetal.，2010)等，但用于描述毛竹林分結(jié)構(gòu)特征因子聯(lián)合分布的研究鮮見報(bào)道。鑒于此，本研究針對(duì)森林結(jié)構(gòu)特征因子常用二元聯(lián)合分布研究方法存在的缺陷，選擇適用條件低、適應(yīng)范圍廣、應(yīng)用價(jià)值大的二元Copula密度(分布)函數(shù)構(gòu)建毛竹胸徑年齡二元聯(lián)合分布模型，以期為精確測(cè)量森林結(jié)構(gòu)特征因子聯(lián)合分布提供參考。

1 材料與方法

1.1 試驗(yàn)材料

浙江省于1979年建立森林資源連續(xù)清查體系，以5年為一個(gè)復(fù)查周期，固定樣地4 250塊，樣點(diǎn)格網(wǎng)大小為4 km×6 km，樣地形狀為正方形，邊長(zhǎng)28.28 m，面積800 m2。本研究利用2009年浙江省177塊毛竹連續(xù)清查固定樣地?cái)?shù)據(jù)進(jìn)行描述。樣地調(diào)查因子包括土層厚度、海拔、坡位、坡度、坡向、樣地毛竹平均胸徑、樣地毛竹株數(shù)以及樣地內(nèi)每株毛竹的胸徑和年齡等。經(jīng)統(tǒng)計(jì)，各樣地毛竹株數(shù)22～897不等，胸徑5～15 cm，年齡1～5度(當(dāng)年生毛竹記為1度竹，2～3年生毛竹記為2度竹，依此類推)，具體數(shù)據(jù)見表1。

表1 浙江省毛竹林胸徑、年齡統(tǒng)計(jì)數(shù)據(jù)Tab.1 Statistical data of DBH and age of moso bamboo forest in Zhejiang Province

1.2 研究方法

1.2.1 基本公式 借用Coupla函數(shù)理論，設(shè)D為胸徑、A為年齡，令u=F(D)和v=G(A)分別為毛竹胸徑和年齡的邊緣分布函數(shù)，則對(duì)于毛竹胸徑年齡聯(lián)合分布函數(shù)H(D,A)，存在著一個(gè)Copula分布函數(shù)C，使得

H(D,A)=C[F(D),G(A)]。

(1)

設(shè)毛竹胸徑年齡的二元Copula分布函數(shù)和密度函數(shù)分別為C(u,v)和D(u,v)，則由式(1)得：

(2)

令f(D)=?u/?D為胸徑D的邊緣概率密度函數(shù)，g(A)=?v/?A為年齡A的邊緣概率密度函數(shù)，則毛竹胸徑年齡的聯(lián)合概率密度函數(shù)h(D,A)為：

(3)

有關(guān)Copula密度函數(shù)更詳細(xì)的介紹參考有關(guān)文獻(xiàn)(Sklar，1959；Fangetal.，2020；Wangetal.，2016)。具體應(yīng)用時(shí)，可根據(jù)毛竹連續(xù)清查數(shù)據(jù)確定其胸徑和年齡的邊緣分布函數(shù)值。

常用Copula密度函數(shù)包括正態(tài)Copula、t-Copula和阿基米德Copula密度函數(shù)，其中Gumbel Copula、Clayton Copula和Frank Copula密度函數(shù)是最常用的3種阿基米德Copula密度函數(shù)(張蕾，2017)。

1.2.2 最優(yōu)Copula密度函數(shù)識(shí)別方法 二元正態(tài)Copula和二元Frank Copula密度函數(shù)適合描述具有對(duì)稱尾部且尾部漸進(jìn)獨(dú)立的二維森林結(jié)構(gòu)特征因子，二元t-Copula密度函數(shù)適合描述具有對(duì)稱尾部且尾部相關(guān)的二維森林結(jié)構(gòu)特征因子，二元Gumbel Copula密度函數(shù)適合描述具有非對(duì)稱尾部且上尾相關(guān)、下尾漸進(jìn)獨(dú)立的二維森林結(jié)構(gòu)特征因子，二元Clayton Copula密度函數(shù)適合描述具有非對(duì)稱尾部且下尾相關(guān)、上尾漸進(jìn)獨(dú)立的二維森林結(jié)構(gòu)特征因子。

赤池信息準(zhǔn)則(Akaike information criterion，AIC)是衡量統(tǒng)計(jì)模型擬合優(yōu)良性(goodness of fit)的一種標(biāo)準(zhǔn)，其以極大似然函數(shù)為基礎(chǔ)評(píng)價(jià)概率密度模型擬合數(shù)據(jù)的優(yōu)良性。本研究采用AIC識(shí)別最優(yōu)Copula密度函數(shù)，計(jì)算公式為：

(4)

式中：m為所選模型參數(shù)個(gè)數(shù)；L為所選分布的極大似然函數(shù)；N為擬合的原始數(shù)據(jù)樣本量；當(dāng)所選模型為二元 Gaussian Copula、二元Frank Copula、二元Gumbel Copula和二元tClayton Copula密度函數(shù)時(shí)，m=1；當(dāng)所選模型為二元-Copula密度函數(shù)時(shí)，m=2；ui,vi(i=1,2,…,N)為原始樣本數(shù)據(jù)第i個(gè)胸徑和年齡的邊緣分布函數(shù)值；θ為所選模型參數(shù)；D(ui,vi;θ)為所選模型的概率密度函數(shù)值。

對(duì)比所有備選密度函數(shù)的AIC，AIC越小，表明密度函數(shù)擬合原始數(shù)據(jù)的性能越好。

1.2.3 毛竹胸徑年齡二元Copula密度模型 為更好說明二元Copula密度函數(shù)的性能，利用2009年浙江省177塊毛竹連續(xù)清查固定樣地?cái)?shù)據(jù)，采用Matlab軟件編寫極大似然算法，應(yīng)用上述5個(gè)二元Copula密度函數(shù)對(duì)毛竹胸徑年齡的聯(lián)合密度函數(shù)進(jìn)行研究，并根據(jù)AIC從5個(gè)二元Copula密度函數(shù)中優(yōu)選出1個(gè)最優(yōu)二元Copula函數(shù)，進(jìn)而構(gòu)建浙江省毛竹胸徑年齡二元Copula密度模型。

1.2.4 最優(yōu)Copula密度函數(shù)與常用二元分布函數(shù)的對(duì)比分析 從測(cè)量精度、參數(shù)個(gè)數(shù)和參數(shù)迭代收斂性3方面，對(duì)毛竹胸徑年齡最優(yōu)Copula密度函數(shù)、二元Weibull分布函數(shù)和二元Johnson’s Sbb密度函數(shù)進(jìn)行深入分析。由于二元廣義β概率密度函數(shù)參數(shù)較多，采用已有參數(shù)估算方法(Lietal.,2002)多次嘗試后迭代始終不收斂,故本研究未選用。

采用Matlab軟件編寫程序并進(jìn)行繪圖。

2 結(jié)果與分析

2.1 最優(yōu)二元Copula密度函數(shù)的選取

2.1.1 毛竹胸徑年齡累計(jì)分布二維頻數(shù)直方圖 根據(jù)表1得出毛竹胸徑和年齡的累計(jì)分布值，繪制毛竹胸徑年齡累計(jì)分布的二維頻數(shù)直方圖(圖1)。令u、v分別為胸徑(D)和年齡(A)的累計(jì)分布函數(shù)值，則圖1即(u,v)聯(lián)合頻數(shù)直方圖，可作為(u,v)二元Copula密度函數(shù)圖像的形狀?？梢钥闯?，毛竹胸徑年齡密度函數(shù)具有不對(duì)稱尾部，上尾相關(guān)性大、下尾漸進(jìn)獨(dú)立。

圖1 毛竹胸徑年齡累計(jì)分布的二維頻數(shù)直方圖Fig.1 The bivariate frequency histogram of cumulative distribution of DBH and age of moso bamboo

2.1.2 毛竹胸徑年齡的二維Copula密度函數(shù)圖 根據(jù)表1得出毛竹胸徑和年齡的累計(jì)分布值，繪制毛竹胸徑年齡各二元Copula聯(lián)合密度圖如圖2所示。可以看出：1)二元Gumbel Copula和二元t-Copula聯(lián)合密度圖與圖1接近，其他二元Copula聯(lián)合密度圖與圖1存在較大差異，故優(yōu)先選擇二元Gumbel Copula和二元t-Copula密度函數(shù)測(cè)量毛竹胸徑年齡聯(lián)合分布；2)二元正態(tài)Copula、二元t-Copula和二元Frank Copula密度函數(shù)具有對(duì)稱尾部，無(wú)法描述毛竹胸徑年齡二元結(jié)構(gòu)的非對(duì)稱尾部相關(guān)性；二元Gumbel Copula密度函數(shù)具有非對(duì)稱尾部，曲線呈J形，說明該函數(shù)對(duì)毛竹胸徑年齡二元結(jié)構(gòu)的上尾部變化較為敏感，能更好描述下尾低、上尾高的尾部相關(guān)性；二元Clayton Copula密度函數(shù)也具有非對(duì)稱尾部，曲線呈L形，說明該函數(shù)對(duì)胸徑年齡二元結(jié)構(gòu)的下尾部變化較敏感，能更好地描述下尾高、上尾低的尾部相關(guān)性。

2.1.3 基于AIC的Copula函數(shù)優(yōu)選 根據(jù)式(4)，計(jì)算二元t-Copula和二元Gumbel Copula密度函數(shù)的AIC分別為 -14.310 4和-19.519 6，確定系數(shù)(R2)分別為0.967 6和0.984 1，可見二元Gumbel Copula密度函數(shù)描述毛竹胸徑年齡二元聯(lián)合分布效果最好。

2.2 二元Gumbel Copula分布(密度)函數(shù)與常用二元分布(密度)函數(shù)的對(duì)比分析

以胸徑和年齡的累計(jì)分布為自變量，采用極大似然法擬合二元Gumbel Copula密度函數(shù)參數(shù)α為1.055 0，R2=0.984 1，以胸徑和年齡為自變量，擬合二元Weibull分布函數(shù)與二元Sbb密度函數(shù)參數(shù)分別為(a1,b1,c1,a2,b2,c2,r)=(-0.163 7,8.917 0,4.211 9,-1.192 6,3.573 4,3.506 2,1.003 2)和(ρ,ξ1,ξ2,λ1,λ2,γ1,γ2,δ1,δ2)=(0.082 5,1.903 1,-2.530 1,15.146 4,7.530 1,0.610 2,1.582 0,1.386 2,0.072 5)，R2分別為0.990 1和0.736 2。毛竹胸徑年齡實(shí)測(cè)概率以及各概率密度函數(shù)的估計(jì)值如表2所示。

表2 實(shí)測(cè)概率以及二元Gumbel Copula(GumCp)密度函數(shù)、二元Weibull分布函數(shù)和二元Sbb密度函數(shù)估計(jì)值Tab.2 The measured probability and estimated values of the bivariate Gumbel Copula density function, the bivariate Weibull distribution function and the bivariate Sbb density function

胸徑DBH/cm5度 5du實(shí)測(cè)概率Measured probabilityGumCp估計(jì)值GumCp estimatedvaluesSbb估計(jì)值Sbb estimatedvaluesWeibull估計(jì)值Weibull estimatedvalues50.001 50.001 50.001 40.002 360.002 20.001 70.002 60.002 270.004 00.002 50.003 60.003 380.004 20.003 10.004 20.004 090.004 20.003 50.004 10.004 0100.004 60.003 90.003 60.003 5 110.002 90.003 20.002 70.002 4 120.001 60.003 80.001 70.001 2 130.000 50.005 20.000 80.000 4 140.000 30.007 00.000 30.000 21500.000 60.000 10

采用柯爾莫哥洛夫檢驗(yàn)各分布(密度)函數(shù)的擬合優(yōu)度(王巖等，2006)，此時(shí)二元Gumbel Copula密度函數(shù)、二元Weibull分布函數(shù)和二元Sbb密度函數(shù)累計(jì)值的最大偏差分別為D55,Copula=0.015 8、D55,Weibull=0.007 0、D55,Sbb=0.078 1，在0.05顯著性水平下的閾值為D55,0.05=0.179 8，可見毛竹胸徑年齡均服從這3種分布(密度)函數(shù)。

從表2可以看出：1)二元Weibull分布函數(shù)的測(cè)量精度比二元Gumbel Copula密度函數(shù)略高，二元Sbb密度函數(shù)測(cè)量精度最低；2)二元Weibull分布函數(shù)有7個(gè)參數(shù)，二元Sbb密度函數(shù)有9個(gè)參數(shù)，而二元Gumbel Copula密度函數(shù)只有1個(gè)參數(shù)，因參數(shù)越多越難收斂，故二元Gumbel Copula密度函數(shù)在估測(cè)精度非常高的情況下參數(shù)最易收斂；3)二元Weibull密度函數(shù)由分布函數(shù)求二階偏導(dǎo)才能得到，導(dǎo)出的密度函數(shù)非常復(fù)雜，用二元Weibull分布函數(shù)估算森林結(jié)構(gòu)特征因子二元概率密度時(shí)，首先要將原始數(shù)據(jù)變換成聯(lián)合分布值，然后再將估算出的二元分布值變換成相應(yīng)的聯(lián)合密度值，聯(lián)合密度值估算較繁瑣，而用二元Gumbel Copula密度函數(shù)可直接得到森林結(jié)構(gòu)特征因子二元概率密度的估計(jì)值，估測(cè)簡(jiǎn)單快捷；4)二元Gumbel Copula密度函數(shù)對(duì)小徑階、低齡級(jí)與大徑階、高齡級(jí)的測(cè)量值較實(shí)測(cè)值偏大，這也是該Copula密度函數(shù)本身固有的缺陷，混合Copula密度函數(shù)可克服這一缺陷，有關(guān)二元混合Copula函數(shù)的建模有待進(jìn)一步研究。

2.3 毛竹胸徑年齡二元概率密度圖與二元Gumbel Copula分布圖

根據(jù)表1得出毛竹胸徑和年齡的邊緣分布值，由式(1)計(jì)算毛竹胸徑年齡二元Copula分布函數(shù)值，繪制圖形如圖3所示?？梢钥闯觯弘S著胸徑和年齡累計(jì)分布值增大，毛竹胸徑年齡聯(lián)合分布值也增大，并在兩邊緣分布值取1和0時(shí)聯(lián)合分布值達(dá)到最大和最小，與實(shí)際相符。

圖3 毛竹胸徑年齡二元Gumbel Copula分布函數(shù)圖Fig.3 The joint distribution between DBH and age of moso bamboo using the bivariate Gumbel Copula distribution function

根據(jù)式(2)，可由毛竹胸徑年齡二元Gumbel Copula分布函數(shù)推導(dǎo)出相應(yīng)的二元密度函數(shù)，再根據(jù)式(3)，可計(jì)算毛竹胸徑年齡二元密度值，以胸徑和年齡邊緣分布值為橫軸和縱軸，胸徑年齡二元密度值為豎軸，繪制圖形如圖4所示?？梢钥闯觯?)由于立地條件、人為經(jīng)營(yíng)等因素影響，毛竹胸徑和年齡邊緣分布值具有一定分布規(guī)律；2)毛竹胸徑年齡二元密度值并不隨胸徑和年齡邊緣分布值增大而增大，當(dāng)胸徑累計(jì)概率為0.664 6和0.999 8、年齡累計(jì)概率為0.824 2和0.973 8時(shí)，胸徑年齡二元密度值達(dá)到最大和最小。

3 討論

3.1 二元Copula分布(密度)函數(shù)與常用二元分布(密度)函數(shù)的比較

1)二者自變量不同：常用二元概率分布(密度)函數(shù)，如二元Weibull、二元Sbb、二元β分布函數(shù)等，其自變量為森林結(jié)構(gòu)特征因子，因變量為對(duì)應(yīng)的聯(lián)合分布(密度)值，而二元Copula分布(密度)函數(shù)的自變量為森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)值，因變量為對(duì)應(yīng)的Copula連接函數(shù)值。

2)擬合參數(shù)所需條件不同：常用二元概率分布(密度)函數(shù)的參數(shù)擬合必須是森林結(jié)構(gòu)特征因子實(shí)測(cè)數(shù)據(jù)的聯(lián)合分布(密度)值，而二元Copula分布(密度)函數(shù)的參數(shù)擬合只需森林結(jié)構(gòu)特征因子實(shí)測(cè)數(shù)據(jù)的邊緣分布值，前者參數(shù)擬合所需條件比后者多，從側(cè)面說明二元Copula分布(密度)函數(shù)適用性更廣、實(shí)用性更強(qiáng)。

3)胸徑與年齡邊緣分布函數(shù)的類型：應(yīng)用二元β分布(密度)函數(shù)擬合胸徑年齡聯(lián)合分布時(shí)，要求胸徑和年齡的邊緣分布函數(shù)類型必須為一元β分布(密度)函數(shù)，二元Weibull分布函數(shù)、二元Sbb密度函數(shù)也有類似要求，如果胸徑和年齡邊緣分布函數(shù)類型不同，則常用二元分布(密度)函數(shù)不適用，而二元Copula分布(密度)函數(shù)則適用于任意類型的胸徑和年齡邊緣分布函數(shù)。

4)聯(lián)合分布函數(shù)的建立：常用二元分布(密度)函數(shù)不是由對(duì)應(yīng)的邊緣分布(密度)函數(shù)推導(dǎo)建立的(葛宏立等，2008；劉恩斌等，2010；Hafleyetal.，1976)，而無(wú)論在理論還是在實(shí)踐中，森林結(jié)構(gòu)特征因子邊緣分布的類型眾多且較容易確定，但其聯(lián)合分布的價(jià)值較邊緣分布更大，因此如何根據(jù)森林結(jié)構(gòu)特征因子的邊緣分布確定對(duì)應(yīng)的聯(lián)合分布意義重大。二元Copula分布(密度)函數(shù)通過連接胸徑和年齡的邊緣分布函數(shù)構(gòu)建毛竹胸徑年齡二元聯(lián)合分布(密度)，為邊緣分布推算對(duì)應(yīng)的聯(lián)合分布提供了一種可行方法。

5)模型的檢驗(yàn)：實(shí)測(cè)森林結(jié)構(gòu)特征因子的聯(lián)合分布(密度)值是擬合常用二元概率分布(密度)函數(shù)參數(shù)必需的數(shù)據(jù)，為了檢驗(yàn)實(shí)測(cè)數(shù)據(jù)是否服從指定的分布(密度)函數(shù)，需進(jìn)行柯爾莫哥洛夫或卡方檢驗(yàn)，而二元Copula分布(密度)函數(shù)的參數(shù)擬合只需森林結(jié)構(gòu)特征因子的邊緣分布函數(shù)值[不需要其聯(lián)合分布(密度)值]，且在很多情況下，對(duì)于有些樣本的實(shí)測(cè)聯(lián)合分布(密度)值是得不到的，故AIC常被用來(lái)優(yōu)選與檢驗(yàn)二元Copula分布(密度)函數(shù)。

6)實(shí)現(xiàn)功能的比較：常用二元分布(密度)函數(shù)的功能是由兩森林結(jié)構(gòu)特征因子得到對(duì)應(yīng)的聯(lián)合分布(密度)值，二元Copula密度函數(shù)的功能是由兩森林結(jié)構(gòu)特征因子的分布函數(shù)值得到對(duì)應(yīng)的Copula密度函數(shù)值，故常用二元分布(密度)函數(shù)圖描述的是兩森林結(jié)構(gòu)特征因子的結(jié)構(gòu)特征，而二元Copula密度函數(shù)圖描述的是兩森林結(jié)構(gòu)特征因子的相關(guān)性特征(圖2)。經(jīng)變量相互轉(zhuǎn)換，二元Copula分布(密度)函數(shù)可實(shí)現(xiàn)常用二元分布(密度)函數(shù)的功能，如設(shè)D為胸徑，A為年齡，u=F(D)為胸徑的邊緣分布函數(shù),v=G(A)為年齡的邊緣分布函數(shù)，則由反函數(shù)D=F-1(u),A=G-1(v)可得對(duì)應(yīng)的胸徑和年齡值，由式(1)得胸徑年齡二元Copula聯(lián)合分布值就是毛竹胸徑年齡二元聯(lián)合分布值，與常用二元分布函數(shù)得出的胸徑年齡聯(lián)合分布值類似，以胸徑和年齡為自變量，繪制毛竹胸徑年齡的二元聯(lián)合分布圖如圖5所示。可以看出：隨著胸徑和年齡增大，毛竹胸徑年齡的累積分布值也增大，當(dāng)胸徑為5和15 cm、年齡為1和5度時(shí)，累計(jì)分布值分別取最小值和最大值，與實(shí)際相符。

以胸徑和年齡為自變量，由式(3)得胸徑年齡二元密度函數(shù)h(D,A)，繪制毛竹胸徑年齡的二元聯(lián)合密度圖如圖6所示。可以看出：1)胸徑年齡二元密度函數(shù)值不隨胸徑和年齡取值增大而增大，當(dāng)胸徑為9和15 cm、年齡為3和2度時(shí)，其值為最大和最小；2)雖然對(duì)單株毛竹而言，其胸徑和年齡的相關(guān)性不大，但由于立地條件和人為干擾等因素使省域尺度毛竹胸徑年齡具有某種分布規(guī)律。

圖6 毛竹胸徑年齡的二元聯(lián)合密度圖Fig.6 The bivariate joint density of DBH and age of moso bamboo

3.2 森林結(jié)構(gòu)特征因子的相關(guān)性

相關(guān)系數(shù)和偏相關(guān)系數(shù)常用來(lái)描述兩森林結(jié)構(gòu)特征因子之間的相關(guān)性(范葉青等，2013)，但其只能從線性角度進(jìn)行描述，即兩森林結(jié)構(gòu)特征因子尾部與中間的相關(guān)性相等，當(dāng)森林結(jié)構(gòu)特征因子間的關(guān)系不能用線性函數(shù)描述時(shí)，則上尾部與下尾部的相關(guān)性可能不同，如從圖2可知浙江省毛竹胸徑和年齡的相關(guān)性上尾部較大、下尾部較小且漸趨于0，此時(shí)(偏)相關(guān)系數(shù)得出的結(jié)果是有誤的，這也說明二元正態(tài)Copula密度函數(shù)不能描述毛竹胸徑年齡分布的原因。

二元Gumbel Copula密度函數(shù)是毛竹胸徑年齡的最優(yōu)Copula密度函數(shù)，在此基礎(chǔ)上，得出胸徑和年齡Kendall秩相關(guān)系數(shù)為0.052 1、下尾相關(guān)系數(shù)為0、上尾相關(guān)系數(shù)為0.070 9、Pearson線性相關(guān)系數(shù)為0.000 0。

因毛竹獨(dú)特的生物學(xué)特性，對(duì)單株毛竹而言，其胸徑和年齡互不相關(guān)，但對(duì)省域尺度的毛竹而言，由于受立地因子、頻繁經(jīng)營(yíng)等諸多不確定因素影響，使得省域尺度毛竹的胸徑和年齡結(jié)構(gòu)發(fā)生改變，從而使毛竹胸徑和年齡存在某種復(fù)雜關(guān)系。Pearson相關(guān)系數(shù)的使用前提是假定毛竹胸徑和年齡滿足線性關(guān)系，顯然不能反映胸徑和年齡存在的這種復(fù)雜關(guān)系，故得出毛竹胸徑和年齡不相關(guān)的結(jié)論不可靠。Kendall秩相關(guān)系數(shù)能反映兩變量的變化趨勢(shì)是否一致，故可以度量毛竹胸徑和年齡的非線性復(fù)雜關(guān)系，毛竹胸徑與年齡的Kendall秩相關(guān)系數(shù)為0.052 1，說明毛竹胸徑和年齡的變化趨勢(shì)一致，即隨著毛竹林分年齡增大，林分胸徑也在增大，這從上尾相關(guān)系數(shù)為0.070 9也能得到體現(xiàn)。

綜上所述，二元Copula密度函數(shù)既能描述兩森林結(jié)構(gòu)特征因子的聯(lián)合分布，也能客觀反映兩森林結(jié)構(gòu)特征因子之間的相關(guān)性。

3.3 森林結(jié)構(gòu)特征因子聯(lián)合分布研究展望

本研究應(yīng)用二元Copula分布(密度)函數(shù)探討毛竹胸徑年齡的二元聯(lián)合分布，Copula分布(密度)函數(shù)對(duì)森林結(jié)構(gòu)特征因子聯(lián)合分布具有很大的應(yīng)用價(jià)值，主要體現(xiàn)在以下4方面：1)使森林結(jié)構(gòu)特征因子聯(lián)合分布動(dòng)態(tài)模型的建立成為可能，二元Copula分布(密度)函數(shù)的參數(shù)少，容易建立函數(shù)參數(shù)與時(shí)間t的函數(shù)關(guān)系；2)建立森林結(jié)構(gòu)特征因子聯(lián)合分布的混合概率模型，有些森林結(jié)構(gòu)特征因子因受諸多因素影響，使得因子間存在復(fù)雜的尾部相關(guān)性，而不同二元Copula分布(密度)函數(shù)可以反映不同的尾部相關(guān)性，Copula分布(密度)函數(shù)混合即可全面描述森林結(jié)構(gòu)特征因子間的這種復(fù)雜關(guān)系；3)使森林結(jié)構(gòu)特征因子聯(lián)合分布的研究范圍大大拓寬，由于條件所限，只能獲得兩森林結(jié)構(gòu)特征因子的實(shí)測(cè)邊緣分布值，無(wú)法得到其實(shí)測(cè)聯(lián)合分布(密度)值，在這種情況下，應(yīng)用常用二元分布(密度)函數(shù)無(wú)法建立其聯(lián)合分布模型，但使用二元Copula分布(密度)函數(shù)仍可研究其聯(lián)合分布；4)本研究方法也可建立其他森林結(jié)構(gòu)特征因子的聯(lián)合分布(密度)模型，即，首先根據(jù)實(shí)測(cè)數(shù)據(jù)得森林結(jié)構(gòu)特征因子的邊緣分布值，然后確定最優(yōu)Copula密度函數(shù)并擬合相應(yīng)Copula密度函數(shù)的參數(shù)，最后建立森林結(jié)構(gòu)特征因子的聯(lián)合分布(密度)模型。

4 結(jié)論

1)二元Weibull分布函數(shù)的測(cè)量精度最高，R2=0.990 1，但其參數(shù)較多，迭代不易收斂，且在估算毛竹胸徑年齡聯(lián)合密度值時(shí)需進(jìn)行數(shù)據(jù)相互轉(zhuǎn)換，較為繁瑣；二元 Sbb密度函數(shù)的測(cè)量精度最低，參數(shù)最多，R2=0.736 2；二元Gumbel Copula密度函數(shù)的測(cè)量精度也非常高，R2=0.984 1，只有1個(gè)參數(shù)，迭代過程非常容易收斂，估算聯(lián)合密度時(shí)操作簡(jiǎn)單方便。

2)基于Copula分布(密度)函數(shù)建立毛竹胸徑年齡二元分布模型的方法，適合于胸徑和年齡的任意邊緣分布函數(shù)，具體應(yīng)用時(shí)，不需要確定邊緣分布函數(shù)的類型，由實(shí)測(cè)數(shù)據(jù)得胸徑和年齡的邊緣分布函數(shù)值即Copula分布(密度)函數(shù)的自變量。

3)二元Gumbel Copula密度函數(shù)的 AIC為-19.519 6，為5個(gè)常用Copula密度函數(shù)中AIC最小的密度函數(shù)，是毛竹胸徑年齡聯(lián)合分布的最優(yōu)Copula密度函數(shù)。

4)二元Gumbel Copula密度函數(shù)、二元Weibull分布函數(shù)和二元Sbb密度函數(shù)累計(jì)值的最大偏差分別為 0.015 8、0.007 0、0.078 1，在0.05顯著性水平下的閾值為0.179 8，故毛竹胸徑年齡聯(lián)合分布均服從3種分布(密度)函數(shù)。

5)二元Copula分布(密度)函數(shù)較常用二元分布(密度)函數(shù)有較低的適用條件與較廣的應(yīng)用范圍。