臺灣水稻褐飛蝨族群發生動態之預測模式

臺灣水稻褐飛蝨族群發生動態之預測模式

賴信宏

1

_、黃守宏

2

_、鄭清煥

2

_、蔣國司

1,

_*

1_{國立中興大學農藝學系} 2_{行政院農委會農業試驗所嘉義分所植物保護系}

摘要

褐飛蝨(Nilaparvata lugens, Stål)為臺灣 最重要水稻害蟲，如防治不當常導致二期作 水稻產量及品質顯著下降。為達到水稻害蟲 經濟安全有效的管理，建立早期預警，作為 適時防治之依據，是為解決此一問題之重要 課題。本研究使用 1988 至 2012 年在嘉義分 所溪口農場所設置誘蟲燈每日捕捉褐飛蝨之 數據，利用 1988-2003 年之數據分析，發現 資料在時間序列上存在高度的自我相關。首 先嘗試運用三情況門檻自我迴歸分析建立模 式，並以近 9 年 2004-2012 的資料來驗證預 測模式之準確性。模式之驗證分為長期與短 期預測，長期預測為進行長時間(約 120-160 日)每日蟲數之預測，可評估二期作水稻生育 期間每日蟲數。結果顯示，族群發生動態與 實際趨勢大部分一致；短期預測之發展為評 估第1 至 14 日後預測值之比較，結果顯示從 預測起始日第 7 日後之預測值具有相當高的 穩健性，說明此模式可有效運用於 7 日後之 族 群 發 生 動 態 趨 勢 。 期 望 能 以 此 研 究 的 結 果，提供未來運用於防治決策的參考。 關鍵詞︰褐飛蝨、族群發生動態、預測、門 檻自我迴歸模式。

Forecasting Model for Brown

Planthopper Population Fluctuation

and Its Effects on Rice Production

in Taiwan

Sin-Hong Lai1_{, Shou-Horng Huang}2_, Ching-Huan Cheng2 _{and Kuo-Szu} Chiang1,_*

1 _{Department of Agronomy, National Chung}

Hsing University, Taichung City 40227, Taiwan ROC

2 _{Department of Plant Protection, Chiayi}

Agricultural Experiment Station, Taiwan Agricultural Research Institute, Chiayi 60044, Taiwan ROC

ABSTRACT

The brown planthopper (Nilaparvata lugens, Stål) is an important pest insect which affects rice production in Taiwan. Incorrect control strategies will reduce the quality and yield of rice production. The brown planthoppers can migrate to Taiwan every year from neighboring areas. Moreover, the immigration time and population abundance are often changeable. By a system of long-term monitoring of the insect, a time-series population fluctuation of brown planthopper can be recorded. The forecasting system for the outbreak time of brown planthoppers will provide early warning and information on chemical application for safe rice production. In this study, population fluctuations based on daily data collected from paddy fields and traps were monitored in Chiayi County from 1988 to 2012. Owing to the autocorrelation of the data, we analyzed them by using the three-regime threshold autoregressive (TAR) statistical model of time series. Firstly, the data from 1988 to 2003 was used to establish the prediction model. Secondly, the data from 2004 to 2012 were employed to test the validity of the predicted model. A long-term forecast provided 120-160 days of prediction after the first prediction date was used to estimate daily forecasting data in the second crop season. Results showed that most of the forecasting trends are near the trends of the * 通信作者, [email protected]

投稿日期：2014 年 4 月 27 日 接受日期：2014 年 5 月 6 日

作物、環境與生物資訊 11:65-79 (2014)

Crop, Environment & Bioinformatics 11:65-79 (2014) 189 Chung-Cheng Rd., Wufeng, Taichung 41362, Taiwan ROC

observed data. For short-term forecasting, we used the results of one-day forecasting to those of fourteen-day forecasting to describe the precision of the forecasting model. The results indicated that the trend of seven-day forecasting is recommended. That is, the forecasting model could effectively estimate population fluctuations seven days in advance. Accordingly, the results of this study are applicable to be included in the plant protection measures.

Key words: Brown planthopper, Nilaparvata lugens,

Population fluctuations, Forecasting, Threshold autoregressive model.

前言

褐飛蝨(Nilaparvata lugens, Stål)為東亞地 區最重要的水稻害蟲，大量族群聚集吸食水稻 汁液，往往在短時間內即可致使水稻植株枯 萎，形成「蝨燒」，嚴重影響水稻產量與稻米品 質。褐飛蝨為長距離遷移性水稻害蟲之一，但 其族群在臺灣可越冬，惟越冬族群密度甚低， 一般在一期作對水稻生產並不構成影響；在二 期稻，由於發生族群密度受7 至 8 月間由境外 一千公里範圍內遷入之數量所左右。其中，以 中國華南、北越及菲律賓呂宋島為主要來源地 區，但每年遷入的蟲源地區、遷入時間及數量， 受氣象條件所影響。褐飛蝨遷入數量、時期及 危害程度，可能因此產生變動，對於臺灣水稻 安全生產之影響帶來潛在威脅。因此，如何有 效地預測及偵(監)測害蟲族群的發生，成為經 濟有效管理褐飛蝨的首要課題(Cheng 1990, Cheng and Huang 2004, Cheng and Lu 1990, Huang et al. 2010a, Huang et al. 2010b, Otuka et al. 2012)。 害蟲預測模式之建立，有賴於對害蟲發 生動態之長期監測，藉由其田間實際發生之 變動，發展有效之模式。在臺灣褐飛蝨於第 一、二期稻都可發生34 個世代，而每 25 至 30 日左右可完成一個世代。其族群立足後快 速 增 長 ， 最 高 族 群 密 度 經 常 發 生 於 第 三 世 代，斯時已屆水稻繁殖生長期，植株老化， 氮碳比顯著下降，所發育之成蟲絕大部分為 長翅型，往外遷徙，導致第一期稻之褐飛蝨 族群顯著下降。挨次期作，於第二期水稻之 生育初期，褐飛蝨長翅型成蟲又陸續由海外 遷入，遷入之蟲數與二期稻期間褐飛蝨之族 群密度密切相關(Cheng 1990, Huang et al. 2009)。過去在臺灣褐飛蝨的預測研究，主要 有利用逐步複迴歸分析方法，利用 7 至 8 月 遷入蟲源，以及逐步分析可能影響發生動態 之各種具生物生態學意義之環境因子，據以 建立預測模式。又擇相關性較高者、經適合 度 試 測 ， 以 準 確 度 較 高 的 環 境 因 子 來 建 立 中 、 短 期 之 發 生 預 測 模 式(Cheng 1984, Cheng 1990)。為了預防褐飛蝨可能造成之經 濟危害損失，Cheng (1978)建議在族群入侵 稻田繁殖後之第 2 世代若蟲期，因其族群密 度達所設定之經濟危害基準密度，約在毎一 叢稻平均達5~10 隻若蟲範圍，宜先予以適時 防治，可以降低與延緩後續族群之建立及可 能造成之危害損失。 在亞洲其他地區水稻害蟲發生預測之 研究，主要是評估每年害蟲族群之發生與其 可能造成之危害損失為主。例如：Ma and Cheng (2001)分析 1986 到 1998 年中國太湖 地區6 月至 11 月每四天之監測褐飛蝨數據， 使 用 混 沌 模 式 (chaos) 與 時 間 規 模 (time-scale)上之關係式來量化褐飛蝨發生系 統。Yan et al. (2010)利用中國貴州地區 1977 年到 2007 年褐飛蝨與白背飛蝨的誘蟲燈數

據，使用馬可夫鏈(Markov chain theory)來 建立預測模式，驗證時採用前五年的數據來 預 測 第 六 年 的 發 生 程 度 ， 結 果 有 效 地 預 測 1982 年到 2007 年的每年發生程度。在日本， Yamamura et al. (2006)提出使用廣義加法模 式(generalized additive model; GAM)與向 量 自 我 迴 歸 分 析 (vector autoregressive analysis; VAR)建立長時間的三種水稻害蟲 (二化螟蟲，Chilo suppressalis；僞黑尾葉蟬， Nephotettix cincticeps 及斑飛蝨，Laodelphax striatellus)之動態模式，探討每年蟲數變化與 氣候變遷的關係。

在生態的相關研究中，因常牽涉到時間 效應，故統計上時間序列的分析方法常被運 用。Lin and Liu (1984)蒐集 1975 至 1982 年 褐飛蝨每週田間數據(僅水稻生殖期間)，使用 線性的時間序列分析的自我迴歸整合移動平 均 (autoregressive integrated moving average; ARIMA)模式來建構族群密度之預 測 模 式 ， 並 利 用 殘 差 分 析 說 明 預 測 之 準 確 性。Su et al. (2003)使用 1994 至 2003 年東方 果實蠅的每旬數據，亦是採用 ARIMA 模式 進行分析，並作為短期預測模式。此外，在 中國大陸有多篇報告使用時間序列分析方法 的 ARIMA 模式，分析族群消長與長期運動 規律，預測褐飛蝨、瘤野螟及栗夜盜等害蟲 之發生預測(Hu and Jiang 1997, Yang et al. 2010)。 當不同時期族群表現出的下降與上升急 遽變動，屬於非對稱型(asymmetry)的週期性 變 動 ， 此 時 應 採 用 非 線 性 時 間 序 列 分 析 方 法。在處理這類的非線性時間序列當中，門 檻 自 我 迴 歸 分 析(threshold autoregressive model; TAR) 常 被 使 用 在 族 群 動 態 的 研 究 中，門檻自我迴歸分析於1980 年由 Tong and Lim 首先提出，針對有限制的循環與週期性 變動發展出量化的方法，以解決具上述生態 特性族群之預測。例如在英國綿羊族群動態 研 究 中 ， 根 據 每 年 的 族 群 大 小 來 找 尋 門 檻 值，然後使用TAR 進行量化建立模式，並依 照 當 年 之 族 群 數 量 來 估 算 隔 年 族 群 之 變 動 (Grenfell et al. 1998)；另外在加拿大山貓族 群之研究中，其族群消長具有10 年之循環變 動，其週期性變動亦與族群豐度有關，使用 TAR 來 量 化 其 族 群 豐 度 (Stenseth et. al. 1998)。 本研究之目的分為兩部份，第一，將分 析褐飛蝨族群長期監測誘蟲燈之每日數據， 以門檻自我迴歸模式來發展臺灣水稻褐飛蝨 誘蟲燈誘捕蟲數之長期預測模式，說明利用 誘蟲燈蟲數所建立的水稻褐飛蝨疫情警報系 統。其中，長期預測在 7 月下旬或 8 月下旬 開始進行，評估在 8 月下旬到水稻二期作收 獲這期間的族群動態，說明這期間誘蟲燈誘 捕蟲數的可能發生趨勢，期望在二期作水稻 移植前有效掌握褐飛蝨接下來每日之族群變 動。第二，短期預測則在8 月下旬開始進行， 預先評估短期之族群發生動態，提供更即時 及準確之資訊，未來供褐飛蝨預警上之參考。

材料與方法

褐飛蝨監測數據 於嘉義縣溪口鄉之農業試驗所嘉義分所 溪口農場(20 ha)，設立吸式誘蟲燈(suction light trap)兩具。吸式誘蟲燈以 30 燭光環型 日光燈(FCL-30W，中國電器，台南市新營區) 為 光 源 ， 燈 下 6 cm 接 吸 入 式 電 風 扇 (CEN-920，聖力儀器，台中市大里區)及蟲網 (32 網目，尼龍，慶聚企業股份有限公司，台 中市清水區)，誘蟲燈之柱高為 1.8 m，蟲網 為雙套式以防捕獲蟲之逃逸。兩盞誘蟲燈相 距約 250 m，各面對數十公頃之稻田。誘蟲 燈之開啟(下午 6 時)及關閉(翌日早上 9 時)均 由定時器所控制，每日於早上 9 時前收集所 誘捕蟲之害蟲，以防止蟲體飛逸逃出。每日 收集所捕獲蟲體(1988 年 1 月 1 日至 2012 年 12 月 31 日)，經烘箱 45℃烘乾後 3-5 日後，以 顯微鏡鏡檢分類記錄所捕獲之害蟲種類及數量 等，以偵(監)測每日褐飛蝨遷入族群之情況。 另於溪口農場內設立 2 塊觀測田，每觀 測田至少0.5 公頃，種植感蟲品種(如臺農 67 號等)，依農民慣行管理方式於水稻生育期中 不施用任何殺蟲劑，記錄每年兩期作水稻移 植後之天數(days after transplanting; DAT) 等 資 訊 ， 以 調 查 水 稻 本 田 中 之 水 稻 生 育 時 期。預測模式建立前將數據分成兩個部份， 第一部分為前 16 年(1988 至 2003 年)之數 據，作為建立模式之數據(training datasets) 來估計模式內之參數；第二部份為最近 9 年 (2004 年至 2012 年)之數據，作為模式驗證之 數據(validating datasets)，來驗證預測模式 之準確性。

褐飛蝨從海外入侵到形成危害，通常在 二期作水稻生育時期，在預測模式內首先建 立兩個解釋變數，分別為二期作水稻移植後 的天數(以X1,t表示，下標t 為時間點 t)與期作 變數(以X2,t表示)，來量化水稻生育期與褐飛 蝨消長之關係。其中X1,t與X2,t兩變數皆進行 等級(level)之轉換，X1,t變數以每10 日為一個 等級，以水稻移植後之天數作轉換(levels = 0, 1, 2, 3,.. , 11)，移植後的第 1 至 10 日等級為 1，移植後的第 11 至 20 日等級為 2，以此類 推，移植後天數超過100 日的等級為 11，收 穫後的等級為0；X2,t變數則以每年1 月 1 日 至二期作移植前等級為0，二期作移植後至年 底(12 月 31 日)等級為 1，可作為模式截距項 之調整。 觀察1988 至 2003 年誘蟲燈誘捕褐飛蝨 之整年數據，每年之蟲數累積情形大多從 7 月開始，才有逐漸上升之趨勢，但每年之總 蟲量(褐飛蝨整年誘捕累積量)與水稻二期作 期間的累積蟲量並非相近，從過去16 年所收 集之資料，大體分為三類：第一，量多之年 份大多具有捕獲高峰，資料顯示有兩個時間 會出現，即 7 月褐飛蝨遷入之高峰與 9 月下 旬盛行之高峰，則該年之總蟲量會高於1,000 隻，這類年份 7 月至 8 月底的累積蟲量會在 100-1,500 隻範圍內，且 9 月至 11 月的累積 蟲量會在450-2,000 隻範圍內；第二是該年僅 有9 月底盛行之高峰，則總蟲量亦會高於 800 隻，此類年份 7 月至 9 月 6 日的累積蟲量會 高於100 隻(1992 與 2003 年即為此類型)，且 10 月至 11 月的累積蟲量亦會高於 700 隻；第 三，該年趨勢皆無上述之兩高峰，該年之總 蟲量會低於 500 隻，此類年份 7 月至 8 月底 的累積蟲量會低於100 隻，10 月至 11 月的累 積蟲量會低於 270 隻。據此將第一類與第二 類定義為「量多」之年份，即以該年7 月至 9 月 6 日的累積蟲量超過 100 隻時，屬於量多 之年份(等級為 1)；第三類定義為「量少」之 年份，即該年 7 月至 9 月 6 日的累積蟲量低 於100 隻時，屬於量少之年份(等級為 0)，將 年份為量多或量少之資訊(年份變數，以X3,t 表示)帶入預測模式中，來幫忙解釋族群動態 之變化。 門檻自我迴歸分析 嘉義分所收集的長期觀測資料的數據， 其中誘蟲燈為每日調查數據，可提供褐飛蝨 每日發生之豐度，利用時間序列分析方法來 處理資料在時間上的自我相關，目的來量化 其發生動態以達到預測之目標，希望能提供 準確之未來每日預測蟲數。依照嘉義分所誘 蟲燈蟲數，每年逐日間之族群變化，可分成 三個結構之族群動態，第一部分的結構在每 年1 月至 6 月之間，即於第一期作水稻期間， 蟲數極少且較零散(大多數為 0)；第二部份的 結構為每年 7 月開始至 8 月下旬有較多的蟲 數且數量變動高低起伏較頻繁；第三部分的 結構為9 月下旬到 10 月下旬之高峰；第一部 分為從一期作水稻生育期至海外少量遷入前 的族群動態(前半年)，後兩部分則是海外遷入 盛期後至二期作水稻收穫前(7 月下旬到 11 月 下旬)之族群動態。數據具有明顯斷層，屬於 非線性(non-linear)之結構，因此採用非線性 的時間序列分析方法來量化之。 藉由已知的每年褐飛蝨族群發生動態， 利用褐飛蝨遷入時間、蟲量與二期作水稻生 育時期等因子，建立門檻自我迴歸模式來量 化褐飛蝨之整年趨勢變化，特別是針對二期 作水稻生育時期之發生變動，進而發展長期 與 短 期 之 預 測 模 式 ； 在 本 研 究 中 模 式 之 建 立，參數之估計及模式之驗證，皆使用 R 語

言進行統計分析(Cryer and Chan 2008)。利 用門檻變數(threshold variable, Zt-d)與門檻 值(c)將數據分成 k 個情況(regime)，d 稱為門 檻延遲參數(threshold delay parameter)，依

照 k 個情況用自我迴歸模式(autoregressive

model; AR)來量化，自我迴歸分析是屬於線 性的時間序列分析方法，但依照門檻變數切 割數據而產生了非線性之結構；自我迴歸係 數之估計，採用最小平方法可獲得可靠的結

果，依照數據來改變門檻變數的形式，最簡 易的門檻變數為 d 天前之數值(yt d ) (Tsay 1989, Fan and Yao 2003) ， 先 以 簡 單 的 3-regime 的 TAR 模式(k = 3)為例： 1 2 3 (1) (1) (1) 0 1 1 (2) (2) (2) 0 1 2 1 (3) (3) (3) 0 2 1 if if if p i t i t t d i p t i t i t t d i p i t i t t d i y y c y y c y c y y c                                     







(1) 其中

y

_t為在時間點 t 上之自變數，上標符號

(1)_{表示lower regime model；上標符號}(2)_表示 middle regime model ； 上 標 符 號(3)_{表 示} upper regime model。門檻延遲參數為yt d 中

的 d ， 自 我 迴 歸 的 order 分 別 為p p p1, 2, 3 (dp)，門檻值為c1和c2，0為模式之截距 項， 1, ,...,2 p為自我迴歸模式之係數，t為 模式中之殘差項，稱之白噪音(white noise)。 為了降低數據高低起伏與急遽高峰之變 異性，誘蟲燈褐飛蝨蟲數數據(Yt)需進行對數 轉換(ytlog (10 Yt0.5))，褐飛蝨逐日蟲數與數 天 前(yt d , d1 ~ 20)之 數 據 作 自 我 相 關 分 析，原始數據的自我相關與時間延遲(d)具高 度相關，無法將時間之因子量化，因此分析 時使用第一階差分( yt ytyt1)來進行選擇 模式，量化時間上之自我相關性，達到分析 時具有穩定性(stationary)之前提假設。 TAR 的分析在進行選模步驟時，將X1,t、 2,t X 與X3,t這三解釋變數設定為因子(factor)， 門檻變數、門檻值與各情況的自我迴歸之順 序(order, p)須同時考量，通常使用 Akaike information criterion (AIC)來 作 決 定 ， 當 AIC 值越低表示配適結果越好，AIC 值之估 式為  



 



  1 AIC , , , s 2ln j , , , 2 1 j p q d s L p q d s k kp vq    _    _ (2) 其中Lj是第 j 個情況的概度函數(likelihood function)，p 是自我迴歸之順序，d 是延遲參 數，q 是解釋變數的個數，k 為情況數，s 和 v 分別為時間序列內變數 Y 與變數 X 的維度 數(Tsay 1998)。但各種變數與參數的組合一 起討論，往往無法快速得到較準確之資訊， 通常會先選取門檻變數(Zt-d)，在此以數天前 之差分(yt d yt d yt d 1, d =1, 2, 3, 4, 5)作 為選擇；接著是門檻值之選擇(

c

₁和

c

₂)，第 一階差分之數據本身具有 67%的 0，在 0 以 下(包含 0)則佔全部數據約有 80.67%，剩下為 超過0 佔全部比例為 19.33%，找尋門檻值之 限制，為三情況內的數據不得過少，至少各 情況需佔全部數據6%以上；此時三情況的順 序(order)皆設定為(p p p1, 2, ) (15, 15, 15)3  ，因 門檻變數介於−1.40 到 1.50 之間，所以設定 1 1.4 c   來找尋

c

2，此時

c

2為−1.20, −1.19…, 1.50，估算其 AIC 值，具最低 AIC 值的

c

₂為 最佳門檻值；再來設定

c

₂



1.5

找尋

c

1，此時 1

c

為−1.40, −1.39, …, 1.30，具最低 AIC 值的 1

c

為最佳門檻值；最後採用p1, 2,...,15來選 擇自我迴歸之順序(p p p1, 2, 3)，亦是比較在各 情況(k)下之 AIC 值。利用建模數據僅能先決 定門檻變數、門檻值與自我迴歸之順序，截 距項、自我迴歸模式之係數與解釋變數之係 數皆非一固定數值，使用驗證數據分析時， 參數估計會隨著數據增加而變動。 預測起始日期 建立模式之後為了要驗證模式之預測能 力，一般模式內的解釋變數需先設定數值， 才能讓預測模式準確估計，解釋變數皆與日 期有關，故需要先估算出預測起始日期來設 定所需之數值。首先在上半年褐飛蝨捕捉到 的數量稀少，且六月下旬一期作水稻剛收穫 完，此時若有褐飛蝨遷入，在無適當寄主的 情況下，大部份族群亦會馬上遷出，因此在 驗證模式時，從7 月 1 日開始進行驗證預測。 每年水稻移植往往在某一時期(日期)，X1,t變 數(移植後的天數)與X 變數(兩期作)分級的2,t 時間點皆以二期作水稻移植前後為準則，驗 證時不論是在移植之前或在移植之後找到預測 起始日期，皆假設已知水稻移植的日期，即兩 變數所設定之數值皆為原本設定之數值。

從 7 月開始執行預測驗證，X3,t(年份變 數)無法在 7 月至 8 月得知該年是否為量多或 量少之年份，需依照建模數據之每月累積量 得知判定方式，預測驗證時採用每年 7 月 1 日之後的褐飛蝨累積量來區分。若從 7 月 1 日至9 月 6 日之累積量未達到 100 隻，則該 年為量少之年份，年份變數等級為0；若超過 (含)100 隻，則該年為量多之年份，年份變數 等級為 1。因此，預測起始日期分為兩種情 形，第一種情形為量多年份，預測起始日期 應為累積量超過100 隻之隔日(此時年份變數 等級為 1)；第二種情形為量少年份，從 7 月 1 日至到 9 月 6 日之累積量低於 100 隻(此時 年份變數等級為0)，則 9 月 7 日為預測起始 日期。一般二期作水稻種植在 7 月下旬至 8 月上旬，當 9 月 7 日開始進行預測，水稻正 值分蘗盛期，得知該年可能是量少之年份， 褐飛蝨族群在田間數量可能是偏低的情況， 可提供在孕穗期至乳熟期防治之參考。 長期預測 預測起始日期確定後，方可進行長期預 測之驗證，例如2004 年的預測起始日期為 7 月 4 日，此時水稻尚未移植，須假設已知水 稻移植日期，從 7 月 4 日開始進行每日之預 測(包含 7 月 4 日)，則 7 月 4 日的預測值為利 用 upper regime model 所進行估算出之數 值，再往下一個日期進行預測，此時依照其 門檻變數的值(yt1yt1yt2)，此數值包含

觀測值與預測值，進行下一個日期預測值之 估計。若yt1低於c1，則使用lower regime model 來進行估計；若yt1介於 c1與 c2之 間，則使用middle regime model 進行估計； 當yt1高於 c2，則用 upper regime model 進行估計。如此反覆執行下一個日期之預測 值估算，預測約120-160 日之數值(從插秧前 30 多日加上水稻生育期約 120 日，共計最多 約160 日)，並估算每 1 日的 95%信賴區間， 這些預測值皆未使用預測起始日期之後的觀 測值。而 2005-2012 年之預測如同上述 2004 年之步驟，模式參數估計會隨著每年數據增 加而有變動，每年皆有一百多日之預測蟲數 資料來進行驗證模式。 驗 證 時 先 比 較 預 測 值 曲 線 與 觀 測 値 曲 線 ， 說 明 其 族 群 趨 勢 變 動 與 關 係 ， 並 利 用 LOWESS (locally weighted scatterplot smoothing) 代 表 平 滑 觀 測 值 曲 線 ， 調 整 LOWESS 平滑曲線的參數為 f，當f設定越 大 越 平 滑 ， 但 越 會 偏 離 實 際 資 料 。 最 後 將 2004-2012 年 的 預 測 值 與 平 滑 的 觀 測 值 (LOWESS) 進 行 迴 歸 分 析 與 一 致 性 分 析 (agreement analysis)，迴歸分析時會估算迴 歸參數，並且與 45 度角的直線(yt1yt)作 一 比 較 ； 一 致 性 分 析 則 使 用 相 關 的 一 些 指 標，如估算兩者相關係數(r)，兩數值差的平 均( ( 1 2) /  1 2 ) ， 兩 變 方 之 比 例 (v v1/ 2)，以及利用 2 2 /( 1/ ) b C  v v 來說明 兩數值之誤差程度(Lin et al. 2002)。 短期預測 短期預測的目的在尋求模式短時間(1-14 日)內蟲口密度之估計，如果預測起始日期時 間為t，yt1為隔天(t+1)的預測值，表示為第

1 期預測值(one step forecasting)；依序yt2

為隔兩天(t+2)的預測值，表示為第 2 期預測 值(two step forecasting)；以此類推，yt14為

隔兩週(t+14)的預測值，表示為第 14 期預測 值(fourteen step forecasting)。在執行模式短 期 預 測 驗 證 時 ， 每 日 預 測 資 料 包 含 yt1,

_{t 2}

y_{ ,…,} y_t14。為了評估模式在短期預測之準

確性，本研究使用如下的指標，將第 1 期預

測值至第 14 期預測值分別與相對應之觀測

值 ， 計 算 均 方 誤 差 值(mean square error, MSE)，公式如下： 





2 , , 1 1 MSEi n t i j t i j j y y n    



 _{, i1, 2,...,14 (3)} 其中i 為第 i 期之預測，j 為驗證數據的樣本 數(j1, 2,...,n)，n 為總樣本數。驗證時根據 量 多 與 量 少 年 份 有 所 區 別 ， 若 為 量 多 年 份 (2004-2007)，為了讓 MSE 計算時在不同情況

下的數據樣本數(n)相同，因此從預測起始日 期之後隔14 日進行每日之預測(2004-2006 年 共140 日，2007 年共 110 日)；若為量少年份 (2008-2012)則是從預測起始日期之後進行每 日之預測(共 120 日)。

結果

預測模式之建立 三 情 況 的 門 檻 自 我 迴 歸 模 式 (three-regime threshold autoregressive model)之建立步驟，首先要決定門檻變數， 使用式(2)來比較不同門檻變數的 AIC 值，結 果以



y

_t1數值最低，作為本模式之門檻變 數。接著是決定門檻值，結果顯示門檻值

c

₁與 2

c

皆在0 的時候具有最低的 AIC 值(Fig. 1)； 水稻二期作生育後期因褐飛蝨族群後續世代 大量繁殖而有捕獲高峰出現，為了解其族群 之變動，找尋

c

2比較重要，則

c

1



0

；各情 況 數 據 內 需 佔 全 數 據 比 例 6%以上，

c

2從 0.301 (86%)至 0.698 (94%)開始找尋，結果顯 示c20.602最佳。最後找尋各情況之自我迴 歸的順序，結果自我迴歸的順序(p p p1, 2, 3)依 序為 15, 14, 11。最終預測模式如下： 15 (1) (1) (1) (1) (1) (1) 0 1 1, 2 2, 3 3, 1 1 14 (2) (2) (2) (2) (2) (2) 0 1 1, 2 2, 3 3, 1 1 11 (3) (3) (3) (3) (3) ( 0 1 1, 2 2, 3 3, 1 if 0 if 0 0.602 i t i t t t t t i t i t i t t t t t i i t i t t t t i y X X X y y y X X X y y X X X                                                       3) 1 if yt 0.602            (4) 其中



₀為模式之截距項， 1, ,...,2 15為自我迴 歸模式之係數，  1, 2, 3為模式中兩個解釋變 數之係數，



t為模式中之殘差項，稱之白噪 音 (white noise) ， 門 檻 值 為 c10 , 2 0.602 c  ， 上 標 符 號(1)_{表 示 lower regime} model (yt10)，為蟲數逐漸減少或沒變化 時 之 族 群 變 動 ； 上 標 符 號(2)_{表 示 middle} regime model (0 yt10.602)，為蟲數增加 較少時之族群變動；上標符號(3)_{表示 upper} regime model (yt10.602)，為蟲數增加較 多時之族群變動。

Fig. 1. Akaike information criterion (AIC) plot of different lower thresholds (c₁ 1.40, 1.39,...,1.30 ) when

2 1.50

c  and different upper thresholds (c₂ 1.20, 1.19,...,1.50 ) when c₁ 1.40 for the autoregressive model. The solid, dashed, and dotted lines represent the AIC values of the lower, middle, and upper regime models, respectively, and the dotdash line represents the AIC values of TAR model.

預測起始日期 在 量 多 年 份(2004-2007) 之 預 測 起 始 日 期，大多在7 月間蟲數累積量就超過 100 隻， 因此可在七月開始進行長短期預測，僅2007 年 7 月捕獲褐飛蝨蟲數，與其他量多年份比 較相對較少，因此預測起始日期相對較晚(8 月 18 日)；在量少之年份(2008-2012)，蟲數 之累積量要到10 月過後才會超過 100 隻，因 此預測起始日期皆為9 月 7 日。 長期預測 長期預測是從預測起始日期找到之後開 始進行驗證，比較在預測起始日期後一百多 天(n)之預測值(yˆt i )與實際觀測值(yt i )的族 群變化情形(i = 1, 2, …, n)，結果顯示驗證數 據之年份中，其中 2004-2007 年之趨勢相近

(Fig. 2、Figs. 3A-C)，這四年預測蟲數之曲線 (以實線表示)與觀測値(以黑圈表示)距離較 短；利用 LOWESS 估算平滑觀測値曲線(以 虛線表示)，其中平滑參數

f

設定為 0.2，較 能把劇烈變動的觀測値平滑，但又不會太過 於失去實際資料的特質，結果發現預測曲線 與平滑曲線的上升時間與幅度皆相近，表示 預測結果與實際趨勢相近，這 4 年預測之曲 線在 8 月會有下降之趨勢，然後到 9 月初開 始有上升之現象，到10 月下旬會達到高峰之 後才開始下降，與實際數值相比在 8 月下旬 至 10 月下旬略有高估之趨勢；並利用預測的 95%信賴區間進行驗證，結果得知觀測値皆落 在95%的信賴區間內(Fig. 2、Figs. 3A-3C)。

2008-2012 年屬於量少之年份，僅預測未 來 120 天之族群可能變化，探討其發生動態 之趨勢，可發現9 月到 10 月中旬數量皆相當 少，而在10 月下旬的預測有略為高估，但整 年趨勢偏低。2008 與 2009 年預測蟲數曲線之 趨勢與前面提到 4 個年份(2004-2007)比較是 偏低的，當與實際蟲數曲線相比，僅有10 月 下旬與 11 月上旬有數天的發生較多之蟲數 外，其他趨勢極為相近(Figs. 2、3D、4A)；

Fig. 2. The transformed light trap catches (white circles) of BPH, the weighted scatter plot smoothing (dashed line) of BPH, the predicted catches of TAR (solid line), and 95% confidence interval of the predicted catches of TAR (double-dashed line) for 160 days after the first prediction date for 2004. BPH: brown planthopper. TAR: threshold autoregressive model.

Fig. 3. The transformed light trap catches (white circles) of BPH, the weighted scatter plot smoothing (dashed line) of BPH, the predicted catches of TAR (solid line), and 95% confidence interval of the predicted catches of TAR (double-dashed line) for 160 days after the first prediction date for 2005-2006 (A-B), 130 days after the first prediction date for 2007 (C), and 110 days after the first prediction date for 2008 (D). BPH: brown planthopper. TAR: threshold autoregressive model.

Fig. 4. The transformed light trap catches (white circles) of BPH, the weighted scatter plot smoothing (dashed line) of BPH, the predicted catches of TAR (solid line), and 95% confidence interval of the predicted catches of TAR (double-dashed line) for 110 days after the first prediction date for 2009-2012 (A-D). BPH: brown planthopper. TAR: threshold autoregressive model.

2011 年 10 月下旬到 11 月上旬有一高峰，此 時接近收穫期，預測曲線已經開始下降，此 時稍有低估現象(Fig. 4C)；2010 與 2012 年則 是高估較多，從9 月上旬到 10 月底，預測曲 線皆在平滑觀測値曲線上方，僅在2012 年 11 月上旬，平滑觀測値曲線高於預測曲線(Figs. 4B、4D)；並探討觀測値是否落在預測的 95% 信賴區間內，結果得知觀測値大多落在 95% 的信賴區間內(Fig. 4)，唯有 2008 年 10 月下 旬到11 月上旬，平滑觀測値曲線有一高峰， 沒在95%信賴區間內(Fig. 3D)。 將2004-2012 年的預測值與 LOWESS 值 進行迴歸分析，再與直線yt1yt做比較， 從散佈圖看出迴歸直線略高於直線yt1yt， 預 測 值 有 高 估 之 現 象 ； 估 得 的 迴 歸 模 式 為 yt₁ 0.2498 0.9853yt lowess_,    _{， 此 模 式 之 斜 率} (0.9853) 與 1 在 統 計 上 並 無 顯 著 之 差 異 (p-value = 0.3023)；且一致性分析中，兩數 值 相 關 係 數(r) 為 0.819， 兩 變 方 之 比 例 為 1.203，表示兩者變動相當接近，兩者差平均 (  ) 為 0.398 ， 兩 數 值 之 誤 差 程 度 (Cb) 為 0.912，數值離 1 很接近，表示兩者差距小。 短期預測 結果顯示在2004-2012 年的第 1 期 MSE 值(MSE1)最小，第 14 期 MSE 最大，而第 2~7 期與第1 期之 MSE 差值低於 0.20 (MSE2-MSE1, MSE3-MSE1, …, MSE7-MSE1)，在 2004 年第 8-11 期 MSE 比第 7 期 MSE 甚至還低，在 2006、2009、2010 與 2012 年第 1 期至第 14 期的MSE 差距較小(Table 1)，在 2009-2010 年第1 期 MSE 與第 7 期 MSE 差距低於 0.02， 在2006 與 2012 年第 1 期 MSE 與第 7 期 MSE 差距低於0.07。以 2004 年為例，利用第 1, 3, 5, 7 期之預測結果來說明短期預測之敏感性 (Figs. 5A-D)，在 2004 年 7 月下旬之後的預 測 值 會 隨 著 觀 測 值 下 降 而 預 測 出 下 降 之 趨 勢，觀測值在 8 月有小的族群高峰，預測值 亦會跟著上升，然後在 9 月上旬開始有上升 之趨勢，第1, 3, 5, 7 期預測值之曲線與觀測 值趨勢皆很接近，並且與觀測值相同在10 月 下旬之後有緩慢下降之趨勢。 對 於 短期 預測 ， 如預 測日 期 越接 近(例 如：前 1-2 日)，從 Table 1 獲知預測的結果 MSE 如預期會較小，即第一期 MSE 值最小， 但此種情形在實際上效益不大，因為並沒有 足夠之時間來作蟲害之預報或警示；同時從 Table 1 亦獲得第 4-7 期的 MSE 值頗為相近， 因此我們使用yt7來進行繪圖，目的評估其與 觀測值配適之程度。當使用第 7 期預測值 (yt7)來作為短期預測評估，在 2005 年第 7 期預測值在 7-10 月皆相當接近觀測值，第 7 期預測值在11 月上旬與觀測值比較其高峰稍 有延遲(Fig. 6A)；在 2006、2009、2010 與 2012 年第7 期預測值與觀測值趨勢相當接近(Figs. 6B、7A、7B 及 7D)，尤其是 2006 年的第 7 期預測值與觀測值的距離極短，顯示這幾年 預測結果相當準確；在 2007 年 10 月中旬與 觀 測 值 相 差 較 大 ， 因 觀 測 値 高 峰 變 動 相 當 大 ， 這 段 期 間 之 外 ， 其 餘 都 相 當 接 近(Fig. 6C)；在 2008 與 2011 年 9 月的蟲數皆相當偏 低(近乎 0)，第 7 期預測結果相當準確，10 月下旬第 7 期預測值比觀測值高峰較晚發生 (Fig. 6D 及 Fig. 7C)。

討論

過去褐飛蝨族群發生變動之研究，使用 逐步複迴歸分析建立各因子與水稻二期作蟲 數之關係(Cheng 1984, Cheng 1990)，實際應 用時僅能獲知水稻二期作之蟲數，皆未使用 到每日數據的資訊，僅以某個時期之蟲數累 積 量 作 為 依 據 。 亦 有 使 用 時 間 序 列 分 析 的 ARIMA 的方法來進行分析，乃數據為水稻移 植後之田間每週數據(Lin and Liu, 1984)，每 年數據時間上不具有連續性之性質。本研究 使用從 1988 至 2012 年之每日誘蟲燈數據， 首次採用門檻自我迴歸分析量化其數據，以 門檻變數與門檻值將隨著時間變化的資料分 成三部份，分別說明在不同情況下之族群發 生動態，並以近9 年(2004 至 2012 年)之數據

Table 1. Mean square errors (MSE) of one day forecasting to those of fourteen day forecasting for short-term predictions in 2004-2012.

Year MSE of one day forecasting to those of fourteen day forecasting 1 i i2 i3 i4 i5 i6 i7 2004 0.194904* _{0.230225 0.269148 0.298262 0.318670 0.323939 0.364981} 2005 0.117492 0.151046 0.186347 0.209596 0.229904 0.253605 0.272849 2006 0.078950 0.097037 0.117670 0.129135 0.137306 0.140860 0.145740 2007 0.081025 0.117864 0.149786 0.193886 0.214362 0.245650 0.278464 2008 0.118968 0.141387 0.165658 0.165900 0.198714 0.191186 0.226479 2009 0.117231 0.118324 0.119287 0.130972 0.137273 0.134737 0.128738 2010 0.175026 0.161581 0.189643 0.195480 0.202758 0.196813 0.195787 2011 0.129613 0.156779 0.192018 0.206594 0.239702 0.233350 0.263574 2012 0.050462 0.062093 0.070193 0.080799 0.073873 0.083349 0.086030 8 i i9 i10 i11 i12 i13 i14 2004 0.340133 0.356793 0.341983 0.340191 0.361230 0.403273 0.433703 2005 0.292777 0.311476 0.327427 0.341298 0.349973 0.366072 0.387296 2006 0.145740 0.153222 0.166526 0.164750 0.170401 0.181880 0.191578 2007 0.298044 0.315536 0.340794 0.351862 0.362580 0.382725 0.409653 2008 0.213491 0.226095 0.225564 0.220125 0.220773 0.213020 0.202997 2009 0.141375 0.129423 0.138685 0.133470 0.134240 0.133883 0.133103 2010 0.201539 0.210239 0.197735 0.202329 0.203599 0.204848 0.198063 2011 0.242625 0.266673 0.266673 0.261636 0.253905 0.245516 0.239759 2012 0.081977 0.085602 0.081917 0.086270 0.081431 0.084829 0.082948

* _{The given value is the average of the predicted values. Also, there are predicted values of 140 and 110}

days for 2004-2006 and 2007-2012, respectively.

Fig. 5. The transformed light trap catches (white circles) of BPH and the weighted scatter plot smoothing (dashed line) of BPH. A: one day forecasting catches of TAR (solid line) for 140 days after the prediction date for 2004; B: three day forecasting; C: five day forecasting; D: seven day forecasting; BPH: brown planthopper; and TAR: threshold autoregressive model.

Fig. 6. The transformed light trap catches (white circles) of BPH, the weighted scatter plot smoothing (dashed line) of BPH, and seven day forecasting catches of TAR (solid line) for 140 days after the prediction date for 2005-2006 (A-B) and for 110 days after the prediction date for 2007-2008 (C-D). BPF: brown planthopper. TAR: threshold autoregressive model.

Fig. 7. The transformed light trap catches (white circles) of BPH, the weighted scatter plot smoothing (dashed line) of BPH, and seven-day forecasting catches of TAR (solid line) for 110 days after the prediction date for 2009-2012 (A-D). BPH: brown planthopper. TAR: threshold autoregressive model.

作驗證。結果顯示，大部份年份具有一致性 之趨勢，甚至在 2004-2006 年可在二期作水 稻移植前預測下半年之發生動態，成功建立 褐飛蝨在誘蟲燈上誘捕蟲數之預測模式，可 在水稻二期作移植前後獲知在整個水稻生育 期的族群發生趨勢。 誘蟲燈數據的突增與突減之現象，可作 為褐飛蝨遷入與遷出之初步參考，例如從在7 月褐飛蝨海外遷入後，數據會有高峰產生， 該遷入褐飛蝨族群一般常經過繁殖 2-3 世代 後才可對水稻造成嚴重危害；褐飛蝨遷入後 第一世代成蟲於 9 月中、下旬出現，其族群 密度對於在水稻幼穗形成期到孕穗期發生之 第二世代若蟲之發生密度具密切相關(Cheng 1990)。因此，在建立模式時加入年份變數， 以 7 月至 9 月 6 日的誘蟲燈累積蟲量是否達 到 100 隻作為閥值，目的來估計水稻二期作 褐飛蝨族群發生趨勢。 Hu et al. (2014)發表之文章，針對褐飛蝨 之生態作了些描述，褐飛蝨遷移行為大多數 在傍晚起飛，少數在黎明起飛，而誘蟲燈數 據大多收集到夜間(趨光誘捕)的成蟲資料，以 誘蟲燈數據變動來說，褐飛蝨的遷移行為其 實是變動的，不論是從外地遷入或是本地遷 出；然而田間數據為日間調查所獲得，其中 包含有不同時期的若蟲與成蟲之數據。在中 國探討長江流域的褐飛蝨大發生情況中，在 大發生與發生較少的年份，田間數據與誘蟲 燈數據在水稻生育後期具有顯著的關連性， 即在田間收集到大量數據的期間，誘蟲燈亦 會收集到大量數據，此情況可推估是從本地 遷出的情況居多。由此可以得知，在水稻分 蘗期之後，誘蟲燈捕獲成蟲數之數據與田間 發生數據之趨勢具密切關連，是故誘蟲燈可 作為反映該區域害蟲之豐度，作為預測水稻 二期作毎日可能發生豐度與日期之依據。 嘉義分所收集的長期觀測資料的數據， 其中誘蟲燈為每日調查數據，可提供褐飛蝨 毎日發生之豐度；田間為毎週目視調查之數 據(每 100 或 120 叢)，提供水稻生育期間的田 間族群密度，大多從水稻二期作水稻分蘗到 孕穗期(9 月下旬至 10 月下旬)才有開始急遽 上升。由於收集的時間間隔與方式不同，田 間調查的族群密度包含有若成蟲，若蟲可能 是數天前才誕生，而成蟲會死去，因此田間 數據雖是每週才調查一次，但可反映數天的 族群變動。當兩種數據比較時，須配合與田 間數據所收集的時間一致來進行分析，將誘 蟲燈毎日數據進行累積轉換成毎週數據，結 果發現在驗證的年份內(2004-2012 年)水稻分 蘗期至糊熟期頗具相關性(r = 0.65)，過去文 獻也有使用此方式來解釋誘蟲燈與田間資料 之關連性(Wang and Yen 1984, Hu et al. 2014)。根據上述所言，雖然誘蟲燈數據無法 完全代表田間族群狀況，惟其確是田間族群 一個重要指標；然而在誘蟲燈所捕獲到長翅 型成蟲的資料，與田間實際發生族群(還包含 短翅型成蟲與若蟲)之間的確切關係，尚待進 一步研究探討。 若使用田間數據建立模式雖較可貼近實 際田間發生情形，但預測單位則會從每日變 成每週，無法獲得時間間隔較小的族群發生 動態之資訊；且田間數據大多從水稻分蘗到 孕穗期才開始有上升的趨勢，在此之前的族 群發生動態資訊不足，例如田間數據在水稻 移植前無法得知褐飛蝨何時遷入與數量等資 訊。因此誘蟲燈數據具有預警之功能，可預 先得知褐飛蝨在該地區的可能族群變動，提 供該地區防治策略之參考。 在亞洲其他地區以時間的觀點來探討水 稻害蟲發生預測研究，大多使用每年數據來 進行分析(Yamamura et al. 2006, Yan et al. 2010)，提供每年發生程度之評估；在探討季 節性的族群變動，可能使用每月或每季之數 據，來評估下一月或下一季發生變動之評估 (Hu and Jiang 1997)，這些評估方法大多只 對未來一個週期(年、季、月)進行預測評估， 僅少數有進行多個週期之預測，原因在於受 限於數據之特性。本研究分析之誘蟲燈數據 屬於每日數據，可提供較多資訊，雖然每年

上半年的蟲數稀少，亦可提供參數估計之參 考，而且能將水稻移植前的褐飛蝨遷入情況 納入模式估計中，幫助我們了解後期褐飛蝨 族群可能發生動態之趨勢。 水稻褐飛蝨防治時期以9 月下旬到 10 月 上、中旬為主，水稻相當於幼穗生長期至乳 熟 期 間 ， 此 時 遭 受 危 害 對 產 量 影 響 最 大 (Cheng 1978, 1979, Cheng and Lu 1990, Cheng and Huang, 2009)。本研究發展出的 預測模式可以提前獲知在水稻二期作的褐飛 蝨 發 生 趨 勢 ， 在 長 期 預 測 的 迴 歸 分 析 結 果 中，斜率參數接近於1，表示預測値與觀測値 趨勢相當接近，因此，此模式之長期預測驗 證結果可得到良好的趨勢量化，期望提供在 水稻二期作生育期間之防治參考；並且針對 水稻分蘗期之後進行短期預測，期望能與田 間實際操作配合，來達到規劃噴灑藥劑之時 機，減少成本與環境汙染等問題。

誌謝

本研究承蒙行政院農業委員會動植物防 疫檢疫局經費補助(100 農科-9.1.2-檢-B2(3)與 101 農科-7.1.1-檢-B1(1))，特致謝忱。

引用文獻

Cheng CH (1978) Determination of the economic injury levels of the rice brown planthopper in Taiwan. I. Control of the brown planthopper at various population density in relation to grain yield ant net return. (in Chinese with English summary) J. Agric. Res. China 27:229-236. Cheng CH (1979) Determination of the

Economic-injury levels of the brown planthopper in Taiwan. II. The population levels of Nilaparvata lugens in relation to yield loss of rice. Natl. Sci. Counc. Mon. 7:1103-1114. Cheng CH (1984) Studies on integrated control of brown planthopper, Nilaparvata lugens (Stål) in Taiwan. (in Chinese) p.149-167. In: Ecology and Control of the Brown Planthopper (YI Chu, ed.). Plant Protection Center, Taichung, Taiwan.

Cheng CH (1990) Studies on population dynamics and forecasting of population abundance of

brown planthopper, Nilaparvata lugens in Chia-nan area. (in Chinese with English summary) Chinese J. Entomol. 10:1-25.

Cheng CH, JL Lu (1990) Detection of the trans-oceanic immigration of rice planthoppers,

Nilaparvata lugens Stål and Sogatella furcifera

Horváth to the southwestern Taiwan and their relative weather-conditions. (in Chinese with English summary) Chinese J. Entomol. 10:301-324.

Cheng CH, SH Huang (2004) Population fluctuations and forecasting of the white-backed planthopper, Sogatella furcifera on rice in Chiayi region, Taiwan. (in Chinese with English abstract) Plant Prot. Bull. 46:-315-332. Cheng CH, SH Huang (2009) Reconsideration of

the rice insect pests control in Taiwan. (in Chinese with English summary) p.65-82. In: Proceeding of Symposium on Achievement and Perspectives of Rice Protection in Taiwan. 160 pp. Published by Department of Plant Protection, Chiayi Agriculture Experimental Station, ARI, Chiayi, Taiwan.

Cryer JD, KS Chan (2008) Threshold Models. p.383-422. In: Time Series Analysis with Applications in R. Springer, New York.

Fan J, Q Yao (2003) Parametric nonlinear time series models. p.125-143 In: Nonlinear Time Series: Nonparametric and Parametric Methods. Springer, New York.

Grenfell BT, K Wilson, BF Finkenstädt, TN Coulson, S Murray, SD Albon, LM Pemberton, TH Clutton-Brock, MJ Grawley (1998) Noise and determinism in synchronized sheep dynamics. Nature 394:674-677.

Hu BH, RZ Jiang (1997) Long-term Occurrence Rules and Forecasting of Crops' Disease and Pest Insects. (in Chinese) China Agricultural Press, Beijing. 193pp.

Hu G, F Lu, BP Zhai, MH Lu, WC Liu, F Zhu, XW Wu, GH Chen, XX Zhang (2014) Outbreaks of the Brown Planthopper Nilaparvata lugens (Stål) in the Yangtze River Delta: Immigration or Local Reproduction? PLoS One 9:1-12. Huang SH, CH Cheng, CN Chen, WJ Wu (2009)

The trend of occurrence and prospective control measures of rice insect pests in Taiwan. (in Chinese with English summary) p.131-147.

In: Proceedings of Symposium on Achievement

and Perspectives of Rice Protection in Taiwan. Chiayi Agricultural Experiment Station,

Taiwan Agricultural Research Institute. Chiayi, Taiwan.

Huang SH, CH Cheng, WJ Wu (2010a) Possible impacts of climate change on rice insect pests and management tactics in Taiwan. (in Chinese with English abstract) Crop Environ. Bioinform. 7:269-279.

Huang SH, CH Cheng, CN Chen, WJ Wu, A Otuka (2010b) Estimating the immigration source of rice planthoppers, Nilaparvata lugens (Stål) and Sogatella furcifera (Horváth) (Homoptera: Delphacidae), in Taiwan. Appl.

Entomol. Zool. 45:521-531.

Lin L, AS Hedayat, B Sinha, M Yang. (2002) Statistical methods in assessing agreement: Models, issues, and tools. J. Amer. Stat. Assoc. 97:257-270.

Lin TL, TH Liu (1984) A forecasting model for brown planthopper population density. (in Chinese with English abstract) Chinese J.

Entomol. 4:77-81.

Ma F, X Cheng (2001) Chaos and predictable time-scale of brown planthopper Nilaparvata

lugens (Stål) occurrence system. J. Asia-Pacific

Entoml. 4:67-74.

Otuka A, SH Huang, S Sanada-Morimura, M Matsumura (2012) Migration analysis of

Nilaparvata lugens (Hemiptera Delphacidae)

from the Philippines to Taiwan under typhoon-induced windy conditions. Appl.

Entomol. Zool. 47:263-271.

Stenseth NC, W Falck, KS Chan, ON Bjornstad, M O’Donoghue, H Tong, R Boonstra, S Boutin, CJ Krebs, NG Yoccoz (1998) From patterns to process: phase and density dependencies in the Canadian lynx cycle. Proc. Natl. Acad. Sci.

USA 95:15430-15435.

Su WY, CN Chen, EY Cheng, YB Hwang (2003)

The geographical distribution and statistical forecasting of Oriental fruit flies in Taiwan. (in Chinese with English abstract) Plant Pro. Bull.

Spec. Publ. New 5 (Proceedings of the

Workshop on Plant Protection Management for Sustainable Development: Technology and New Dimension):67-110.

Tong H, KS Lim (1980) Threshold autoregression, limit cycles and cyclical. J. R. Statist. Soc. B 42:245-292.

Tsay RS (1989) Testing and modeling threshold autoregressive processes. J. Amer. Stat. Assoc. 84:231-240.

Tsay RS (1998) Testing and modeling multivariate threshold models. J. Amer. Stat. Assoc. 93:1188-1202.

Wang SS, FC Yen (1984) Analysis of the data of the brown planthopper [Nilaparvata lugens (Stål.)] with light trap in Taiwan. Chinese J.

Entomol. 4:13-22.

Yamamura K, M Yokozawa, M Nishimori, Y Ueda, T Yokosuka (2006) How to analyze long-term insect population dynamics under climate change: 50-year data of three insect pests in paddy fields. Popul. Ecol. 48:31-48.

Yan XH, H Liu, JJ Wang, ZM Zhao, F Huang, DF Cheng (2010) Population Forecasting forecasting Model model of Nilaparvata lugens and Sogatella furcifera (Homoptera: Delphacidae) Based based on Markov Markov Chain chain Theorytheory. Environ. Entomol. 39:1737-1743. Yang LW, EG Wang, LW Ding, CL Gan, JG Chen, YT Dai (2010) Long-term movement cycles of brown planthopper on late rice in Tiantai county of Zheijing province. (in Chinese with English abstract) Chinese Agric. Sci. Bull. 26:310-313.