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

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

(南京財(cái)經(jīng)大學(xué)2，南京 210023)

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

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

1 材料與方法

1.1 模擬儲(chǔ)藏實(shí)驗(yàn)材料及樣品處理

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

1.2 試劑和儀器

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

1.3 實(shí)驗(yàn)方法

1.3.1 模擬儲(chǔ)藏及檢測周期

將不同水分含量(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)箱中模擬儲(chǔ)藏180 d，每10 d取樣一次，檢測儲(chǔ)糧真菌生長數(shù)量。

1.3.2 稻谷實(shí)倉情況及取樣方法

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

圖1 實(shí)驗(yàn)倉扦樣點(diǎn)位置示意圖

1.3.3 儲(chǔ)糧真菌檢測方法

參照LS/T 6132《糧食檢驗(yàn) 儲(chǔ)糧真菌的檢測 孢子計(jì)數(shù)法》。

1.3.4 稻谷含水量測定方法

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

1.3.5 糧溫檢測方法

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

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

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

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

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

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

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

2 結(jié)果與討論

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

2.1 多元線性回歸建模

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

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

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

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

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

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

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

(1)

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

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

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

(2)

2.2 預(yù)測模型的效果分析

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

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

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

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

2.2.2 線性關(guān)系顯著性檢驗(yàn)(F檢驗(yàn))

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

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

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

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

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

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

2.2.4 殘差分析

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

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

2.3 實(shí)倉檢測數(shù)據(jù)驗(yàn)證

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

表1 儲(chǔ)糧安全評價(jià)參考表

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

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

圖3 第一次殘差分析結(jié)果

表2實(shí)倉檢測數(shù)據(jù)與預(yù)測數(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日取樣點(diǎn)真菌孢子數(shù)/105個(gè)/g實(shí)測值預(yù)測值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

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

3 結(jié)論

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

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

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

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

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

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

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