基于系統動力學模型的淡水漁業捕撈時間和喂養策略的優化分析

劉梅　原居林　倪蒙　張雷鳴　顧志敏

摘要：從系統動力學角度出發，著眼于動態系統的整體統籌，綜合考慮環境、生態和經濟3方面要素，利用系統模擬軟件Stella 9.1.3構建我國內陸淡水漁業池塘養殖系統動力學模型，再結合Berkeley Madonna優化軟件，以水產養殖戶總利潤最大化為目標函數對模型進行優化。獲取最佳捕撈時間和最優喂養方案，并反推此時對應的最佳養殖容量，同時利用擾動法對所建模型進行參數敏感性分析，以進一步有效提高養殖容量，指導養殖生產，實現淡水漁業經濟生態互利共贏局面，以期實現我國內陸淡水養殖業的可持續發展。

關鍵詞：最佳捕撈時間;最優喂養方案;池塘養殖系統;系統動力學;Stella模型軟件

中圖分類號：S964

文獻標志碼： A

文章編號：1002-1302（2020）02-0164-06

收稿日期：2018-08-17

作者簡介：劉 梅（1988—），女，安徽臨泉人，博士，助理研究員，主要從事養殖水環境研究。E-mail：liumei@zju.edu.cn。

通信作者：顧志敏，研究員，主要從事水產養殖研究。Tel：（0572）2043911;E-mail：guzhimin2006@163.com。

我國是世界淡水漁業養殖大國，養殖產量已連續多年穩居世界第一[1]，淡水漁業總產量約占國內漁業產量的40%[2-3]。除了內陸湖泊、水庫等天然養殖區域外，池塘養殖在我國淡水漁業中占據主導地位。但隨著池塘養殖業的快速發展，養殖區的環境污染問題也越來越受到人們的關注[4-6]。由于養殖過程中片面追求高產，結果不僅導致養殖水體對周邊水體富營養化影響日趨加劇，而且造成養殖對象免疫力下降、病害頻發，造成巨大經濟損失[7-8]。這些問題都嚴重影響我國水產養殖業健康可持續發展。因此，為了緩解養殖水體水質污染，保證淡水水產養殖業的可持續發展，一方面需要優化養殖結構，另一方面須加強對養殖容量的研究。

關于養殖容量的概念，不同研究領域的學者有不同的理解，目前尚無完整統一的定義。1988年，Carver等將貝類養殖的養殖容量定義為產量最大化的同時對生長率不產生負面影響的放養密度，該定義只考慮產量，未考慮生態環境因素的影響，具有明顯的缺陷[9]。董雙林等把養殖容量定義為單位水體內在保護環境、節約資源和保證應有效益的各個方面都符合可持續發展要求的最大養殖量[10]。楊紅生等把淺海貝類養殖業的經濟效益與社會效益、生態效益結合起來，將養殖容量定義為對養殖海區的環境不會造成不利影響，又能保證養殖業可持續發展并有最大效益的最適產量[11]。劉劍昭等將養殖容量定義為特定的水域、單位水體養殖對象在不危害環境、保持生態系統相對穩定、保證經濟效益最大，并且符合可持續發展要求條件下的最大產量[12]。可見，養殖容量的含義在不斷地被充實、豐富，越來越趨于完善。

魚塘生態系統是一種復雜的人工生態系統，受人類控制和影響非常大。若想實現單位水體養殖對象在不危害環境、保持生態系統相對穩定、保證經濟效益最大的同時，符合可持續發展要求條件下的最大產量，只有控制魚塘參數達到最佳時才能實現，如捕撈時間及餌料的輸入等。可以看出，魚塘生態系統涉及環境、生態和經濟社會等多學科多領域，是一個復雜的系統，具有多層次、高階非線性、動態性、自反饋性等特征[8，13]，一般數學方法難以對其進行量化描述和分析。鑒于系統動力學建模軟件（structure thinking experimental learning laboratory with animation，Stella）具有處理動態性、非線性和高階次復雜問題的功能，運用該軟件建立綜合考慮環境、生態和經濟之間相互作用的養殖生態系統動力學模型，然后以水產養殖總利潤最大化為目標函數，利用Berkeley Madonna（http：//www.berkeleymadonna.com）優化Stella所建池塘動態系統模型，以確定最佳捕撈時間及最優喂養策略，并反推此時對應的最佳養殖容量。然后推廣到實際生產中，實現經濟生態互利共贏局面，以期實現我國內陸淡水養殖業的可持續發展。

1 建模工具介紹

系統動力學（system dynamics，簡稱SD）于1956年由美國麻省理工學院福瑞斯特教授開發創建。本研究采用Isee Systems公司的一款可視化系統動力學模擬軟件Stella 9.0.1作為模型構建的平臺，即圖形導向的系統動力學模型發展軟件，它與Vinsim以及Matlab的Simulink軟件包相似，可以提供一個實現動態系統建模、仿真的集成環境[14-15]，動態展現和表達復雜系統和概念實際是如何進行運作的。Stella在研究復雜系統的行為和在處理高度非線性、高階次、多變量、多重反饋問題方面具有優勢，隨著系統動力學理論和方法的不斷發展與完善以及計算機技術的改進，目前已被廣泛應用于規劃學、生物學、生態學和環境科學等領域[16-18]。但是Stella軟件并不能對所建動態系統模型的參數進行優化，以使系統達到所期望的最佳狀態。

而Berkeley Madonna是一款優秀的動態系統建模分析優化軟件，且與Stella兼容，可以運行優化Stella所建模擬系統。由著名的美國加利福尼亞大學伯克利分校研究人員開發，具體介紹詳見網址http：//www.berkeleymadonna.com，其算法優秀，可以在幾秒鐘之內完成數百萬個方程式的計算，可做參數校準（曲線擬合）和優化，有更多解決常微分方程問題的數值方法。因此，本研究利用該軟件對所建魚塘生態系統模型參數進行優化，以方便相關管理部門以及養殖戶選取最佳養殖策略，獲取經濟生態共贏的局面。

2 魚塘生態系統模型的構建

構建系統模型的指導思想是用盡可能簡單的模型來描述系統結構，而構造高層結構圖的指導思想則是簡單、概括，即抓住系統的關鍵性流程，建立系統內的高層關系。由于魚塘是一個匯聚了生物、社會、生態、經濟等各個領域交互作用的動態系統，其包含的因素較多，為便于分析，假設這是一個完全人工飼養的池塘，魚類生存的一個重要條件是在池塘中有足夠的溶解氧。池塘中溶解氧的消耗主要分為2個部分，一部分是魚的呼吸消耗，另一部分則是被用來分解死有機物質即碎屑[19]。碎屑是池塘生態系統中一個非常重要的因素，隨著其濃度的增加，就會出現缺氧甚至厭氧條件[20]。當氧氣濃度低于一定閾值時，魚就會死光，因此本研究假設耗氧量隨著碎屑濃度的增加而增加，簡單地認為魚的死亡是由高碎屑濃度導致的。基于以上假設分析，可以將魚塘生態系統簡單地劃分為社會子系統、經濟子系統、環境子系統3個子系統，每個子系統都有自己的結構特點和獨特功能，其中1個子系統的輸出是其他子系統的輸入，子系統之間彼此聯系，構成了模型的反饋結構，各子系統的構成及相互關系分析具體見下文。

2.1 系統結構圖分析

建立系統動力學模型需要4個基本構造塊：棧（stock）、流（flow）、轉換器（converter）和連接器（connector）。棧用方框表示，代表事物的積累，即用來代表存儲進入和流出棧的物質;流用帶閥門的箭頭表示，用來描述系統中的活動，連接到棧上的流會引起棧中存量的增加或減少，流上的箭頭表示物質或非物質正向流動方向;轉換器用圓圈表示，用以存儲常量，并把它傳輸到其他模塊中為其他變量所用;連接器用1條帶有箭頭的線段或虛線表示，作用是連接模型中的元素，與流不同的是，連接器傳送的是信息。

由于研究區域生態問題的復雜性，因此，在建立系統流程圖時，選擇影響庫區生態安全的核心因子，最終確定參數，建立系統動力學模型的系統流程，具體如圖1所示。

2.2 系統分析及參數的確立

養殖塘生態系統主要有3個子系統，分別是池塘養殖子系統、社會子系統以及經濟子系統，每個子系統的具體控制參數如下所列。

（1）池塘養殖子系統：池塘漁業資源生物量（Fish）與放漁量（Number）、生長率（Growth）死亡率（Mortality）、魚質量（Weight）有關。其中，死亡率導致漁業資源生物量減少，而生長率及放漁量則構成了漁業資源量的補充量。模型中設置漁業資源生物量（Fish）和放漁量（Number）為狀態變量，生長率（Growth）、死亡率（Mortality）為相應的流率變量。

（2）社會子系統：模型的生態部分將設置餌料生物資源（Feed）和碎屑量（Detritus）為狀態變量，餌料喂養量（Feeding）、餌料損失量（Loss）和碎屑分解量（Decomp）為流率變量。

（3）經濟子系統：模型該部分中，設置總利潤Total_profit為狀態變量。流率變量成本（Cost）、餌料價格（Feed_price）、收入（Revenue）、出售時間（Time_of_scale）為流量變量。其中收入（Revenue）表示流入，它會帶來總利潤的積累和增加，其余4個變量表示流出，會導致利潤的減少。

2.3 模型方程分析

根據系統動力學流圖中的參數關系，本研究構造了32個系統動力學方程，主要模型方程式如下：

d/dt（Fish）=Growth-Mortality;

d/dt（Feed）=-Growth+Feeding-Loss;

d/dt（Detritus）=+Accum-Decomp;

d/dt（Total_profit t）=+Profit;

d/dt（Numbers）=-J3;

Growth=if Feed>0 then C_growth*Feed*Fish_W/（Feed+C_Hs） else 0;

Feeding=if（Time0） then C_feed else 0;

Loss=C_loss*Feed+Growth*0;

Accum=Loss+Mortality;

Decomp=C_decomp*Detritus;

Profit=Revenue-Cost;

Mortality=if（TIME>Time_of_sale+1）then Fish/DT else（C_mort+Detritus^4）/（C_mort_d^4+Detritus^4））*Fish;

C_feed=A*（TIME+B）^2+C;

Fish_Price=10+2*Weight;

Revenue=if Time>Time_of_sale AND Time>Time_of_sale+2 then Fish_Price*Fish else 0;

Cost=Feed_price*Feeding。

3 結果與分析

3.1 養殖塘系統模型的建立

首先選擇一些喂養策略，使模型產生一些合理的結果。當A=0.001、B=-10、C=0.1或0.2時，得到圖2、圖3中的喂養方案和魚及碎屑量的動態變化。由圖2、圖3可以看出，魚類種群逐漸生長，直到某一點即第74天時，由餌料產生的碎屑量超過一定閾值時，導致魚類的大量死亡。魚類種群的崩潰，從而進一步增加了碎屑量。當碎屑分解率增加，碎屑量達到種群崩潰時的閾值也相應提高，則種群崩潰的時間也相應延長。說明系統模型模擬結果合理可靠，可以用來模擬分析養殖塘水生態系統。

3.2 養殖塘系統模型的優化

魚塘生態系統動力學模型建好后，再以總利潤最大為目標函數，用Berkeley Madonna軟件優化模型中的喂養策略因子A、B、C以及捕撈時間 Time_of_sale，優化界面見圖4。

用這個模型進行優化可以得到最大總收益為550萬元，如圖5所示。此時A=0.000 997 764，B=-13.517 6，C=0.171 658，Time_of_sale=75.401 6（圖6）。捕撈的時間定在碎屑即將達到閾值并使魚群處于滅絕的危險時。個別個體質量的增加與越來越多的魚死亡相抵消，魚的生物量迅速下降。因此選在第76天或者提前幾天將養殖塘中的魚進行捕撈后出售， 則能保證水產養殖戶的最大經濟效益。

3.3 模型的敏感性分析

參數靈敏度分析的目的是確定對模型預測結果影響較大的參數，找到這些參數后，預示著模擬預測時必須認真選取這些參數，如取值不當，會導致模擬失敗[21-22]。本研究利用綜合偏導數分別對生長率（C_growth）、死亡率（C_mortality）、餌料殘留率（C_loss）和碎屑分解率（C_decomp）參數進行敏感性分析，結果顯示其敏感度系數分別為0.46、-0.32、-0.35、0.58，其中絕對值越大，參數越靈敏。因此，碎屑分解率最為靈敏，其他依次為生長率、餌料殘留率和死亡率。檢驗的方式則是保持其他參數不變，只改變被檢驗參數，以最靈敏的碎屑分解率為例，以幅度為±20%、±50%、±100%進行擾動分析，觀察總利潤Total_profit的變化。由圖7可以看出，隨著C_decomp的逐漸增加，即池塘新陳代謝加快，池塘碎屑殘留量明顯減少，即水環境得到改善，因而魚的生物量增加，從而提高了養殖效益。

假如當其他參數未變，碎屑分解率從C_decomp=0.1增加到C_decomp=0.2時，即魚塘碎屑分解加快，則優化結果就會發生很大變化。在碎屑分解率較高的情況下，碎屑積累得很慢，很難達到使魚類死亡的臨界值。因此優化只能試圖使魚類的體質量盡可能增加，而消耗最少的飼料，具體追蹤優化過程如圖8所示。而此時優化結果則與之前發生了很大變化，可以看出得到了更高的利潤868（圖9），而最佳捕撈時間則從放養的第75天變為現在的第83天。因此，喂食情景對模型中的碎屑分解率（C_decomp 參數）非常敏感。在實際生產中可以通過池塘底部微孔曝氣方式增加溶解氧和添加微生態制劑等方式[23-24]來提高碎屑分解率，從而提高單位面積產量和經濟效益。

4 結論

本研究演示了系統動力學方法在魚塘生態養殖系統中的應用，可以看出，系統動力學方法直觀性強，在宏觀上將魚塘養殖過程中的生態、經濟、社會3個方面整合在一起，形成一個大的系統，融會貫通，且操作簡便，使用者即便沒有相當好的數學基礎，也可以根據不同的使用目的改變模型，根據實際情況改變相應參數設置。另外，本研究是將系統動力學模型應用于淡水漁業水產養殖領域的一項探索性研究，模型中大多采用簡化處理，部分參數采用估計值，而實際養殖系統中遠遠比本研究中建立的模型要復雜得多，尤其是參數的變動會在比較大的程度上影響結果，所以開展各項基礎調查研究以及如何完善模型以提高其仿真的準確性將是今后工作的重點和難點。

參考文獻：

[1]袁曉初，趙文武. 中國漁業統計年鑒（2015）[M]. 北京：中國農業出版社，2015.

[2]唐啟升，丁曉明，劉世祿，等. 我國水產養殖業綠色、可持續發展保障措施與政策建議[J]. 中國漁業經濟，2014，32（2）：5-11.

[3]宋 超，孟順龍，范立民，等. 中國淡水池塘養殖面臨的環境問題及對策[J]. 中國農學通報，2012，28（26）：89-92.

[4]鮑旭騰，徐 皓，張建華，等. 水產養殖面源污染控制的最佳管理實踐[J]. 南方水產科學，2012，8（3）：79-86.

[5]賈敬德. 淡水漁業環境現狀及保護對策[J]. 淡水漁業，2004，34（5）：59-61.

[6]賈敬德. 保護環境是漁業持續發展的當務之急[J]. 淡水漁業，1997，27（4）：26-29.

[7]Kenneth A，Bert B，Robert C，et al. Economic growth，carrying capacity，and the environment[J]. Science，1995，268（5210）：520-521.

[8]Jiang W M，Gibbs M T. Predicting the carrying capacity of bivalve shellfish culture using a steady，linear food web model[J]. Aquaculture，2005，244（1/2/3/4）：171-185.

[9]Carver C E A，Mallet A L. Assessing the carrying capacity of a coastal inlet in terms of mussel culture[J]. Aquaculture，1988，88：39-53.

[10]董雙林，李德尚，潘克厚. 論海水養殖的養殖容量[J]. 青島海洋大學學報，1998，28（2）：253-258.

[11]楊紅生，張福綏. 淺海筏式養殖系統貝類養殖容量研究進展[J]. 水產學報，1999，23（1）：84-90.

[12]劉劍昭，李德尚，董雙林. 關于水產養殖容量的研究[J]. 海洋科學，2000，24（9）：33-35.

[13]王其藩. 系統動力學[M]. 北京：清華大學出版社，1988：25-27.

[14]杰拉爾德·溫伯格. 系統化思維導論[M]. 張 佐，萬起光，董 菁，譯. 北京：清華大學出版社，2003：50-52.

[15]Ouyang Y，Zhang J E，Lin D，et al. A STELLA model for the estimation of atrazine runoff，leaching，adsorption，and degradation from an agricultural land[J]. Journal of Soils and Sediments，2010，10（2）：263-271.

[16]狄乾斌，徐東升，周樂萍. 基于STELLA軟件的海洋經濟可持續發展系統動力學模型研究[J]. 海洋開發與管理，2012，29（3）：90-94.

[17]于 濤，鐘 非，賀 鋒，等. 基于STELLA的循環水養殖系統池塘總氨氮動態模擬[J]. 農業機械學報，2013，44（7）：199-203.

[18]馮東溥，魏曉妹，降亞楠，等. 基于STELLA和氣候變化情景的灌區農業供需水量模擬[J]. 農業工程學報，2015，31（6）：122-128.

[19]黃歲樑，臧常娟，杜勝藍，等. pH、溶解氧、葉綠素a之間相關性研究Ⅰ：養殖水體[J]. 環境工程學報，2011，5（6）：1201-1208.

[20]孫 耀，陳聚法. 中國對蝦養殖水體中溶解氧的動態收支平衡模式[J]. 水產學報，1999，23（4）：424-428.

[21]邱康俊，溫華洋，王 根，等. 基于高斯隨機擾動的潛在蒸散敏感性分析[J]. 中國農學通報，2015，31（26）：156-160.

[22]何 亮，侯英雨，趙 剛，等. 基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數優化[J]. 農業工程學報，2016，32（2）：169-179.

[23]顧海濤，劉興國，何雅萍，等. 微孔曝氣式增氧機的性能及應用效果[J]. 漁業現代化，2017，44（3）：25-28.

[24]李 敏，段登選，許國晶，等. 大薸-微生態制劑協同凈化養殖池塘富營養化水體的效果[J]. 生態與農村環境學報，2015，31（1）：94-99.王宏玉，裴雪瑩，王 濤，等. 雜交黃顙魚黃優1號體質量和內臟指標的相關分析[J]. 江蘇農業科學，2020，48（2）：170-176.