狀態反饋脈沖控制的福壽螺-水稻防控數學模型

曾夏萍， 梁志清， 龐國萍， 周澤文

（玉林師范學院 數學與統計學院 廣西高校復雜系統優化與大數據處理重點實驗室，廣西 玉林537000）

0 引言

利用數學理論及方法研究有害物防治已經成為有害物種防治的一個重要內容.許多學者通過建立數學模型在有害生物綜合治理研究方面做了大量工作［1-5］.在有害物的實際防控中，常常需要按有害物種的發展狀態實施防控，當有害物種數量達到某一水平時實施殺害，這種殺害方式可用狀態脈沖微分方程的理論進行研究［6］.有關狀態脈沖微分方程的階1 周期解的研究參見文獻［6 -10］.這些學者的工作主要著眼于研究投放天敵、染病害蟲或病毒的生物防治和噴灑農藥殺蟲劑（或微生物制劑）的害蟲防控策略.

在一些特殊的生態環境，如稻田，福壽螺作為入侵生物，已經在農業生產方面產生了一系列的危害和影響，造成了近年來世界范圍內水稻生產難以估量的損失.稻田福壽螺已成為我國局部水稻產區有害生物的重點防控對象，設法有效防治、控制稻田福壽螺的入侵和危害是一項刻不容緩的重大研究任務.防治福壽螺的同時，必須兼顧對水稻及其他稻田生物的保護，盡量避免使用高毒的化學農藥，這方面的生態學者一致認為針對福壽螺生長于稻田的這種特殊生態環境，福壽螺防控應采用生物防控和物理防控［11-16］.物理防治以放茶麩或金腰箭提取物殺滅螺為主，生物防治以稻田中放養吃螺魚類的鴨子為主，防止福壽螺的過度繁殖.文獻［11 -16］針對福壽螺的防控做了許多工作.但對福壽螺防控的研究工作大多數是從生態角度及對有害物的防控手段等方面展開，屬于定性的研究.用數學模型研究防治福壽螺和對水稻保護作動力學方面的研究目前還比較少.運用生物防控和物理防控綜合防治控制福壽螺過度繁殖的數學模型研究目前也少見.

基于以上的生物學背景，本文利用數學理論及方法研究具有狀態脈沖反饋控制的福壽螺-水稻的生態系統.在福壽螺-水稻的捕食與被捕食種群系統的基礎上，通過在稻田中放養鴨子，結合福壽螺發展的“狀態”脈沖噴灑茶麩水控制福壽螺過度繁殖的方法，建構具有脈沖狀態控制的捕食與被捕食種群模型，利用微分方程幾何理論及后繼函數法對模型進行分析，獲得在稻田中將福壽螺控制在水稻可以承受臨界值范圍之內的生物治理策略.

捕食與被捕食關系是學者一直關注的熱點，具有Holing 功能性反應的捕食與被捕食模型形式如下

其中，g（x）為食餌種群的內稟增長率，dy為捕食者種群的死亡率，φ（x）為捕食者的功能性反應函數.

1 模型的建立

在捕食-被捕食系統的基礎上，考慮具有狀態反饋脈沖控制的福壽螺-水稻的生態系統，合理利用福壽螺凈化水體，同時又防止其過度繁殖破壞水稻生長.因為考慮到2 種或多種方法結合使用效果會更理想，所以在模型（2）中引入福壽螺的生物防治，即在稻田中飼養鴨子，鴨子可以吃掉幼螺；物理防治，即當福壽螺超過一定的閾值時，脈沖噴灑一定量的茶麩水.多種方法配合使用，逐漸達到和保持福壽螺和水稻之間的種群動態平衡，取得持續控制福壽螺的結果，得到如下的具有狀態脈沖控制的福壽螺過度繁殖防控數學模型：

捕食-被捕食模型可用于描述福壽螺-水稻的關系，福壽螺為被捕食種群，水稻為食餌種群，是福壽螺的生活資源.在模型（2）中，x表示水稻在時刻t的密度，y表示福壽螺在時刻t的密度.假設在福壽螺生長初期，使用稻田中飼養鴨子抑制福壽螺的增長，單位時間鴨子吃螺率為α（α∈（0，1）），福壽螺數量達到一定數量時會對水稻產生危害，這個數稱之為臨界值，記為h.當福壽螺數量達到時，采用噴灑茶麩水的脈沖方式降低福壽螺數量，殺死的福壽螺數量與種群數量成正比，記比例系數為p（p∈（0，1））.于是以福壽螺發展的“狀態”，考慮生物殺螺和物理殺螺相組合的福壽螺防控措施，得到狀態依賴的脈沖組合模型.

利用后繼函數研究系統（2）周期解的存在性、唯一性和穩定性，尋找放養鴨子和噴灑茶麩水相結合的福壽螺防控策略.

2 預備知識

定義2.1［6］考慮狀態脈沖微分方程稱由狀態脈沖微分方程（3）所定義的解映射所構成的＂動力學系統＂稱為＂半連續動力系統＂，記為（Ω，f，φ，M），規定系統的映射初始點P不能在脈沖集上，P∈Ω =-M｛x，y｝，φ 為連續映射，φ（M）=N，φ 稱為脈沖映射，這里M｛x，y｝和N｛x，y｝為=｛（x，y）∈R2：x≥0，y≥0｝平面上的直線或曲線.M｛x，y｝稱為脈沖集，N｛x，y｝稱為相集.

定義2.2（?x，?y）為狀態脈沖微分方程的平衡態（奇點），若：

1）若（?x，?y）?M｛x，y｝，有f（?x，?y）=0，g（?x，?y）=0；

2）若（?x，?y）∈M｛x，y｝，有f（?x，?y）= 0，g（?x，?y）=0，且α（?x，?y）=β（?x，?y）=0.

1）π（P，0）=P；

2）π（P，t）對P和t均連續；

3）（π（P，t1），t2）=π（P，t1+t2），單參數π（P，t）稱為P的運動軌道.

定義2.4對狀態脈沖微分方程（3）定義的半連續動力系統，映射f（P，t）：Ω→Ω 稱其為自身映射，它包括2 個部分：

1）記下面方程初值為P的Poincarê 映射為π（P，t），有

如果f（P2，t）∩M｛x，y｝=?，則半連續動力系統初值為P的映射為f（P，t）=π（P，t）.

2）如果存在時刻t1使得f（P，T1）=H∈M｛x，y｝，脈沖映射φ（H）= φ（f（P，T1））=P1∈N｛x，y｝，且f（P2，t）∩M｛x，y｝=?，則半連續動力系統初值為P的映射為f（P，t）=π（P，T1）+f（P1，t）.

定義2.5［6］直線M為脈沖集，直線N為相集，如圖1，在相集直線N上A點坐標為（xA，（1 -p）h），設由點A出發的軌線與脈沖集交于一點A′，點A′的脈沖相點A1在相集N上，坐標為A1（xA1，（1 -p）h），定義A1稱為A的后繼點，稱g（A）=xA1-xA為后繼函數.

圖1 脈沖系統（2）的后繼函數Fig. 1 Successor function of（2）

定義2.6階1 周期解若相集N中存在點H，且存在t1使得f（H，t1）=H′∈M｛x，y｝，而且脈沖映射φ（H′）=φ（f（H，t1））=H∈N，則f（H，t1）稱為階1 周期解，其周期為T1，則軌道稱為階1 環.孤立的階1 環為階1 極限環，如圖2.

圖2 脈沖系統（2）的階1 極限環Fig. 2 Order-1 limit cycles of（2）

定義2.7設Γ=f（P，t）是脈沖系統（2）的階1 周期解.如果對任意的ε ＞0，存在δ ＞0 和t0≥0，使得對任意點P1∈U（P，δ）∩N｛x，y｝，當t＞t0時，有ρ（f（P1，t），Γ）＜ε，ρ為半徑，則稱脈沖系統（2）的階1 周期解Γ為軌道穩定的.

引理2. 1后繼函數是微分方程的連續解π（x，t1）與脈沖連續函數φ（x）的復合，這是2 連續函數的復合函數，故后繼函數是連續的.

引理2.2（連續函數的零點定理） 設g（x）是x∈［a，b］上的連續函數，若g（a）·g（b）＜0，則至少存在一點ξ∈（a，b），使得g（ξ）=0.

引理2.3（開區間套定理） 如果：

引理2.4如果相集N上2 點A（xA，（1 -p）h）和B（xB，（1 -p）h），它們的后繼點分別為A1（xA1，（1 -p）h）和B1（xB1，（1 -p）h）.若A和B兩點的后繼函數g（A）=xA1-xA和g（B）=xB1-xB互為異號，則在A和B之間至少存在1 個階1 周期解.

證明由狀態脈沖微分方程Poincarê 映射的定義，A1（xA1，（1 -p）h）是A（xA，（1 -p）h）的后繼點，故存在t1使

于是，后繼函數g（A）=xA1-xA.

又B1（xB1，（1 -p）h）是B（xB，（1 -p）h）的后繼點，故存在t2使

于是，后繼函數g（B）=xB1-xB.

由引理2.1，至少存在一點H∈（A，B）∈N｛x，y｝使g（H）=0，g（H）=c1-c=0，即H點的后繼點為本身，由定義2.4 知，系統至少存在1 個階1 周期解.

3 無脈沖作用下系統的分析

對于系統（2），在無脈沖作用下的子系統為

經簡單的計算，系統（5）有一個平凡平衡點E1（0，0），系統有唯一的正平衡點E*（x*，y*），其中

下面討論正平衡點E*（x*，y*）的穩定性.

在E*（x*，y*）變分矩陣為

因為TrJ（E*）＜0 恒成立，可知E*（x*，y*）局部穩定.

定理3.1系統（5）的正平衡點E*（x*，y*）局部漸近穩定，系統無閉軌，從而沒有極限環.

證明取B=（xy）-1，則

圖3 系統（5）的系統的相圖Fig. 3 Phase diagram of system（5）

4 脈沖系統的定性分析

在沒有脈沖作用下系統（5）的正平衡點E*（x*，y*）是全局漸近穩定的焦點或結點.下面利用微分方程幾何理論及后續函數方法，按正平衡點E*（x*，y*）的類型研究脈沖系統（2）周期解的存在唯一性和穩定性.

4.1 E*（x*，y*）為焦點

情況1h≤y*.

定理4.1若h≤y*，則脈沖系統（2）存在唯一的階1 周期解，且階1 周期解是軌道穩定.

證明當h≤y*時，直線y=h在y=y*下方，如圖4.

這是脈沖集在直角坐標系中為y=h軸在x軸右側的集合.

這是脈沖集在直角坐標系中為y=（1 -p）h軸在x軸右側的集合.

先證階1 周期解的存在性.對N上的點A，過A的軌線交y=h于A′，脈沖到A1∈N，A1為A的后繼點，且xA1＞xA，后繼函數g（A）=xA1-xA＜0.

同理，圖4 中B1為B的后繼點，且xB1＜xB，后繼函數g（B）=xB1-xB＜0.

由引理2.3得，在A與B之間?H，使g（H）=0.H點對應的坐標為（xH，（1-p）h），此時過H的軌線交y=h后脈沖到H點，說明系統（2）存在階1周期解.

圖4 對h≤y*時，脈沖系統（2）階1 周期解的存在性Fig. 4 Existence of order-1 periodic solutions of system（2）for h≤y*

再證唯一性.在A與B之間任意取兩點I，J，不妨設I＜J，則xA＜xI＜xJ＜xB.過I的軌線交直線y=h于I1，脈沖后交直線y=（1-p）h于I+.過J的軌線交直線y=h于J1，脈沖后交y=（1-p）h于J+，如圖5.

圖5 對h≤y*時脈沖系統（2）階1 周期解的唯一性Fig. 5 Uniqueness of order-1 periodic solutions of system（2）for h≤y*

由xI＜xJ得xI+＜xJ+.對后繼函數g有

從而g（I）單調遞減函數，故函數g（I）的根存在唯一，進而脈沖系統（2）的周期解存在唯一.

下面證階1 周期解得軌道穩定性.

設從A1出發的軌線交直線y=h于點A2，脈沖后到直線y=（1 -p）h的點A3，從A3出發的軌線交直線y=h于A4，脈沖后到直線y=（1 -p）h的A5，如此繼續下去得到2 個序列：

同理，從點B1出發的軌線交直線y=h于點B2，脈沖后到直線y=（1 -p）h的點B3，從B3出發的軌線交直線y=h于B4，脈沖后到直線y=（1 -p）h的點B5，如此繼續下去得到2 個序列，見圖6.

上面的序列滿足如下條件：

于是

圖6 序列｛A2n-1｝，｛A2n｝，｛B2n-1｝，｛B2n｝Fig. 6 Sequences｛A2n-1｝，｛A2n｝，｛B2n-1｝，｛B2n｝

在H的左邊附近任取一點Q0（xQ0，（1-p）h），不妨設Q0在A1與H之間，從而xA1＜xQ0＜xH.由｛xA2n+1｝的單調遞增性得：存在n0，使得xA2n0+1＜xQ0＜xA2n0+3，從Q0出發的軌線交脈沖集于Q′0，相點為Q1（xQ1，（1-p）h），易知相點Q1在A2n0+3到A2n0+5之間，從而xA2n0+3＜xQ′0＜xA2n0+5.從Q1出發的軌線交脈沖集于Q′1，相點為Q2（xQ2，（1-p）h），易知相點Q2在A2n0+5到A2n0+7之間，從而xA2n0+5＜xQ′2＜xA2n0+7.如此進行下去得一個點列｛Qk｝，k=0，1，2，…，見圖7，滿足以下條件：

在H的右邊附近任取一點P0（xP0，（1 -p）h）.不妨設P0在H與B1之間，從而xH＜xP0＜xB1.由｛yB2n+1｝的單調遞增性得：存在n0，使得xB2n0+3＜xP0＜xB2n0+1，從P0出發的軌線交脈沖集于，相點為P1（xP1，（1 -p）h），易知相點P1在B2n0+5到B2n0+3之間，從而xB2n0+5＜xPQ0＜xB2n0+3.

從P1出發的軌線交脈沖集于P′1，相點為P2（xP2，（1-p）h），易知相點P2在B2n0+7到B2n0+5之間，從而xB2n0+7＜xP′2＜xB2n0+5.如此進行下去得一個點列｛Pk｝，k=0，1，2，…，見圖7，滿足以下條件：

由兩邊夾定理得

圖7 軌線走向及序列｛Qk｝、｛Pk｝Fig. 7 Sequences｛Qk｝，｛Pk｝

綜上可得，由P、Q的任意性知階1 周期解是軌線穩定的.證畢.

情況2h＞y*，有如下結論.

設直線y=h與過點A的軌線LA相交的一個交點為Ah，記從Ah至A的軌線段為LA hA，由直線y=h、x軸、y軸及LA hA圍成的區域為G，如圖8.

圖8 區域GFig. 8 Region G

定理4.2若h＞y*，則?p*∈（0，1），使得：

1）當0 ＜p＜p*時，軌線最終趨向于E*（x*，y*）；

2）當p≥p*時，脈沖系統（2）存在唯一的階1周期解，且是軌道穩定的.

證明先證p*∈（0，1）的存在性.

圖9 p*的存在性Fig. 9 The existence of p*

從A′出發的軌線LA達到A，脈沖后又回到A′，從而形成周期解.記過A′且平行x軸的直線LA′：y=hA′.由函數y=（1 -p）h關于P的單調連續性，存在?p*∈（0，1），使hA′=（1 -p*）h.

當p＜p*時，直線y=（1 -p）h在y=（1 -p*）h的上方，此時從初值出發的軌線與y=（1 -p）h可能相交也可能不交相.若軌線與y=（1 -p）h沒有相交，即沒有發生脈沖，軌線最終趨向于E*（x*，y*）.若軌線與y=（1 -p）h相交，即發生脈沖，則最后一次與x=xq的交點在A′與A之間，軌線最終趨向于E*（x*，y*），如圖10.

圖10 當p ＜p*時，軌線趨于E*Fig. 10 Trajectory tendency to E*for p ＜p*

當p=p*時，在脈沖的作用下軌線LA在直線y=（1 -p*）h與y=h之間的部分與AA′構成一個階1 周期解.當p＞p*時，從初始點（x（0，y（0））∈G0出發的軌線與脈沖集的交點在點A右側.

在圖11 中，從C出發的軌線交脈沖集于點C′，相點為C1，C1為C的后繼點，且xC1＞xC，后繼函數f（C）=xC1-xC＞0.

同理，圖11 中的B1為B的后繼點，且xB1＜xB，后繼函數f（B）=xB1-xB＜0.

由引理2.3 知脈沖系統（2）存在階1 周期解存在.

圖11 當p ＞p*時，系統（2）階1 周期解的存在性Fig. 11 Existence of order-1 periodic solutions of system（2）for p ＞p*

用定理4.1 的方法可證得階1 周期解的唯一性和穩定性.

4.2 E*（x*，y*）為結點

情況1h≤y*，有如下結論.

定理4.3如果h≤y*，則脈沖系統（2）存在唯一的階1 周期解，且階1 周期解是軌道穩定的.

證明當h≤y*時，

記脈沖集

記相集

如圖12，過A′作y=（1 -p）h的垂線，垂足為A1，記A1坐標為（xA1，（1 -p）h）.

圖12 E*為結點時，系統（2）階1 周期解的存在性Fig. 12 Existence of order-1 periodic solutions of system（2）for the node of E*

過A1的軌線L1交y=h于A2，脈沖后交y=（1 -p）h于A3，記A3坐標為（xA3，（1 -p）h），則A3為A1的后繼點，且xA3＞xA1.于是，后繼函數f（A1）=xA3-xA1＞0.

過B的軌線L2交y=h于B′，脈沖后交y=（1 -p）h于B1，坐標為（xB1，（1 -p）h），則B1為B的后繼點，且xB1＜xB，于是，后繼函數

由連續函數的介值定理，在點A1與B之間存在點H，使f（H）=0.H點對應的坐標為（xH，（1 -p）h）.過H的軌線交y=h后脈沖到H點，說明階1 周期解存在.

用定理4.1 的方法可證得階1 周期解的唯一性和穩定性.

情況2h＞y*.

因E*（x*，y*）是結點，故軌線至多經過若干次脈沖后趨于E*（x*，y*）.如圖13.

圖13 結點的軌線趨勢Fig. 13 Trajectory tendency for the node

5 數值模擬及其生物意義

對脈沖系統（2）進行數值模擬

情形1 取b=1，a=0.4，c=0.1，δ =0.5，β =0.5，則系統的正平衡點E*是穩定的焦點，其中x*=0.5，y*=2.7.

1）取h=2.5，p=0.3，則h≤y*，脈沖系統（2）存在階1 周期解，如圖14.

圖14 對h≤y*，脈沖系統（2）的相位圖及時間序列Fig. 14 Phase portrait and time series of system（2）for h≤y*

其中，x（0）=0.893 5，y（0）=2.5.

2）取h=3.32，則h＞y*，存在p*，p=p*=0.475時脈沖系統（2）存在階1 周期解，如圖15.

圖15 當p=p* =0.475 時，脈沖系統（2）的相位圖及時間序列Fig. 15 Phase portrait and time series of system（2）for p=p* =0.475

其中，x（0）=0.5，y（0）=3.32.

3）取p=0.65，則p＞p*，脈沖系統（2）存在階1 周期解，如圖16.

圖16 對p ＞p*，脈沖系統（2）的相位圖及時間序列Fig. 16 Phase portrait and time series of system（2）for p ＞p*

其中，x（0）=0.87，y（0）=3.32.

4）取p=0.3，則p＜p*，軌線至多脈沖若干次最終趨向E*，如圖17.其中，x（0）=0.85，y（0）=2.

圖17 對p ＜p*，脈沖系統（2）的相位圖及時間序列Fig. 17 Phase portrait and time series of system（2）for p ＜p*

情形2 取b=0.5，a=1，c=1，δ =0.15，β =0.5系統的正平衡點E*是穩定的結點，其中x*=0.5，y*=1.

1）取h=0.9，p=0.3，則h≤y*，系統（2）存在唯一的階1 周期解，如圖18.

圖18 對h≤y*，脈沖系統（2）的相位圖及時間序列Fig. 18 Phase portrait and time series of system（2）for h≤y*

其中，x（0）=0.7，y（0）=0.9.

2）取h=1.2，p=0.4，則h＞y*.

3）系統（2）軌線至多脈沖若干次最終趨向E*，如圖19.

圖19 對h ＞y*，脈沖系統（2）的相位圖及時間序列Fig. 19 Phase portrait and time series of system（2）for h ＞y*

其中，x（0）=1.3，y（0）=1.15.

4）系統（2）沒有脈沖，軌線最終趨向E*，如圖20.

圖20 對h ＞y*，脈沖系統（2）的相位圖及時間序列Fig. 20 Phase portrait and time series of system（2）for h ＞y*

其中，x（0）=0.6，y（0）=0.6.

6 結論及生物意義

1）無脈沖系統有全局漸進穩定的正平衡點E*（x*，y*），其中E*滿足條件

2）周期T的長度由y=（1 -p）h與y=h的寬度λ決定，p越大，λ越大，階1 周期解越長.用后繼函數的方法獲得系統階1 周期解的存在性、唯一性和軌道漸近穩定性，只要h≤y*，脈沖系統（2）均存在唯一階1 周期解且軌道穩定.特別地，當正平衡點E*是穩定的焦點時，存在p*，當p＞p*時，脈沖系統存在唯一階1 周期解.脈沖系統的階1 周期解是軌道穩定的.階1 周期解的周期T隨著p增大而增大，如圖22.

結論表明：在稻田中插一些毛竹引誘成年螺集中附著產卵，然后對卵進行集中銷毀，噴射茶麩水可有效防控福壽螺，茶麩水可通過脈沖控制噴射強度來控制殺螺的時間間隔.

圖21 平衡點與β 的關系Fig. 21 Relation between equilibrium and β

圖22 T與p的關系Fig. 22 Relation between T and p

在實際應用中，根據福壽螺的生長規律，觀測和記錄福壽螺的數量，采取不同的管理策略，從理論上可預測周期時間，從而按周期時間防治福壽螺.

致謝玉林師范學院重點項目（2015YJZD02、G20192K18）和廣西高校復雜系統優化與大數據處理重點實驗室項目（2017CSOBDP0201、2017CSOBDP0202、2016CSOBDP0201）對本文給予了支持，謹致謝意.