應用MATLAB計算昆蟲生命表矩陣模型

陶士強+吳福安

摘要：內稟增長率是生命表的一個重要參數，它綜合了昆蟲發育、存活、產卵等因素，反映種群在一定環境下的數量增長能力。內稟增長率的計算有近似與精確兩種方法，外文較多的采用精確算法，但較為復雜、不易實現。該文介紹應用MATLAB軟件實現計算內稟增長率的精確值、種群的繁殖力值、穩定的年齡分布的方法，簡單快捷，便于推廣使用。

關鍵詞：內稟增長率；MATLAB；矩陣模型

中圖分類號：TP393 文獻標識碼：A 文章編號：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

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

1 材料與方法

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

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

2 結果與分析

依據轉移矩陣A計算得到的特征值見表1 ，由表1看出：有7個值是數字0，有22個虛數和1個負數（-1.0627），得到唯一的正實數是序號第7的1.4116，依據Caswell描述種群的周限增長率為最大的正實數[9]，該值1.4116即為種群的周限增長率。采用公式rm=log（1.4116）計算得到種群內稟增長率的精確值0.3447。該計算結果同文獻[6] 計算結果一致，說明采用上述計算方法結果可靠。

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

3 討論

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

參考文獻：

[1] 吳福安，周金星，余茂德，等.不同桑樹品種上朱砂葉螨實驗種群內稟增長率的統計推斷[J].昆蟲學報，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] 丁巖欽.昆蟲數學生態學[M].北京：科學出版社，1994.

[4] 陶士強，吳福安.應用Jackknife技術統計昆蟲生命表參數變異的VFP實現[J].昆蟲知識，2006，43（2）：262-265.

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

[6] 陶士強，吳福安，余茂德，等.朱砂葉螨在桑品種育71-1上的實驗種群生命參數研究[J].蠶業科學，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