應(yīng)用MATLAB計(jì)算昆蟲生命表矩陣模型

陶士強(qiáng)+吳福安

摘要：內(nèi)稟增長(zhǎng)率是生命表的一個(gè)重要參數(shù)，它綜合了昆蟲發(fā)育、存活、產(chǎn)卵等因素，反映種群在一定環(huán)境下的數(shù)量增長(zhǎng)能力。內(nèi)稟增長(zhǎng)率的計(jì)算有近似與精確兩種方法，外文較多的采用精確算法，但較為復(fù)雜、不易實(shí)現(xiàn)。該文介紹應(yīng)用MATLAB軟件實(shí)現(xiàn)計(jì)算內(nèi)稟增長(zhǎng)率的精確值、種群的繁殖力值、穩(wěn)定的年齡分布的方法，簡(jiǎn)單快捷，便于推廣使用。

關(guān)鍵詞：內(nèi)稟增長(zhǎng)率；MATLAB；矩陣模型

中圖分類號(hào)：TP393 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2018）01-0257-02

Abstract： The intrinsic rate of increase is an important parameter in the life table， it synthesizes the factors of insect development， survival， oviposition， reflect the increased capacity of population in a certain environment. The calculation of the intrinsic rate of increase have approximate and precise methods， the exact method is more acceptable for foreign literature， but it is more complex. The application of MATLAB software was introduced to realize the precise value ， the intrinsic rate of increase in the calculation is fast and simple， and is convenient to popularize.

Key words： Intrinsic rate of increase； MATLAB； Matrix model

內(nèi)稟增長(zhǎng)率（intrinsic rate of increase，簡(jiǎn)寫rm）是生態(tài)學(xué)研究中生命表的一個(gè)重要參數(shù)，它綜合考慮了昆蟲發(fā)育存活、發(fā)育、產(chǎn)卵等因素，反映種群在一定環(huán)境條件下的數(shù)量增長(zhǎng)能力，因此被廣泛應(yīng)用于種群生態(tài)學(xué)研究[1-2]。計(jì)算rm值，有近似與精確兩種方法。近似值的計(jì)算采用rm=ln（R0）/T，精確值的計(jì)算采用公式用迭代法求出[3]，由于計(jì)算復(fù)雜，一般要編制計(jì)算機(jī)程序來(lái)實(shí)現(xiàn)，如陶士強(qiáng)等采用VFP程序設(shè)計(jì)語(yǔ)言編制了求內(nèi)稟增長(zhǎng)率精確值的程序[4]，但應(yīng)用需有一定的VFP語(yǔ)言基礎(chǔ)，不易推廣。Matlab是一種廣泛應(yīng)用于工程計(jì)算及數(shù)值分析領(lǐng)域的新型高級(jí)語(yǔ)言，經(jīng)過(guò)多年的發(fā)展，現(xiàn)已成為國(guó)際公認(rèn)的最優(yōu)秀的工程應(yīng)用開發(fā)軟件[5]。應(yīng)用MATLAB只需簡(jiǎn)單的命令就能求出內(nèi)稟增長(zhǎng)率精確值，簡(jiǎn)單可靠，值的推廣應(yīng)用。

1 材料與方法

采用朱砂葉螨在桑樹育711品種上的生命表數(shù)據(jù)，詳見文獻(xiàn)[6]。x為以ld為單位時(shí)間間隔，li表示任一個(gè)體在i期間得以存活概率，mi是在i期間內(nèi)平均每雌螨產(chǎn)雌卵數(shù)。根據(jù)文獻(xiàn)[7-9]介紹，采用公式Pi=l（i+1）/l（i）， Fi=mi 分別計(jì)算存活率Pi和生殖率Fi，代入轉(zhuǎn)移矩陣A，采用函數(shù)eig（A）求得矩陣A的特征值，特征值中最大的正實(shí)數(shù)即為周限增長(zhǎng)率（Finite rate of increase）。

然后采用公式rm=log（求出rm精確值。采用[v，d]=eig（A）和[w，d]=eig（A）分別計(jì)算矩陣A的左右特征向量。

2 結(jié)果與分析

依據(jù)轉(zhuǎn)移矩陣A計(jì)算得到的特征值見表1 ，由表1看出：有7個(gè)值是數(shù)字0，有22個(gè)虛數(shù)和1個(gè)負(fù)數(shù)（-1.0627），得到唯一的正實(shí)數(shù)是序號(hào)第7的1.4116，依據(jù)Caswell描述種群的周限增長(zhǎng)率為最大的正實(shí)數(shù)[9]，該值1.4116即為種群的周限增長(zhǎng)率。采用公式rm=log（1.4116）計(jì)算得到種群內(nèi)稟增長(zhǎng)率的精確值0.3447。該計(jì)算結(jié)果同文獻(xiàn)[6] 計(jì)算結(jié)果一致，說(shuō)明采用上述計(jì)算方法結(jié)果可靠。

同時(shí)利用eig函數(shù)可以計(jì)算出對(duì)應(yīng)于特征值1.4116的左右特征向量，左特征向量代表了種群的繁殖力值（Reproductive value），右特征向量代表了種群穩(wěn)定的年齡分布（Stable age distribution）。經(jīng)歸一化處理后可得到種群處于穩(wěn)定狀態(tài)下不同年齡組的分布情況。圖1表示該種群未成熟期所占比重大。繁殖力值從第1天開始逐漸升高，第9天達(dá)到高峰（0.467），然后逐漸降低。

3 討論

種群內(nèi)稟增長(zhǎng)率的計(jì)算，外文采用精確算法的較多，一般采用編制專業(yè)的軟件來(lái)實(shí)現(xiàn) [2，10]，由于軟件對(duì)數(shù)據(jù)格式的特定要求，不易操作，有的還涉及著作版權(quán)問題。本文根據(jù)種群生命表數(shù)據(jù)，構(gòu)建矩陣模型，利用MATLAB的2個(gè)命令就能快速計(jì)算內(nèi)稟增長(zhǎng)率的精確值，有利于推動(dòng)生命表技術(shù)的應(yīng)用，有助于不同處理之間的比較。利用MATLAB中的Eig函數(shù)還可以計(jì)算種群的繁殖力值和穩(wěn)定的年齡分布（圖1），對(duì)于掌握種群的特征有著重要意義。為比較不同種群內(nèi)稟增長(zhǎng)率的差異，必須計(jì)算內(nèi)稟增長(zhǎng)率的變異程度。內(nèi)稟增長(zhǎng)率變異程度的計(jì)算有學(xué)者提出采用jackknife技術(shù)來(lái)實(shí)現(xiàn)[11]，如何利用MATLAB實(shí)現(xiàn)內(nèi)稟增長(zhǎng)率精確值變異程度的快速計(jì)算，有待于進(jìn)一步研究。

參考文獻(xiàn)：

[1] 吳福安，周金星，余茂德，等.不同桑樹品種上朱砂葉螨實(shí)驗(yàn)種群內(nèi)稟增長(zhǎng)率的統(tǒng)計(jì)推斷[J].昆蟲學(xué)報(bào)，2006，49（2）：287-294.

[2] Alinejad M， Kheradmand K， Fathipour Y. Sublethal effects of fenazaquin on life table parameters of the predatory mite Amblyseius swirskii （Acari：Phytoseiidae）[J]. Experimental and Applied Acarology.2014（64）：361-373.endprint

[3] 丁巖欽.昆蟲數(shù)學(xué)生態(tài)學(xué)[M].北京：科學(xué)出版社，1994.

[4] 陶士強(qiáng)，吳福安.應(yīng)用Jackknife技術(shù)統(tǒng)計(jì)昆蟲生命表參數(shù)變異的VFP實(shí)現(xiàn)[J].昆蟲知識(shí)，2006，43（2）：262-265.

[5] 王小玉.MATLAB計(jì)算方法[M].北京：清華大學(xué)出版社，2012.

[6] 陶士強(qiáng)，吳福安，余茂德，等.朱砂葉螨在桑品種育71-1上的實(shí)驗(yàn)種群生命參數(shù)研究[J].蠶業(yè)科學(xué)，2005，31（3）：344-347.

[7] Caswell H. Perturbation analysis of nonlinear matrix population models[J]. Demographic Research，2008，18，59-116.

[8] Devillers J. Ecotoxicology Modeling[M].2009，ISSN：1868-1344.

[9] Caswell H. Reproductive value， the stable stage distribution， and sensitivity of the population growth rate to changes in vital rates[J]. Demographic Research. 2010 ，33 ，531-548.

[10] Vangansbeke D， Schrijver LD， Spranghers T， et al. Alternating temperatures affect life table parameters of Phytoseiulus persimilis，Neoseiulus californicus（Acari： Phytoseiidae） and their prey Tetranychus urtica e（Acari：Tetranychidae）[J]. Experimental and Applied Acarology. 2013（61）：285-298.

[11] Maia ADH， Luiz AJB， Campanhola C. Statistical inference on associated fertility life table parameters using jackknife technique： computational aspects[J]. Journal of Economic Entomology，2000，93（2）：511-518.endprint