王鵬杰 祁智慧 張海洋 田 琳 高瑀瓏 唐 芳

(國家糧食和物資儲備局科學研究院1，北京 100037)

(南京財經(jīng)大學2，南京 210023)

稻谷儲藏期間，霉變是造成儲糧損失的一個重要因素。近幾十年來，國內(nèi)外有關稻谷儲藏霉變的研究報道很多。Laca等[1]研究發(fā)現(xiàn)霉菌主要分布在稻谷表面。Genkawa等[2]研究了不同水分的稻谷儲藏過程中霉菌數(shù)量的變化，低水分稻谷無霉菌生長，其儲藏的效果同低溫儲藏一樣有效。周建新等[3-5]研究了不同貯藏條件下稻谷真菌的變化情況，環(huán)境溫度和稻谷含水量對稻谷儲糧真菌生長有直接影響，儲藏溫度和時間與霉菌量呈顯著的二元線性關系。金昌福[6]等利用近紅外檢測技術，建立了貯藏環(huán)境溫度和水分含量對貯藏稻谷表面霉菌菌落總數(shù)影響規(guī)律的多元線性回歸模型。唐芳等[7]通過對稻谷儲藏過程中水分和溫度的變化與真菌生長關系的研究，得出稻谷儲存水分、溫度與真菌起始生長時間的關系曲線。上述研究都是側重于一個或兩個因素的影響，而真菌生長是受多個因素的影響，且相互耦合共同對真菌生長產(chǎn)生影響，因此，要更真實地反映真菌生長狀況，需基于大量數(shù)據(jù)進行統(tǒng)計分析，建立多影響因素的真菌生長預測模型。

多元線性回歸是多元統(tǒng)計分析中一種重要方法，可以解決多個因素對同一變量影響的相關性分析，進而建立相關性的統(tǒng)計模型。Green等[8]基于室內(nèi)溫度、濕度等變量建立了室內(nèi)有害微生物的多元線性回歸預測模型，預測模型精度可達97%。南英華等[9]建立了泉流量與大氣降水的多元線性回歸模型。本文基于水分、溫度和儲藏時間三個主要影響因素，定期采集真菌生長數(shù)據(jù)，采用多元線性回歸建立儲糧真菌生長預測模型，并對模型進行了F檢驗、t檢驗和殘差分析，結合稻谷實倉數(shù)據(jù)對預測模型進行了初步驗證。

1 材料與方法

1.1 模擬儲藏實驗材料及樣品處理

取自黑龍江的粳稻樣品，采用噴霧加無菌水方式，將樣品水分含量調(diào)至目標水分，密封于4 ℃低溫均衡水分。目標水分含量較高的樣品，需多次加水，低溫均衡水分時間至少30 d，直至稻谷含水量含量達到目標水分并通過水分均勻性檢驗。

1.2 試劑和儀器

HPS-250生化培養(yǎng)箱；PL3002-IC電子分析天平；SMART顯微鏡。

1.3 實驗方法

1.3.1 模擬儲藏及檢測周期

將不同水分含量(13.0%、13.5%、14.0%、14.5%、15.0%、15.5%、16.0%、16.5%、17.0%、17.9% )的稻谷樣品，密封并分別置于不同溫度的(10、15、20、25、30、35 ℃)生化培養(yǎng)箱中模擬儲藏180 d，每10 d取樣一次，檢測儲糧真菌生長數(shù)量。

1.3.2 稻谷實倉情況及取樣方法

選擇華北地區(qū)某糧庫稻谷倉為實驗倉。稻谷倉基本信息：2017年1月入庫，產(chǎn)地黑龍江，容量6 000 t，入倉水分為14.5%。取樣位置的選擇與糧情檢溫電纜表層布點位置相重合，具體分布如圖1所示。取樣深度為距離糧面0.4 m深度。度夏期間采用手動扦樣器定期取樣，每次取樣200 g，置于無菌袋中保存，取回后4 ℃低溫保存，以備檢測。

圖1 實驗倉扦樣點位置示意圖

1.3.3 儲糧真菌檢測方法

參照LS/T 6132《糧食檢驗 儲糧真菌的檢測 孢子計數(shù)法》。

1.3.4 稻谷含水量測定方法

參照GB/T 5497 糧食、油料檢驗水分測定法 105 ℃烘干法。

1.3.5 糧溫檢測方法

參照糧庫內(nèi)糧情檢測系統(tǒng)表層糧溫數(shù)據(jù)。

1.4 數(shù)據(jù)處理方法

1.4.1 真菌生長數(shù)據(jù)預處理方法

本實驗設置了10個水分梯度的稻谷樣品，6個模擬儲藏溫度，180 d中，對每個梯度的樣品每10 d取一次樣，檢測儲糧真菌生長數(shù)量，共取19次，總計1 140組數(shù)據(jù)，每次取樣做雙實驗，兩次結果偏差超過50%時，重新取樣檢測，檢測結果取平均值。

真菌孢子數(shù)在103～107個/g范圍，且微生物生長為非線性，根據(jù)微生物生長模型處理經(jīng)驗[10]，對真菌孢子數(shù)進行對數(shù)變換，不僅利于線性回歸，而且可以縮小數(shù)據(jù)的絕對值以方便數(shù)據(jù)處理，降低樣本的異方差程度[11]。真菌檢測方法的檢出限為3×104個/g，對于未檢出真菌生長的樣品，并不代表糧食不攜帶真菌孢子，而是孢子濃度未達到該方法檢出限，依據(jù)環(huán)境及正常糧食真菌帶菌量經(jīng)驗數(shù)據(jù)[12]，將該方法未檢出的稻谷樣品帶菌量補為1×103個/g，更符合稻谷帶菌量的實際情況。

1.4.2 數(shù)據(jù)分析軟件

文中數(shù)據(jù)結果由MATLAB2016軟件處理獲得。

2 結果與討論

糧食儲糧真菌的生長受多個因素的影響。多元線性回歸模型，可以實現(xiàn)多個自變量的最優(yōu)組合共同來預測或估計因變量。其主要步驟如下：第一步，確定自變量和因變量；第二步，根據(jù)研究現(xiàn)象，依據(jù)相應的理論和經(jīng)驗，設定模型并加以確定；第三步，參數(shù)估計；第四步，模型的檢驗和修正，常用的檢驗有擬合優(yōu)度檢驗、線性回歸模型F檢驗、參數(shù)的t檢驗以及殘差分析。第五步，模型的運用[13]。下面重點介紹建模、模型檢驗及應用。

2.1 多元線性回歸建模

多元線性回歸是一種數(shù)理統(tǒng)計方法，設因變量為y，自變量為x，自變量有p個，它們的n組觀測值為(x1i，x2i，xpi，yi)(i=1，2，…，n)，其多元線性回歸表達式如下所示：

β是多元線性回歸方程的系數(shù)。

在數(shù)據(jù)分析過程中，將所得實驗數(shù)據(jù)導入MATLAB，分別繪制出取對數(shù)后的真菌孢子數(shù)與環(huán)境溫度，稻谷含水量和儲藏時間的散點圖，見圖2。

從圖2中所示的散點圖可看出，取對數(shù)后的真菌孢子數(shù)分別與環(huán)境溫度、稻谷含水量和儲藏時間大致呈線性關系。因此，可以利用MATLAB對其進行多元線性回歸。

本研究對上述實驗采集的1 140組實驗數(shù)據(jù)，將稻谷儲藏溫度、含水量、儲藏時間作為自變量，取對數(shù)處理的真菌孢子數(shù)作為因變量，用MATLAB軟件regress語句對這些實驗數(shù)據(jù)進行擬合，得到的多元線性回歸模型如下：

圖2 取對數(shù)后真菌孢子數(shù)與環(huán)境溫度、儲藏時間及稻谷含水量的散點圖

lny=-20.787+0.157×x1+0.027×x2+1.683×x3

(1)

式中：y表示真菌孢子數(shù)，單位為個/g；x1表示環(huán)境溫度，單位為℃；x2表示儲藏時間，單位為d；x3表示稻谷含水量，以%表示。

式(1)為以lny為因變量的多元線性方程，對方程(1)進行指數(shù)變換，得到以y為因變量的指數(shù)函數(shù)方程。

y=exp(-20.787+0.157×x1+0.027×x2+1.683 ×x3)

(2)

2.2 預測模型的效果分析

根據(jù)實驗數(shù)據(jù)初步建立了多元線性回歸模型，還須通過模型的擬合優(yōu)度、方程線性關系的顯著性、系數(shù)的顯著性等統(tǒng)計量的檢驗后，才可以用于解釋、分析實際問題[14]。

2.2.1 擬合優(yōu)度檢驗

擬合度是用于檢驗回歸方程對樣本觀測值的擬合程度。多元線性回歸的擬合程度，使用多重判定系數(shù)，其定義為：

式中: SSR為回歸平方和，SSE為殘差平方和，SST為總離差平方和。R2表示因變量與所有自變量之間的線性相關程度，實際反映的是樣本數(shù)據(jù)與預測數(shù)據(jù)間的相關程度，R2越接近于1，方程擬合優(yōu)度越高[15]。R2的平方根稱為復相關系數(shù)(R)，本模型的相關系數(shù)R= 0.878，判定系數(shù)R2= 0.770，這說明取對數(shù)后的儲糧真菌孢子數(shù)的變化有77%可由環(huán)境溫度、稻谷含水量和儲藏時間這三個因素的變動來解釋。

2.2.2 線性關系顯著性檢驗(F檢驗)

在建立多元線性回歸模型之后，還必須對因變量與多個自變量間的線性關系的假設進行顯著性檢驗。F統(tǒng)計量定義為：平均的回歸平方和與平均的殘差平方和之比，對于多元線性回歸方程：

式中: SSR為回歸平方和;SSE為殘差平方和;n為樣本數(shù);k為自變量個數(shù)。F統(tǒng)計量服從第一自由度為k、第二自由度為n-k-1的F分布，即F～F(k，n-k-1)。通過查詢F分布分位數(shù)表，我們可以得出特定顯著度條件下F檢驗的臨界值。當統(tǒng)計量F的值大于臨界值時，即可認為在總體上自變量與因變量呈顯著的線性關系。F值越大，線性回歸效果越顯著。

在本模型的F檢驗中，總平方和為15 150.32，回歸平方和為11 667.82，殘差平方和為3 482.50。本研究顯著性檢驗的顯著性概率值為0，顯然滿足P<α=0.05。F統(tǒng)計的觀測值為1 219.552，查表得F0.05(3，1 136)= 2.61，F(xiàn)值遠大于F0.05(3，1 136)，表明因變量與自變量之間整體線性關系顯著，即取對數(shù)后的真菌孢子數(shù)與環(huán)境溫度，稻谷含水量和儲藏時間存在顯著多元線性關系。

2.2.3 回歸系數(shù)顯著性檢驗(t檢驗)

在多元線性回歸方程中，回歸方程顯著性F檢驗和回歸系數(shù)t檢驗是不等價的，還需進行方程回歸系數(shù)的t檢驗，以判斷各相關因子之間是否存在共線性現(xiàn)象。統(tǒng)計量t定義為：

給定一個顯著性水平α，可以從t分布分位數(shù)表查得tα/2(n-k-1)。當|ti| >tα/2(n-k-1)，P<α=0.05，接受備擇假設H1。本模型各項變量的t顯著性概率均為0，均小于0.05，環(huán)境溫度、儲藏時間和稻谷含水量的t統(tǒng)計量分別為24.785、27.175、46.607，均大于t0.025(1 136)= 1.96。由此可見，稻谷含水量、儲藏時間和環(huán)境溫度對取對數(shù)后的儲糧真菌孢子數(shù)均有顯著影響。

2.2.4 殘差分析

殘差是實際觀察值與回歸估計值的差。從殘差圖可以看出數(shù)據(jù)的殘差離零點的遠近，當殘差的置信區(qū)間均包含零點，說明回歸模型能較好的符合原始數(shù)據(jù)，否則視為異常點。

對1 140組實驗數(shù)據(jù)進行殘差分析后的殘差圖如圖3所示，由此可以確定殘差落在其置信區(qū)間內(nèi)的大致位置，也可以觀察殘差的分布變化的趨勢，殘差圖越散亂代表模型的適配越好[16]。經(jīng)分析，異常點主要分布在35 ℃、13%水分下儲藏5個月后和35 ℃、13.5%水分下儲藏4個月后的數(shù)據(jù)點，以及高溫高水分儲藏條件下的儲藏初期數(shù)據(jù)點，共計44組數(shù)據(jù)。這些異常數(shù)據(jù)中，高溫低水分區(qū)(13.0%和13.5%)真菌孢子數(shù)量基本處于檢測方法的最低檢出限附近，檢測結果誤差偏大，高溫高水分區(qū)(17.9%)，儲糧真菌生長速度較快，10 d的取樣頻率，無法檢測到真菌孢子數(shù)量逐漸升高的趨勢，因此初期檢測數(shù)據(jù)波動較大，影響模型的準確率。經(jīng)過一次殘差分析，剔除異常點后，模型的擬合度R2從0.72上升到0.77。若經(jīng)過多次篩選完全剔除異常點后，擬合優(yōu)度能達到0.84，但此時的結果是為了達到數(shù)學上的最優(yōu)化而剔除了一部分符合實際的實驗數(shù)據(jù)，可能影響模型的真實性，因此選擇第一次剔除異常點后的數(shù)據(jù)點進行擬合建模。

2.3 實倉檢測數(shù)據(jù)驗證

為了便于對真菌生長預測模型在實際應用中預測效果進行評價，根據(jù)標準“LS/T 6132 糧食檢驗 儲糧真菌的檢測 孢子計數(shù)法”[17]附錄C“儲糧安全評價參考表”，將危害真菌孢子檢出數(shù)量分為4個級別，詳見表1。在實際儲藏中，儲糧樣品真菌孢子檢出數(shù)量基本可代表糧食樣品霉變程度，依據(jù)參考表中級別對儲糧安全狀況進行初步判定，進而指導相應處置措施。

表1 儲糧安全評價參考表

實倉檢測數(shù)據(jù)來源于華北地區(qū)某糧庫，實驗倉詳細儲糧信息及采樣方法見方法1.3.2。檢測時間從6月到10月，對應糧食儲藏時間為150 ～ 255 d，溫度范圍為19.0～ 27.4 ℃，水分范圍為12.8% ～ 15.4%，共計84組實倉檢測數(shù)據(jù)，根據(jù)實倉測得的溫度、稻谷含水量和儲藏時間，采用式(2)對實倉儲糧真菌孢子數(shù)進行預測，實測數(shù)據(jù)與預測數(shù)據(jù)結果見表3。

由于微生物生長檢測方法本身誤差很大，對真菌孢子數(shù)量進行預測，預測值與實測值在同一個數(shù)量級以內(nèi)，不影響儲糧安全級別評價，可指導合理的應急處置措施，即認為可以接受。實測數(shù)據(jù)與預測數(shù)據(jù)結果如表2所示，84組實倉數(shù)據(jù)中，有70組數(shù)據(jù)在同一個數(shù)量級內(nèi)，對儲糧安全評價級別在同一級內(nèi)，即對各點儲糧真菌危害程度預測正確率為83.3%。由此可見，通過多元線性回歸方法得到的儲糧真菌生長數(shù)量預測模型為實倉儲糧安全狀況預測提供了一個新的方法和途徑。

圖3 第一次殘差分析結果

表2實倉檢測數(shù)據(jù)與預測數(shù)據(jù)比較

日期6月1日7月15日8月1日8月12日8月26日9月21日10月7日6月1日7月15日8月1日8月12日8月26日9月21日10月7日取樣點真菌孢子數(shù)/105個/g實測值預測值10.00.00.00.30.30.00.30.10.20.50.41.31.21.720.00.00.30.00.00.30.90.20.61.20.91.91.93.630.00.30.31.83.00.02.70.10.30.60.72.50.72.840.00.00.01.51.811.413.80.41.73.83.614.211.620.353.37.211.78.77.27.814.71.53.56.45.75.09.217.660.00.00.00.31.59.39.00.52.02.72.98.25.98.670.00.60.90.05.717.117.70.42.34.13.815.015.122.980.90.30.61.22.72.72.71.22.75.05.04.34.59.890.00.30.00.93.011.719.20.62.43.51.39.99.720.5100.91.53.34.29.624.919.81.35.49.617.632.717.538.0117.88.418.37.57.26.97.82.64.97.010.112.512.616.3121.51.22.75.46.312.015.92.49.816.827.051.959.381.9

本研究通過模擬儲藏條件得到的儲糧真菌孢子數(shù)量靜態(tài)預測模型，但在實際倉儲環(huán)境中，由于糧堆體積較大，不同季節(jié)環(huán)境溫濕度的變化，與實驗室模擬的儲藏條件存在較大差異，模型預測結果仍存在一定偏差。因此，有待于通過大量實倉檢測數(shù)據(jù)反饋驗證，優(yōu)化模型算法形成新的動態(tài)預測模型，模型預測的準確率會進一步提高。

3 結論

以多元統(tǒng)計分析作為基礎與前提，利用MATLAB軟件對1 140組模擬儲藏稻谷的溫度、水分、儲藏時間及取對數(shù)后的真菌孢子數(shù)進行多元線性回歸分析，檢驗、分析多個自變量對因變量的綜合線性影響的顯著性，建立了多元線性回歸模型應用于儲糧真菌生長數(shù)量預測。

經(jīng)殘差分析剔除異常點后得到多元線性模型：

lny=-20.787+0.157×x1+0.027×x2+1.683×x3

多元線性擬合優(yōu)度R2=0.77，F(xiàn)檢驗的統(tǒng)計量值為1 219.552，且顯著性概率為0，結果表明回歸模型總體效果顯著；環(huán)境溫度、儲藏時間和稻谷含水量的t統(tǒng)計量分別為24.785、27.175、46.607，各項參數(shù)對取對數(shù)后的儲糧真菌孢子數(shù)有顯著影響。由多元線性回歸模型的系數(shù)可知，對稻谷儲糧真菌孢子數(shù)量影響最大是稻谷含水量，其次是儲藏溫度和儲藏時間。

通過對多元線性回歸方程進行指數(shù)變換，可得到儲糧真菌孢子數(shù)的預測模型:

y=exp(-20.787+0.157×x1+0.027×x2+1.683 ×x3)

為了探究真菌孢子數(shù)量模型在實際應用中的預測效果，利用實倉數(shù)據(jù)對模型進行初步驗證，以儲糧安全等級作為評價標準，儲糧真菌危害程度預測的正確率為83.3%。但仍有待于經(jīng)過大量實倉檢測數(shù)據(jù)驗證，優(yōu)化算法建立新的動態(tài)預測模型，進一步提高儲糧安全狀況預測的正確率。