月相和溫度對棉紅鈴蟲成蟲種群動態的影響

余 純，許 冬，黃 維，雷明洪，萬 鵬，*

(1.江西財經大學 統計學院，江西 南昌 330013； 2.農業農村部華中作物有害生物綜合治理重點實驗室/農作物病蟲草害防控湖北省重點實驗室/湖北省農業科學院植保土肥研究所，湖北 武漢 430064； 3.湖南省安鄉縣植保植檢站，湖南 安鄉 415600)

昆蟲的發生與生理行為受環境影響，不同種類的昆蟲，對外界環境條件中的溫度、濕度、光照甚至月相等因子的適應性各不相同。如棗鐮翅小卷蛾(Ancylissativa)和二化螟(Chilosuppressalis)在高溫時求偶持續時間縮短[1]，異色瓢蟲在低溫下的飛行能力較強但搜索能力較弱[2]，棉鈴蟲(Helicoverpazea)雌成蟲在高溫或低溫條件下性信息素釋放量均會減少[3]。甜菜夜蛾成蟲壽命隨濕度的提高而增加[4]，濕地松粉蚧若蟲的存活率隨降雨強度的增加顯著下降[5]。長光照能促進異色瓢蟲的求偶交配欲[6]，短日照能減少稻縱卷葉螟的南遷[7]，三化螟在滿月期的撲燈數顯著少于月虧期，而黑尾葉蟬的誘集量在滿月期顯著高于其他月相期[8]。因此，研究環境因子對各種昆蟲生理行為的影響，有助于深入認識該昆蟲與環境之間的互作規律，為其在生產上的治理提供理論依據。

紅鈴蟲為長江流域棉花上的重要害蟲，能嚴重危害棉花的產量與品質。有研究表明，月相和溫度均能顯著影響紅鈴蟲的發生[9]。為此，本研究采用性信息素誘集的方式，監測了湖南安鄉棉區的紅鈴蟲種群動態，并運用非凸懲罰的穩健線性回歸模型，分析了月相與溫度對紅鈴蟲發生量的影響，以期能加深對紅鈴蟲發生規律的認識，為該地區紅鈴蟲的綜合治理提供理論依據。

1 材料與方法

1.1 基于性誘劑的紅鈴蟲種群動態監測

2012—2015年，在紅鈴蟲的盛發期，采用紅鈴蟲信息素誘集器(北京中捷四方生物科技股份有限公司)在湖南省安鄉縣的棉田內進行紅鈴蟲成蟲誘集。按照說明書組裝好誘集器各組件，將其放置于棉田內，使誘集器高于棉花植株0.5 m，每天收集所誘成蟲并記數。為保證誘集效果，每月更換1次誘芯。誘集時間從每年的5月下旬至9月底。

1.2 紅鈴蟲監測點氣象數據的獲取

紅鈴蟲監測點的氣象數據從安鄉縣氣象局獲取。

1.3 數據處理與分析

本研究采用穩健的回歸分析——基于非凸懲罰的線性回歸模型，研究紅鈴蟲種群與月相、溫度之間的關系，以避免觀察數據中極端值或離群值對模型參數估計的影響。在回歸分析前，從3個方面進行考慮：一是完全用陰歷的日期來代替月相，分析不同月相時紅鈴蟲的發生量；二是根據月相，將陰歷的日期拆分成新月、弦月、滿月等時間段，分別賦值，分析其與紅鈴蟲發生量之間的關系；第三是除月相外，還需考慮光線的影響，即夜間光線的強弱并結合月相來對陰歷的每天進行賦值，再由此分析其與紅鈴蟲發生量之間的關系。具體賦值情況如下：雨天，0；陰天，1；多云，2；晴天則依月相的不同而賦值不同，其中新月、殘月(農歷初1、初2、29、30)，3；蛾眉月(農歷初3、初4、27、28)，4；蛾眉月(農歷初5、初6、25、26)，5；弦月(農歷初7、初8、23、24)，6；弦月(農歷初9、初10、21、22)，7；凸月(農歷11、12、19、20)，8；凸月(農歷13、14、17、18)，9；滿月(農歷15、16)，10。

1.3.1 線性回歸模型

對于給定的相互獨立的、具同樣分布特征(independent and identically distributed， iid)的觀測值(xi，yi)，i=1，2，3，…，n。為了研究因變量(y)和自變量(x)之間的關系，通常用以下的線性回歸模型：

y=Xβ+，～N(0，σ2I)。

(1)

模型中y=(y1，…，yn)T，是n個觀測值yi組成的向量，如果有p(p≥ 1)個自變量，X就是維度為n×(p+1)的矩陣。β=(β0，β1，…，βp)T是維度為(p+1)×1的未知系數向量。=(1，…，n)T是隨機的誤差向量，假定其服從多元正態分布時，對于第i個觀測值(xi，yi)，yi和xi滿足如下關系：i。最常用的估計系數向量β的方法是最小二乘法，其工作原理就是使得殘差平方和最小。

Yu等[10]指出，當數據中出現極端值或離群值時，最小二乘法不能準確估計系數向量β。在此情況下，需要運用穩健的線性回歸模型來處理數據并檢測離群值。

1.3.2 基于非凸懲罰的穩健線性回歸模型

She等[11]提出了基于非凸懲罰的穩健線性回歸模型：

y=Xβ+γ+，～N(0，σ2I)。

(2)

這個穩健的線性回歸模型(2)是在傳統的線性回歸模型(1)中引入一個均值漂移參數γ=(γ1，…，γn)T，n是樣本容量。當第i個觀測值為離群值時，γi值為非零；否則，γi為0。模型(2)通過對均值漂移參數γ的準確估計，可達到同時估計參數β和檢測離群值的目的。可通過使如下的目標函數達到最小來估計參數β和γ。

(3)

式(3)中，pλ(|γi|)是關于γi的懲罰函數，λ是調諧參數。

1.3.3 懲罰函數與閾值法則

上文所述目標函數(3)中的懲罰函數可以有多種選擇，比較普遍的有l0范數和l1范數，以及SCAD等懲罰函數[12-15]。每個懲罰函數都有其對應的閾值法則(thresholding rule)，比如l0范數的懲罰函數對應有hard閾值法則，l1范數的懲罰函數對應于soft閾值法則，SCAD的懲罰函數對應于SCAD閾值法則。本研究中我們采用l0范數的懲罰函數：

式中，I(·)為指示函數(indicator function)。l0范數的懲罰函數對應于以下的hard閾值法則：

(4)

式(4)中，Θ為閾值法則函數，ξ為變量，λ為調諧參數。

1.3.4 參數估計的具體算法

1.3.5 調諧參數λ最優取值的選擇

因為懲罰函數和閾值法則中包含調諧參數λ，很顯然λ控制γi的估計值，所以調諧參數λ的取值對準確估計γi起到很重要的作用。經驗的做法是在(λmin，λmax)范圍內取100個λ值，其中，λmin是λ的最小取值，使得向量γ=(γ1，…，γn)T中50%的γi值非零，λmax對應的是λ的最大取值，使得向量γ中所有的γi為0。顯然100個不同的λ值對應100組不同的向量γ=(γ1，…，γn)T的估計值，我們用BIC(bayesian information criterion)的標準來選擇最優的調諧參數λ取值。選定最優的調諧參數λ取值后，其對應的一組γ估計值即為最終的參數γ估計值。

通過對參數γ和β的準確估計，穩健線性回歸模型(2)可以實現同時進行離群值檢測和模型參數估計。

2 結果與分析

2.1 月相對誘蛾量的影響

以月相作為自變量(X)，誘蛾量作為因變量(Y)，分別對2012年、2013年、2014年和2015年的數據做誘蛾量與月相之間的線性回歸分析。同時，為了增加數據的代表性，增大樣本容量，2012—2015年的數據也被合并進行了線性回歸分析，以期擬合出更準確的回歸模型。考慮到誘蛾量的觀測值在(0，1 100)之間，變量取值的范圍和方差太大會給線性回歸的參數估計帶來較大的誤差，故在數據分析時以誘蛾量除以10作為因變量Y的取值。如上文所述，我們建立穩健的線性回歸模型為：Y=β0+β1X+γ+。模型中β0、β1和γ為3個未知的參數，其中β0為截距，即當月相取值為0(下雨天氣)時的平均燈誘蛾量；β1為斜率，即當月相每增加1個單位時，平均誘蛾量的變化。均值漂移參數γ=(γ1，…，γn)T是離群值指示器。如果估計的γi為非零，則第i個觀測值是離群值。除了參數估計，本研究還進行了模型顯著性t檢驗，用于檢驗的2個假設分別為零假設H0:β1=0和備擇假設Ha:β1≠0。

圖1展示了2012年、2013年、2014年和2015年和2012—2015年匯總數據中誘蛾量與月相之間關系的散點圖。為了使散點圖更直觀地呈現模型的特征，各散點圖中添加了線性回歸趨勢線。圖1顯示誘蛾量與月相存在負相關關系，這一共性在各年份和匯總數據中都有明顯的體現。各年份的數據和2012—2015年匯總數據中都存在離群值，檢測出的離群值用星號標記。而且，模型擬合出的回歸趨勢線并沒受離群值的影響，體現出所用線性回歸模型的穩健性。

表1顯示了對于誘蛾量(Y)與月相(X)的穩健線性回歸模型中非零的γ估計值。從中可以看出：2012年的數據中，第100、101、104、105和114個觀測值為離群值。2013年的數據中，第43和45個觀測值為離群值。2014年的數據中，第1、100和101個觀測值為離群值。2015年的數據中，第6個觀測值為離群值。2012—2015年的匯總數據中，第100和320個觀測值為離群值。其他的γ估計值為0，說明其對應的觀測值為正常值。

表2為誘蛾量(Y)與月相(X)的穩健線性回歸分析的參數估計和假設檢驗結果。從表中可以看出，在模型顯著性檢驗中，各年份誘蛾量與月相的線性回歸模型都在5%顯著性水平下表現顯著。在各年份的數據和匯總數據中誘蛾量(Y)與月相(X)之間的線性回歸方程非常接近。而且，由各年份的數據和匯總數據得出的95%置信區間相互重疊，說明誘蛾量(Y)與月相(X)之間的線性關系在不同的年份表現出相同的模式，且月相對誘蛾量的抑制程度在不同的年份沒有表現出顯著區別。以2012—2015年的匯總數據為例，晚上沒有月亮時，平均誘蛾量為169頭。月相在0～10的取值范圍內，月相每增加1個單位，平均誘蛾量下降11頭。

圖中星號代表離群值。下同。The asterisks (*) represent the outliers. The same as bellow.圖1 2012、2013、2014、2015年和2012—2015年匯總數據中誘蛾量與月相之間關系的散點圖Fig.1 Scatter plots of moth yield and lunar phase in year 2012， 2013， 2014， 2015 and 2012-2015 combined data

表1誘蛾量與月相的穩健線性回歸模型參數γ的估計

Table1Estimation of parameterγof robust linear regression model with moth yield and lunar phase

年份Year樣本容量Sample size離群值數Number of outliers非零的γ估計值Estimated nonzero γ20121165γ100=90.47,γ101=68.63,γ104=58.73,γ105=58.93,γ114=75.932013832γ43=51.31,γ45=47.5720141153γ1=40.32,γ100=28.39,γ101=30.0220151221γ6=112.022012—20154362γ100=91.46,γ320=111.69

表2誘蛾量與月相的穩健線性回歸分析結果

Table2Results of robust linear regression analysis for the model with moth yield and lunar phase

年份Year數據量Sample size估計的回歸方程Estimated regression equationβ的95%置信區間95% confidence intervalt檢驗P值P value of t test20122013201420152012-201511683115120436Y=19.107-1.888XY=23.364-1.500XY=13.280-1.083XY=15.639-1.079XY=16.900-1.178X(-3.294,-0.481)(-2.581,-0.418)(-1.844,-0.323)(-2.077,-0.081)(-1.722,-0.634)0.0090.0070.0060.034<0.001

2.2 溫度對誘蛾量的影響

把日平均溫度作為自變量(X)，誘蛾量作為因變量(Y)，分別對2012年、2013年、2014年和2015年的數據以及2012—2015年匯總數據做誘蛾量與日平均溫度之間的穩健線性回歸分析。這部分所建立的線性回歸模型、參數含義和模型顯著性檢驗和2.1節所述一樣。

圖2報告了各年份和2012—2015年匯總數據中誘蛾量與日平均溫度的散點圖。和圖1展示的一樣，在各圖中加入線性回歸趨勢線。圖2顯示誘蛾量與溫度存在負相關的關系，這一共性在各年份和匯總數據中都有明顯的體現。各年份的數據和2012—2015年匯總數據中都存在離群值，檢測出的離群值用星號標記。而且，和前文所述一樣，模型擬合出的回歸直線并沒有受離群值的影響，體現出所用的線性回歸模型的穩健性。

圖2 2012、2013、2014、2015年和2012—2015年匯總數據誘蛾量與日平均溫度的散點圖Fig.2 Scatter plots of moth yield and daily mean temperature in year 2012， 2013， 2014， 2015 and 2012-2015 combined data

表3報告了對于誘蛾量(Y)與溫度(X)的穩健線性回歸模型中非零的γ估計值。從表3可以看出：2012年的數據中，第100、101、104、105和114個觀測值為離群值；2013年的數據中，第43、45、53、54和55個觀測值為離群值；2014年的數據中，第1、99和100個觀測值為離群值；2015年的數據中，第6個觀測值為離群值；2012—2015年的匯總數據中，第100和320個觀測值為離群值。

表4報告了對于誘蛾量(Y)與溫度(X)的穩健線性回歸分析結果。從表4可以看出，各年份和匯總數據擬合的模型在5%顯著性水平下都表現為顯著。在各年份數據中誘蛾量(Y)與溫度(X)之間的線性回歸方程非常接近。而且，由各年份數據和匯總數據得出的β1的95%置信區間相互重疊，說明誘蛾量(Y)與溫度(X)之間的線性關系在不同年份表現出相同的模式，且溫度對誘蛾量的抑制程度在不同年份沒有表現出顯著區別。以2012—2015年的匯總數據為例，當日平均溫度為15～35 ℃時，溫度每增加1 ℃，平均誘蛾量減少10頭。

2.3 月相和溫度對誘蛾量的影響

在分析誘蛾量與單個因素(月相或溫度)之間關系的基礎上，我們還考慮在誘蛾量(Y)與月相(X1)之間的線性回歸模型中加入第2個變量——溫度(X2)。得到如下穩健二元線性回歸模型：Y=β0+β1X1+β2X2+γ+。在對二元線性回歸模型進行顯著性假設檢驗中，月相(X1)在2012、2013、2014和2015年均表現為不顯著；月相和溫度在2012—2015年匯總數據中均表現為顯著。所以，我們只在2012—2015年匯總數據中同時考慮誘蛾量(Y)與月相(X1)和溫度(X2)以建立線性回歸模型。

圖3展示了2012—2015年的匯總數據中誘蛾量(Y)與月相(X1)和日平均溫度(X2)三維散點圖，在圖中2個離群值用星號標注出來。從圖中的長方形趨勢平面來看，月相和溫度都與誘蛾量呈負相關，說明月相和溫度都對誘蛾量有抑制作用。對2012—2015年的匯總數據中誘蛾量(Y)與月相(X1)和日平均溫度(X2)建立穩健的線性回歸模型進行參數估計，非零的估計為：γ100=92.90，γ320=113.49，顯示第100和320個觀測值為離群值。在對模型的顯著性檢驗中，二元線性回歸模型在5%顯著性水平下表現顯著(F檢驗的P值為2.957e-07<0.05)。而且月相(X1)和溫度(X2)2個變量在單獨顯著性檢驗中都表現為顯著(對X1和X2的t檢驗P值分別為0.020 2和0.000 1)。估計出的線性回歸模型方程為Y=39.066 - 0.680X1-0.842X2。當溫度保持恒定時，月相每增加1個單位，平均誘蛾量減少7頭；當月相保持恒定時，日平均溫度每增加1 ℃，平均誘蛾量減少8頭。

表3誘蛾量與溫度的穩健線性回歸模型參數γ的估計

Table3Estimation of parameterγof robust linear regression model with moth yield and temperature

年份Year數據量Sample size離群值數Number of outliers非零的γ估計值Estimated nonzero γ20122013201420152012—20151168311512243655312γ100=94.65, γ101=74.14, γ104=59.22, γ105=55.24, γ114=71γ43=57.37, γ45=53.12, γ53=38.48, γ54=41.09, γ55=48.84 γ1=42.83, γ99=29.82, γ100=29.99γ6=116.78γ100=93.28, γ320=115.05

表4誘蛾量與溫度的穩健線性回歸分析結果

Table4Results of robust linear regression analysis for the model with moth yield and temperature

年份Year樣本容量Sample size估計的回歸方程Estimated regression equationβ的95%置信區間95% confidence intervalt檢驗P值P value of t test20122013201420152012—201511683115122436Y=59.645-1.674XY=36.644-0.6973XY=29.627-0.725XY=40.234-1.004XY=42.484-1.023X(-2.403,-0.944)(-1.208,-0.186)(-1.188,-0.263)(-1.879,-0.129)(-1.422,-0.625)0.00010.00810.00240.0249<0.0001

圖3 2012—2015年匯總數據誘蛾量與月相和日平均溫度的三維散點圖Fig.3 Three dimensional scatter plot of moth yield and lunar phase and daily mean temperature for 2012-2015 combined data

3 結論與討論

盛承發等[9]通過圓形分析表明：從全月來看，紅鈴蟲性誘劑誘蛾高峰集中角無統計意義；從分上、下半月來看，紅鈴蟲的性誘蛾高峰數與月相有一定關系，即望、朔蛾峰均少，上弦(初7、初8)、下弦(22～23日)蛾峰均較多。但他僅分析了紅鈴蟲的蛾峰與月相的關系。本研究將紅鈴蟲全年的發生動態與月相之間的關系進行了分析，研究表明，紅鈴蟲發生量在全年內與月相無顯著相關性，但進一步將天氣因素考慮進去后發現，夜間光線對紅鈴蟲發生量的影響更大，即夜間光線越強，紅鈴蟲的誘蛾量就越小，反之則大。該結果表明，月相對紅鈴蟲成蟲活動的影響主要由夜間光線引起，月相本身(月球引力)并不影響紅鈴蟲的發生。

與光線相比較，溫度對昆蟲發生的影響更大。一般而言，影響昆蟲生長發育的是有效積溫。在其發育起點溫度之上時，隨溫度的上升，昆蟲生長發育速率加快，導致昆蟲的種群數量也隨之上升，但溫度過高也會增加其死亡率。在田間條件下，紅鈴蟲的盛發期為每年的6月中下旬至10月上旬，在此期間，如果天氣呈現出高溫高濕條件，就會促進紅鈴蟲的發生，從而導致其發生量大，為害加重。但從本研究結果來看，溫度為15～35 ℃時，隨溫度的上升，紅鈴蟲成蟲的活動顯著下降，即較高的溫度對紅鈴蟲成蟲活動有抑制作用。由于溫度能影響棉鈴蟲的性信息素分泌，如在22～30 ℃條件下，溫度越低，棉鈴蟲性信息素的分泌量越多[16]。因此，在本研究中，高溫對紅鈴蟲誘蛾量的抑制可能與紅鈴蟲的性信息素分泌有關。

線性回歸模型是研究因變量與2個甚至幾個變量(因子)之間關系的主要統計分析方法，通常采用最小二乘法來估計回歸模型中的參數。但是最小二乘法對離群值很敏感，從而影響數據分析的結果。本研究采用了基于非凸懲罰的線性回歸模型來進行分析，該方法具備同時進行離群值檢測和穩健參數估計的功能。分析結果表明，2012—2015年紅鈴蟲的誘蛾量均存在離群值，顯示其不適合用最小二乘法分析。但在處理離群值時，本研究并沒有全部舍棄，而是根據實際情況進行了分析，比如在2012年的數據中，分析出有離群值5個，分別為第100、101、104、105和114。本研究僅舍棄了第100個數據，即9月7日的誘蛾量。原因是：第一，該數值在本年中最高，顯示其為蛾峰值；第二，其實際值與預測值之間的殘差最大。表明這一天的誘蛾量主要受紅鈴蟲成蟲發生高峰的影響，故本研究舍棄了該值。其他幾個離群值雖然實際值與預測值的殘差較大，但它們低于蛾峰值，且其與前后幾天的誘蛾量差異比較小，表明其主要受環境影響所致，故加以保留。同理，對2013年的數據僅舍棄了第43個，2014年的數據沒有舍棄，2015年的數據舍棄了第6個。這些值都為當年的蛾峰值。經此處理后，各年份的回歸模型t測驗的P值均有顯著提升，進一步證明了文中對離群值處理的合理性。同理，在處理溫度與誘蛾量回歸的離群值時，也僅舍棄了蛾高峰期的數據，其他均予以保留。

紅鈴蟲發生除受月相影響外，受降雨的影響也較大。在一定條件下，紅鈴蟲成蟲的活動與雨量呈正相關，即在較小雨量下，弱光和高濕有利于紅鈴蟲成蟲的活動與產卵；但當雨量過大時，紅鈴蟲成蟲的活動受到抑制[17]。說明紅鈴蟲的發生非常復雜，是多種因素共同作用的結果。如能更系統地研究這些因子對紅鈴蟲發生與行為的影響，將有助于深化對紅鈴蟲發生規律的認識，從而在生產上制定合適的防治對策。