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

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

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

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

1 材料與方法

1.1 基于性誘劑的紅鈴蟲種群動態(tài)監(jiān)測

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

1.2 紅鈴蟲監(jiān)測點氣象數(shù)據(jù)的獲取

紅鈴蟲監(jiān)測點的氣象數(shù)據(jù)從安鄉(xiāng)縣氣象局獲取。

1.3 數(shù)據(jù)處理與分析

本研究采用穩(wěn)健的回歸分析——基于非凸懲罰的線性回歸模型，研究紅鈴蟲種群與月相、溫度之間的關系，以避免觀察數(shù)據(jù)中極端值或離群值對模型參數(shù)估計的影響。在回歸分析前，從3個方面進行考慮：一是完全用陰歷的日期來代替月相，分析不同月相時紅鈴蟲的發(fā)生量；二是根據(jù)月相，將陰歷的日期拆分成新月、弦月、滿月等時間段，分別賦值，分析其與紅鈴蟲發(fā)生量之間的關系；第三是除月相外，還需考慮光線的影響，即夜間光線的強弱并結(jié)合月相來對陰歷的每天進行賦值，再由此分析其與紅鈴蟲發(fā)生量之間的關系。具體賦值情況如下：雨天，0；陰天，1；多云，2；晴天則依月相的不同而賦值不同，其中新月、殘月(農(nóng)歷初1、初2、29、30)，3；蛾眉月(農(nóng)歷初3、初4、27、28)，4；蛾眉月(農(nóng)歷初5、初6、25、26)，5；弦月(農(nóng)歷初7、初8、23、24)，6；弦月(農(nóng)歷初9、初10、21、22)，7；凸月(農(nóng)歷11、12、19、20)，8；凸月(農(nóng)歷13、14、17、18)，9；滿月(農(nóng)歷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的未知系數(shù)向量。=(1，…，n)T是隨機的誤差向量，假定其服從多元正態(tài)分布時，對于第i個觀測值(xi，yi)，yi和xi滿足如下關系：i。最常用的估計系數(shù)向量β的方法是最小二乘法，其工作原理就是使得殘差平方和最小。

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

1.3.2 基于非凸懲罰的穩(wěn)健線性回歸模型

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

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

(2)

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

(3)

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

1.3.3 懲罰函數(shù)與閾值法則

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

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

(4)

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

1.3.4 參數(shù)估計的具體算法

1.3.5 調(diào)諧參數(shù)λ最優(yōu)取值的選擇

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

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

2 結(jié)果與分析

2.1 月相對誘蛾量的影響

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

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

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

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

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

表1誘蛾量與月相的穩(wěn)健線性回歸模型參數(shù)γ的估計

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

年份Year樣本容量Sample size離群值數(shù)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誘蛾量與月相的穩(wěn)健線性回歸分析結(jié)果

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

年份Year數(shù)據(jù)量Sample size估計的回歸方程Estimated regression equationβ的95%置信區(qū)間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年的數(shù)據(jù)以及2012—2015年匯總數(shù)據(jù)做誘蛾量與日平均溫度之間的穩(wěn)健線性回歸分析。這部分所建立的線性回歸模型、參數(shù)含義和模型顯著性檢驗和2.1節(jié)所述一樣。

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

圖2 2012、2013、2014、2015年和2012—2015年匯總數(shù)據(jù)誘蛾量與日平均溫度的散點圖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)的穩(wěn)健線性回歸模型中非零的γ估計值。從表3可以看出：2012年的數(shù)據(jù)中，第100、101、104、105和114個觀測值為離群值；2013年的數(shù)據(jù)中，第43、45、53、54和55個觀測值為離群值；2014年的數(shù)據(jù)中，第1、99和100個觀測值為離群值；2015年的數(shù)據(jù)中，第6個觀測值為離群值；2012—2015年的匯總數(shù)據(jù)中，第100和320個觀測值為離群值。

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

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

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

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

表3誘蛾量與溫度的穩(wěn)健線性回歸模型參數(shù)γ的估計

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

年份Year數(shù)據(jù)量Sample size離群值數(shù)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誘蛾量與溫度的穩(wěn)健線性回歸分析結(jié)果

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

年份Year樣本容量Sample size估計的回歸方程Estimated regression equationβ的95%置信區(qū)間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年匯總數(shù)據(jù)誘蛾量與月相和日平均溫度的三維散點圖Fig.3 Three dimensional scatter plot of moth yield and lunar phase and daily mean temperature for 2012-2015 combined data

3 結(jié)論與討論

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

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

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

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