糧堆磷化氫濃度衰減階段數學模型及驗證研究

郭 超 王殿軒 勞傳忠 沙寬寬 溫湛清

(廣東省糧食科學研究所糧食儲藏與害蟲防治研究室1，廣州 510310) (河南工業大學糧油食品學院2，鄭州 450001) (糧食儲藏與安全教育部工程研究中心3，鄭州 450001) (廣州市白云區糧食儲備公司4，廣州 510440) (珠海市三灶糧食收儲公司5，珠海 519040) (糧食儲運國家工程實驗室6，廣州/鄭州 510310/450001)

磷化氫是當前世界范圍內應用最普遍的熏蒸劑，廣泛應用在糧食、煙草、干果、藥材等儲藏物的殺蟲處理。在中國的儲糧熏蒸主要應用則以磷化氫倉外發生器施藥、磷化鋁糧面潮解施藥、磷化鋁布袋潮解施藥、磷化鋁探管施藥等方式下的環流或非環流熏蒸技術應用為主。同磷化氫熏蒸殺蟲技術，因技術掌握程度、應用操作水平、倉房密閉條件等差異，而出現不同的殺蟲效果。磷化氫熏蒸失敗的事例時有所見，二次或多次熏蒸的情況仍有存在[1]。磷化氫熏蒸殺蟲的關鍵主要取決于有效濃度和在有效濃度下的熏蒸時間，實現熏蒸成功必須確保這兩個影響因素超過一定的閾值[2-3，7，9]。熏蒸時間一般為糧堆密閉時間，在實際熏蒸過程易于掌握，而熏蒸中采用和保持合理的磷化氫濃度至關重要[4-5]。分析掌握熏蒸過程中磷化氫的濃度變化規律，對科學使用磷化氫并有效殺滅害蟲具有重要的現實意義。

熏蒸是在密閉環境中進行的，一般施藥后熏蒸劑濃度變化的基本過程分3個階段，即毒氣發生階段、毒氣衰減階段、毒氣排放階段[1]。國內外對這3個階段磷化氫濃度的變化規律報道主要分為基于熏蒸劑對流擴散和吸附特性建立的糧堆內部分布模型，以及糧堆熏蒸劑平均濃度隨時間的變化模型。關于熏蒸劑在糧堆內部分布模型的類似研究也有一些探討，王遠成等[6]建立和驗證了熏蒸劑(磷化氫)的對流擴散和吸附模型，并采用計算流體動力學方法對圓筒倉內谷物熏蒸過程中磷化氫質量濃度進行了數值預測，分析了熏蒸過程中磷化氫的質量濃度分布規律；Isa等[7]采用計算流體動力學方法，對立筒倉強制環流熏蒸和糧面施藥熏蒸中磷化氫的濃度分布進行預測分析，研究發現立筒倉泄漏的部位對環流熏蒸中熏蒸劑的分布影響很大。若泄漏部位在立筒倉的底部，表層的害蟲將殺滅困難，底部泄漏部位不影響糧面施藥后熏蒸劑的分布，但接近立筒倉底層的磷化氫濃度保持低。采用計算流體動力學方法建立的糧堆內部熏蒸劑分布模型在改進熏蒸工藝設計時有一定參考，由于在建模過程多采用微積分，熏蒸實踐中難以推廣應用。張立力等[8]利用Wolfram對大型立筒倉磷化氫外環流熏蒸過程濃度和時間的變化建立模型，其中熏蒸劑發生階段符合二次函數模型，濃度衰減階段和散氣階段均符合冥函數模型；Cryer等[9]分別報道了硫酰氟熏蒸過程中熏蒸劑濃度衰減階段和散氣階段的動力學模型符合C/C0=exp(-kt)(C0，熏蒸劑初始濃度；C，t時的熏蒸劑濃度；k,單位時間內衰減速度常數)；中國儲備糧管理總公司等[1]報道了在不補充熏蒸劑到密閉環境中，熏蒸的濃度和時間呈指數模型變化C=C0exp(-(t-t0)k)或lnC0-lnC=k(t-t0)(C0，t0時刻熏蒸劑濃度；C，t時的熏蒸劑濃度；k,單位時間內衰減速度常數)，其中熏蒸密閉期間熏蒸劑的濃度在初始階段迅速降低，之后以半對數的關系減少。從已報到的研究來看，其中提到的模型種類較多，但一些模型中含有的熏蒸劑初始濃度難以界定，磷化氫初始濃度的精確度影響到模型準確性，由此導致此類模型指導實際工作偏離度較大。

作者研究了在不進行熏蒸劑補充的密閉環境條件下，糧堆內磷化氫濃度隨時間的變化規律，建立了糧堆內磷化氫平均濃度和熏蒸時間的模型，并進行了實倉驗證，以期為解決生產中抗性害蟲防治提供指導。

1 材料與方法

1.1 實驗倉房

建模實驗倉房采用廣州市白云區糧食儲備公司九佛直屬糧庫04-1A01倉，模型驗證倉房采用珠海市三灶糧食收儲公司P1倉、P2倉、P3倉。實驗倉房均為高大平房倉，其儲糧熏蒸基本條件見表1所示。

1.2 實驗方法

1.2.1 倉房氣密性測定

參照GB/T 25229—2010 平房倉氣密性要求。

1.2.2 磷化氫濃度測定

廣州市白云區糧食儲備公司九佛直屬糧庫04-1A01倉和珠海市三灶糧食收儲公司P1倉、P2倉、P3倉均在糧堆表層四角及中央設置5個氣體取樣點，用橡膠管引至倉外磷化氫氣體檢測箱內。糧堆的磷化氫濃度為5個氣體取樣點的平均濃度。熏蒸過程中，九佛直屬糧庫04-1A01倉不實施補藥，珠海市三灶糧食收儲公司P1倉、P2倉、P3倉濃度低于200 mL/m3時立即補藥。

1.2.3 數據處理

模型的擬合和參數計算均采用SPSS 17.0和Excel 2010分析，本研究的分析結果為Excel 2010計算的結果，其中參數計算采用最小二乘法。數據的相關性分析采用SPSS 17.0雙變量相關性分析，相關系數r選取Pearson相關系數。

2 結果與分析

2.1 熏蒸過程中糧堆內磷化氫濃度隨時間變化

九佛糧庫04-1A01倉磷化氫濃度隨時間的變化如圖1所示。由圖1可知，在倉內糧面投放43.5 kg磷化鋁后，磷化氫濃度逐漸升高。經第2 天環流4 h和第8天環流2.5 h，倉內磷化氫濃度在第9 天達到最高值，為692.2 mL/m3。隨著時間的延長，磷化氫濃度逐漸降低，第28天濃度降低為194 mL/m3，第45 天濃度降低至33.6 mL/m3。對磷化氫發生階段和衰減階段濃度C和熏蒸時間t擬合曲線，結果表明磷化氫發生階段濃度C和熏蒸時間t遵循C=2.62t3-44.50t2+264.87t模型(R2=0.960 24)，衰減階段遵循C=1 649.1e-0.080 89t模型(R2=0.977 8)。根據衰減階段模型計算出的磷化氫濃度與實測值相關系數0.990 8，這表明預測值與實測值具有較高的相關性。

圖1 九佛糧庫04-1A01倉磷化氫濃度隨時間的變化

2.2 磷化氫濃度-時間模型的建立及應用

基于九佛糧庫04-1A01倉磷化氫發生階段濃度C和熏蒸時間t遵循C=2.62t3-44.50t2+264.87t模型(R2=0.960 24)，衰減階段遵循C=1 649.1e-0.080 89t模型(R2=0.977 8)，假定熏蒸過程中磷化氫濃度C和熏蒸時間t遵循上述模型，以磷化氫濃度最大值為分界點，建立磷化氫發生階段和衰減階段磷化氫濃度-時間模型，如下：

磷化氫發生階段：C=αt3+βt2+γt(α>0，β<0,γ>0) Ⅰ

磷化氫衰減階段：C=ae-bt(a>0；b>0) Ⅱ

其中C為磷化氫濃度，mL/m3；t為熏蒸時間,d；e為自然指數；α,β,γ,a,b均為常數。

2.2.1 模型Ⅱ擬合所需的數據量分析

分別采用3個點至8個點擬合模型Ⅱ，采用最小二乘法計算獲得a、b，并計算擬合度R2，如表2所示。從表2可知，采用3個點至8個點擬合模型Ⅱ，其數據與模型的擬合度均大于0.97，這表明采用3個或3個以上的數據即可得到擬合度高的模型Ⅱ。同時，采用不同數量的數據所得到的模型，預測磷化氫衰減階段濃度的預測值與實測值相關性r均大于0.98，這表明采用不同數量的數據所得到的模型Ⅱ均可以較好的預測。如九佛糧庫04-1A01倉磷化氫在28 d時濃度為194.6 mL/m3，而采用不同數量的數據所得到的模型Ⅱ分別預測濃度降至194 mL/m3的熏蒸天數如表2所示，結果表明采用3個數據擬合的模型預測的熏蒸天數誤差最大，超過實際天數的5.6 d，但隨著數據量的增多，預測的熏蒸天數由33.6 d逐漸減少并趨近于28 d，這表明隨著數據量的增加，模型Ⅱ的預測的準確性不斷提高。

表2 不同數據量擬合的模型參數變化及預測值和實測值的相關性

2.2.2 計算磷化氫熏蒸最大濃度

根據磷化氫發生階段模型I：C=αt3+βt2+γt(α>0,β<0,γ>0)和衰減階段模型Ⅱ：C=a×e-b×t(a>0;b>0)，在t>0 d時模型分別為增函數和減函數。令αt3+βt2+γt=a×e-b×t，則兩模型在t時刻相等，此刻磷化氫濃度為最大，即當t=9.6 d時，磷化氫濃度最大值Cmax=758.9 mL/m3[10]。磷化氫濃度半數衰減時間HLT為8.6 d。

2.2.3 計算磷化氫衰減階段濃度和熏蒸時間

由C=ae-bt(即C=1 649.1e-0.080 89t)，計算到不同熏蒸時間對應的磷化氫濃度，如表3所示。由表3可知，實驗期間九佛糧庫04-1A01倉磷化氫濃度隨時間的增加逐漸降低，且磷化氫濃度下降的幅度逐漸降低，其中14～21 d濃度降低了229.7 mL/m3，28～35 d濃度降低了74.1 mL/m3，42～49 d濃度降低了23.9 mL/m3，但每隔7 d濃度均以43.2%的比例降低。

表3 基于模型Ⅱ計算出的不同熏蒸時間對應的磷化氫濃度

令磷化氫濃度發生階段模型Ⅰ：C=2.62t3-44.50t2+264.87t和衰減階段模型Ⅱ：C=1 649.1e-0.080 89t濃度C分別為350、300、250、200、150、100 mL/m3，可計算到磷化氫濃度發生階段和衰減階段不同濃度下對應的熏蒸時間，如表4所示。

表4 基于模型Ⅰ和Ⅱ計算出的不同磷化氫濃度對應的有效熏蒸時間

2.3 模型Ⅱ在實倉熏蒸中的相關性與吻合度驗證

2.3.1 P1倉實際熏蒸濃度變化與模型驗證關系

P1倉壓力從500 Pa降低至250 Pa，壓力半衰期為105 s。熏蒸過程中P1倉共補藥2次。由圖2可知，磷化氫衰減階段濃度和時間變化分3個階段，均符合模型Ⅱ，且R2介于0.949 91～0.994 21，這表明模型的擬合度較高。由表5可知，P1倉磷化氫衰減階段階段Ⅰ、階段Ⅱ、階段Ⅲ濃度預測值與實測值的相關性r分別為0.995 6、0.988 5、0.970 4，這表明模型Ⅱ擬合效果好，磷化氫預測濃度接近實測值。

圖2 三灶糧庫P1倉磷化氫濃度隨時間的變化

圖3 三灶糧庫P2倉磷化氫濃度隨時間的變化 表5 三灶糧庫P1倉磷化氫濃度預測值與實測值的相關性

模型實測值/mL/m3預測值/mL/m3r階段Ⅰ:C=1139.3e-0.17001t582.0,517.0,443.5,304.8,253.3602.2,508.1,428.6,305.1,257.40.9956階段Ⅱ:C=1895e-0.11637t468.6,427.8,403.8,305.0,263.0482.8,429.8,382.5,303.1,269.80.9885階段Ⅲ:C=2444.5e-0.08712t442.6,424.5,412.0,380.3,363.0,297.5,201.0477.3,437.4,401.0,367.5,336.8,283.0,217.90.9704

表6 三灶糧庫P2倉磷化氫濃度預測值與實測值的相關性

表7 三灶糧庫P2倉階段I磷化氫實測濃度對應的熏蒸時間

表8 三灶糧庫P2倉階段II磷化氫實測濃度對應的熏蒸時間

表9 三灶糧庫P3倉磷化氫濃度預測值與實測值的相關性

2.3.2 P2倉實際熏蒸濃度變化與模型驗證關系

P2倉壓力從500 Pa降低至250 Pa，壓力半衰期為96 s。熏蒸過程中P1倉共補藥1次。由圖3可知，磷化氫衰減階段濃度和時間變化分兩個階段，均符合模型Ⅱ，且R2介于0.950 81～0.962 28，這表明模型的擬合度較高。由表6可知，P2倉磷化氫衰減階段階段Ⅰ、階段Ⅱ濃度預測值與實測值的相關性r分別為0.958 3、0.984 9，這表明模型Ⅱ擬合效果好，磷化氫預測濃度接近實測值。此外，由表7、表8可知，P2倉磷化氫衰減階段Ⅰ和階段Ⅱ熏蒸時間預測值與實測值的相關性r分別為0.988 6、0.996 2，這表明模型II擬合效果好，磷化氫預測的熏蒸時間接近實測值。

2.3.3 P3倉實際熏蒸濃度變化與模型驗證關系

P3倉壓力從500 Pa降低至250 Pa，壓力半衰期為69 s。熏蒸過程中P3倉共補藥2次。由圖4可知，磷化氫衰減階段濃度和時間變化均符合模型Ⅱ，且R2介于0.924 54～0.999 38，這表明模型的擬合度較高。由表9可知，P3倉磷化氫衰減階段階段Ⅰ、階段Ⅱ、階段Ⅲ濃度預測值與實測值的相關性r分別為0.975 0、0.996 3，0.959 2，這表明模型Ⅱ擬合效果好，磷化氫預測濃度接近實測值。此外，由表10可知，P3倉不同磷化氫濃度下的預測熏蒸時間與圖4吻合，其中第9.5 天、第18.2天時的磷化氫濃度為200 mL/m3。2015年12月8日P3倉開始施藥熏蒸，根據熏蒸方案，此時應對糧堆補藥。實際生產中，P3倉分別于第9 天(12月17日)和第19 天(12月27日)施藥，這表明模型II擬合效果好，不同濃度下的磷化氫預測的熏蒸時間接近實測值。當磷化氫濃度下降至200 mL/m3時，階段Ⅲ模型預測到的熏蒸時間為30 d，即下次補藥應該在2016年1月7日。

圖4 三灶糧庫P3倉磷化氫濃度隨時間的變化 表10 三灶糧庫P3倉不同磷化氫濃度下的預測熏蒸時間

磷化氫濃度/mL/m3熏蒸時間預測值/d階段Ⅰ:C=1316.5e-0.19821t階段Ⅱ:C=3358e-0.15458t階段Ⅲ:C=3000.7e-0.09021t3506.714.623.83007.515.625.52508.416.827.52009.518.230.015011.020.133.210013.022.737.7

3 討論

本實驗以廣州市白云區糧食儲備公司九佛直屬糧庫04-1A01倉為建模實驗倉，在不進行熏蒸劑補充的密閉環境條件下，分析了糧堆內磷化氫濃度隨時間變化規律，構建了糧堆內磷化氫平均濃度和熏蒸時間的關系模型，其中磷化氫發生階段模型Ⅰ：C=2.62t3-44.50t2+264.87t和衰減階段模型Ⅱ：C=1 649.1e-0.080 89t。模型II在珠海市三灶糧食收儲公司P1倉、P2倉、P3倉進行實倉驗證。結果表明，衰減階段磷化氫濃度隨熏蒸時間的增加，磷化氫濃度呈指數下降趨勢。衰減階段模型II能夠用于計算不同磷化氫濃度對應的熏蒸時間、不熏蒸時間對應的磷化氫濃度等指標，結合模型I還可計算磷化氫熏蒸最大濃度、最低熏蒸濃度下的有效熏蒸時間以及磷化氫濃度半數衰減時間HLT。此外，模型II還可用于預測糧堆的補藥時間，指導糧食倉儲行業實施磷化氫熏蒸作業。

九佛糧庫04-1A01倉平均糧溫為23.5 ℃，最低點糧溫為22.3 ℃。根據GB/T 29890—2013 《糧油儲藏技術規范表》E.1，當溫度為20～25℃時不同蟲種不同密閉時間的磷化氫熏蒸最低有效濃度設定要求，磷化氫保持最低有效濃度為350 mL/m3的密閉時間保持超過14 d以上，300 mL/m3時應超過21 d，250 mL/m3時應超過28 d。由表3可知，磷化氫保持最低有效濃度分別為350、300、250 mL/m3時，糧堆密閉的時間分別為17.35、19.58、22.37 d。雖磷化氫濃度為300、250 mL/m3時糧堆密閉時間均低于GB/T 29890—2013 《糧油儲藏技術規范表》E.1規定，但磷化氫濃度為300 mL/m3時密閉時間超過14 d，從理論上講，該糧堆中扁谷盜類(屬)、蛾類、谷蠹、米象、書虱、螨類、赤擬谷盜、米扁蟲及其他抗性蟲種均應被完全殺滅。

三灶糧庫P1倉、P2倉、P3倉倉房壓力半衰期分別為105、96、69 s，而三個倉階段I單位時間內衰減速度常數b分別為0.170 01、0.176 21、0.198 21。從數值上看，隨著壓力半衰期的減少，衰減常數逐漸增大，即倉房氣密性越差，倉內磷化氫衰減速度越快，這與實際情況相一致。同樣，以三灶糧庫P1倉階段Ⅰ、階段Ⅱ、階段Ⅲ為例，單位時間內磷化氫衰減速度常數b分別為0.170 01、0.116 30、0.087 12。從數值上看，P1倉磷化氫衰減速度逐漸減慢。隨著倉內磷化氫濃度的降低，倉內外濃度差逐漸減小，這可能導致磷化氫衰減速度逐漸減慢。在氣密性較好的倉內熏蒸，對大多數熏蒸劑來說吸著將成為影響衰減速率常數的主要部分[12],因此，磷化氫衰減速度逐漸減慢還可能與糧堆吸附磷化氫有關。此外，九佛糧庫04-1A01倉采用負壓法測試倉房氣密性，而三灶糧庫P1倉、P2倉、P3倉采用正壓法測試倉房氣密性。高彬彬等[13]報道的北11倉和南11倉、周長金等[14]報道的安陸45倉和安陸49倉磷化氫熏蒸過程衰減階段的濃度變化均符合本研究建立的模型。以上倉房無論使用負壓法或正壓法，熏蒸過程磷化氫衰減階段濃度變化符合本研究提到的模型，實際上由于倉房壓力半衰期測定主要在新建倉房驗收時使用，在糧食保管過程中測定較少，而基于單位時間內衰減速度常數b可能與倉房壓力半衰期、糧堆吸附磷化氫有關，進一步研究單位時間內衰減速度常數b，可能更好用于指導磷化氫熏蒸。

4 結論

糧堆內磷化氫濃度衰減階段的濃度呈指數下降趨勢。磷化氫濃度衰減階段濃度和時間符合指數模型C=ae-bt(a>0；b>0)(其中C為磷化氫濃度,mL/m3；t為熏蒸時間,d；e為自然指數)。該模型可計算不同磷化氫濃度對應的熏蒸時間、不同熏蒸時間對應的磷化氫濃度、預測糧堆的補藥時間等，用于指導糧食倉儲行業磷化氫熏蒸作業。

致謝：感謝廣州市白云區糧食儲備公司九佛直屬糧庫戴惠良、劉柄和、滕俊濤等和珠海市三灶糧食收儲公司伍燦洪、曹文杰、王山、黃海濤、肖小波、蘇雄彬、陳錦銳等的工作。感謝廣東省糧食科學研究所曾伶、盧木波、張少波、楊永強、冼慶、鄭妙、張小松等對本工作的支持。

[1]中國儲備糧管理總公司,河南工業大學.儲糧磷化氫熏蒸技術區域優化[M].北京:中國農業科學技術出版社,2008 China National Grain Reserves Corporation, Henan University of Technology.Regional optimization of phosphate fumigation technology for sorage grain[M]. Beijing:China Agricultural Science and Technology,2008

[2]CHAYAPRASERT W, NUKHAM K, SUKCHAROEN A. Evaluation of the superposition method for predicting gas leakage rates during fumigations in empty model silos[J].Journal of Stored Products Research, 2015, 64, Part B: 13-20

[3]CHAYAPRASERT W, MAIER D E, SUBRAMANYAM B, et al. Gas leakage and distribution characteristics of methyl bromide and sulfuryl fluoride during fumigations in a pilot flour mill[J].Journal of Stored Products Research, 2012, 50: 1-7[4]王殿軒,卞科.儲糧熏蒸劑的發展動態與前景[J].糧食儲藏, 2004, 32(4): 3-7 WANG D X,BIAN K. Development trend of fumigant used in stored grain[J]. Grain Storage, 2004, 32(4): 3-7

[5]白旭光.儲藏物害蟲與防治[M].北京:科學出版社,2008 BAI X G. Stored product pest and control[M].Beijing: Science Press,2008

[6]王遠成, Thorpe G R, 趙會義,等.儲糧熏蒸過程中磷化氫濃度的分布模型及驗證研究[J].中國糧油學報, 2015, 30(7): 81-84 WANG Y C, THORPE G R,ZHAO H Y,et al. The model and validation of phosphine fumigation in stored grains in silos[J]. Journal of the Chinese Cereals and Oils Association, 2015, 30(7): 81-84

[7]ISA Z M, FARRELL T W, FULFORD G R, et al.Mathematical modelling and numerical simulation of phosphine flow during grain fumigation in leaky cylindrical silos[J].Journal of Stored Products Research, 2016, 67: 28-40

[8]張立力,霍洪源,趙志永,等.立筒庫磷化氫環流熏蒸模型的研究[J].糧食儲藏,2000, 28(1): 27-32 ZHANG Z L,HUO H Y,ZHAO Z Y,et al.Study on silo recirculation fumigation model of phosphine in china[J]. Grain Storage, 2000, 28(1): 27-32

[9]CRYER SA,BARNEKOW DE.Estimating Outside Air Concentrations surrounding Fumigated Grain Mills[J]. Biosystems Engineering, 2006, 94(4): 557-572

[10]Wolfram Alpha LLC.Wolframalpha computational knowledge engine[EB/OL].[2016-03-02].http://www.wolframalpha.com/

[11]REED R C. Managing Stored Grain to Preserve Quality and Value[M].AACC International, 2005

[12]王殿軒,曹陽.磷化氫熏蒸殺蟲技術[M].成都：成都科技大學出版社，1999 WANG D X,CAO Y. Phosphine fumigation technology[M].Chengdu: Chengdu University of Science and Technology Press,1999

[13]高彬彬,沈天翔,金宗銓,等.平房倉磷化氫膜下環流熏蒸技術的應用研究[J].糧食儲藏, 2000,6:41-48 GAO B B,SHEN T X,JIN Z Q,et al.Trials of phosphine recirculation fumigation in grain bulk sealed surface with plastic sheet in horizontal bins[J]. Grain Storage, 2000,6:41-48

[14]周長金,馬明君,李錫偉,等.不同氣密性對環流熏蒸效果的影響[J].糧食儲藏, 2003,5:21-23 ZHOU C J,MA M J,LI X W,et al. Effect of different gastightness on recirculation fumigation[J]. Grain Storage, 2003,5:21-23.