基于糧溫時空相關性的儲糧數量監管方法研究

崔宏偉 吳文福 吳子丹 韓 峰 朱浩天 秦 驍

(吉林大學生物與農業工程學院， 長春 130022)

0 引言

我國糧倉儲量大、分布廣。在我國每次大規模清倉查庫過程中，分散的糧庫使得國家對儲糧數量稽核的工作量較大。由于儲備糧數量較大，若個別糧庫發生虛報庫存、以差換好等違規違法行為，將會造成較大的經濟損失[1-3]。因此儲糧數量監管方法具有一定的研究價值。

目前有關儲糧數量的檢測方法有視頻監控法、稱量法、主動測量法等[1]。視頻監控系統能夠實時監測糧食數量，但監控系統易受供電影響，且監管的視頻采用文件方式管理，安全性不足[1]。稱量法在美國等發達國家的大型糧庫以及我國少量的國家級大型儲備庫中采用[4-5]，但該方法對設備要求較高、維護成本較高[6]。主動測量法包括激光掃描[7-8]等，宋立明等[9]提出了一種借助于激光技術攝像機標定的三維重構法，通過激光測距傳感器測得的距離信息對攝像機進行標定，利用距離信息和激光光斑的二維圖像進行三維坐標計算，從而獲得糧面高度信息。但激光掃描儀價格偏高，大量采購成本高，精度調試也比較困難[10]。除上述方法外，目前本領域學者研究的技術還包括基于壓力傳感器的檢測技術[11-12]、電磁波檢測技術[13]等。張德賢等[14]針對倉內壓強分布隨機性和側面摩擦力的影響，提出一種基于一定分布的糧倉底面和側面壓力傳感器檢測法，其通過引入糧倉數量檢測模型解決了底面壓力傳感器壓強測量值的隨機性的問題，使得系統更加通用，但應用該方法需在糧倉安裝壓力傳感器，增加了檢測成本。秦瑤等[15]提出，利用扦插在糧堆中的電磁波發射天線和接收天線，采用波速層析和衰減層析算法，反演得到糧堆的介電質常數，同時采用速度和衰減聯合層析成像技術，對儲糧高度進行成像，得到實際儲糧體積，從而實現糧倉儲量的估算，但該方法易受糧堆溫度、雜質等因素影響。

針對上述研究存在的問題，本文利用大部分糧庫基本配備的糧情監控系統采集的糧溫數據，通過分析實倉糧溫數據時空序列的相關特征，提出一種基于糧溫時空相關性的糧倉儲量監管方法。根據實倉測溫傳感器分布，提出測溫面、線、點的糧溫自相關與互相關分析方法；對3個不同區域的實倉測溫平面、測溫線、測溫點的糧溫進行自相關與互相關分析，設定異常判定的相關系數閾值；進行儲糧監管試驗驗證。

1 材料與方法

1.1 糧情數據

根據文獻[16]，低溫儲藏為平均糧溫常年保持在15℃及以下，局部最高糧溫不超過20℃的儲藏方式；準低溫儲藏為平均糧溫常年保持在20℃及以下，局部最高糧溫不超過25℃的儲藏方式。文獻[17]定義溫控儲藏包括低溫儲藏、準低溫儲藏和控溫儲藏3種模式(溫控儲藏[18]是指在糧油儲藏期間，采用自然通風、機械通風和制冷通風的降溫技術與倉房氣密壓蓋保溫技術相結合的綠色儲糧工藝)，其中控溫儲藏為平均糧溫不超過25℃，局部最高糧溫不超過30℃的儲藏方式，主要適用于我國南方高溫高濕地區。依據儲藏模式可將我國儲糧區域劃分為：低溫區、中溫區和高溫區[17]。低溫區的冬季時間較長，溫度一般都在0℃以下，每年有5個月(11月—次年3月)平均溫度均在5℃以下；中溫區的冬季時間較短，溫度一般在0℃以上，全年有3個月(12月—次年2月)平均溫度在0～10℃之間；高溫區全年幾乎沒有冬季，平均溫度都在10℃以上，12月—次年1—2月平均溫度在10℃以上。在3個溫度區各選擇一個糧倉近一年的糧溫數據用于方法分析，分別為黑龍江省綏化市某糧倉(簡稱綏化倉)、四川省綿陽市某糧倉(簡稱綿陽倉)、廣東省廣州市花都區某糧倉(簡稱花都倉)。

綏化倉為高大平方倉，尺寸(長×寬×高)為47.5 m×26 m×6 m。倉內裝有78條測溫電纜，共計312個測溫點，測溫點排布為13×6×4，儲藏糧食為玉米，選擇儲糧數據為入倉開始2014-05-06至儲藏中期2015-05-05?；ǘ紓}為高大平方倉，尺寸(長×寬×高)為48 m×26.5 m×6 m。倉內裝有77條測溫電纜，共計308個測溫點，測溫點排布為11×7×4，儲藏糧食為玉米，儲藏時間為2010-01-06至2011-01-31。綿陽倉為小型平方倉，為兩個相鄰廒間，選取其中一個，尺寸(長×寬×高)為24 m×24 m×6 m。倉內設有36條溫度電纜，共計144個測溫點，測溫點排布為6×6×4，儲藏糧食為玉米，儲藏時間為2015-01-15至2016-01-01。糧倉結構與測溫點布置如圖1所示。

圖1 平房倉糧倉結構與測溫點布置示意圖Fig.1 Layout diagrams of temperature measuring points in bungalow

以糧倉底邊頂點作為坐標原點，豎直向上為Z軸正向，糧倉長度方向為X軸正向，寬度方向為Y軸正向，糧倉測溫點坐標系如圖2所示。

圖2 糧倉測溫點坐標系示意圖Fig.2 Diagram of temperature measurement point coordinate system in granary

根據上述坐標系，建立測溫點三維矩陣，矩陣X、Y、Z方向分別有n、m、l個測溫點，則倉內各測溫點坐標(x,y,z)表示為(i,j,k)，其中0

1.2 糧堆異常判定理論基礎與方法

糧堆孔隙小，導熱性能不佳，各測溫點日均溫升較小，因此糧堆中同一位置不同時刻的溫濕度信息存在一定的相關性[19-21]；而糧堆中濕熱信息沿著孔隙傳導，相鄰測溫位置的濕熱信息相互影響，因此糧堆中同一時刻的溫濕度信息也應存在一定的相關性，這種時間與空間的相關性特征稱為糧堆的物理印記，糧倉中傳感器采集的糧堆濕熱信息稱為糧堆的信息印記。結合糧情數據的物理印記與信息印記進行儲糧監管，可使監管方法具有容斷性、容缺性和容錯性。容斷性是指儲備糧的信息印記可以被中斷記錄，但物理印記不會被輕易消除，例如同一批儲備糧的糧溫趨勢保持連貫性，作業會導致糧溫趨勢發生跳躍。容缺性是指儲備糧的信息印記因硬件故障而缺失，可通過插值方法補全。容錯性是指儲備糧的信息印記因硬件或通訊故障而導致錯誤，可以通過插值方法進行修正。根據糧倉濕熱信息物理印記的三大特征，分析糧情濕熱信息在時間與空間序列上的相關性，監測異常的時間與空間位置，即實現儲糧監管。

根據上述相關性原理，本文提出基于糧溫時空相關性的儲糧數量監管方法，具體步驟為：確定糧庫所屬區域，設置異常判定閾值；計算XOY測溫平面糧溫的自相關系數與互相關系數，根據閾值判斷異常發生的日期及平面；計算異常平面測溫線的自相關系數或互相關系數，根據閾值判斷異常測溫線；計算異常線上測溫點糧溫的自相關系數或互相關系數，根據閾值判斷測溫點是否異常；對異常點位進行統計分析，判定異常種類與異常的比例范圍。

參考皮爾遜相關系數計算公式[22]，XOY平面糧溫的互相關系數為

(1)

式中Tijk、Tij(k+1)——第k、k+1層XOY平面糧溫，℃

XOY平面糧溫的自相關系數為

(2)

式中Tijkt1、Tijkt2——第k層XOY平面t1、t2時刻糧溫，℃

XOY平面測溫線糧溫的互相關系數為

(3)

式中Ti(j+1)k——第k層XOY平面內第j+1條測溫線的糧溫，℃

測溫線糧溫的自相關系數為

(4)

測溫點糧溫的互相關系數為

(5)

式中T1t、T2t——相鄰測溫點t時段的糧溫，℃

t——時間序列序號

d——儲糧周期內時間序列的個數

測溫點糧溫的自相關系數為

(6)

圖3 糧倉XOY平面糧溫自相關系數Fig.3 Graphs of autocorrelation coefficient about grain temperature in XOY planes

式中Tt、Tt+1——測溫點在相鄰時間序列t、t+1上的糧溫，℃

2 相關性分析與閾值設定

2.1 傳感器截面糧溫相關性

2.1.1XOY平面糧溫自相關性

糧堆中平行于XOY平面的測溫平面，由上至下分別記為第1、2、3、4層，分別計算糧倉第1、2、3、4層的相鄰兩天糧溫的自相關性。3個糧倉XOY截面糧溫自相關系數變化曲線如圖3a～3c所示。

將圖3a～3c中3個糧倉XOY平面糧溫自相關系數在縱坐標軸方向進行投影，得到自相關系數的分布直方圖，如圖3d所示。由直方圖分析結果統計發現，設置XOY平面自相關系數處于0.8～1區間的個數占自相關系數總數的比例大于97.5%，由此設置XOY平面自相關系數閾值為0.8。

為驗證設置閾值0.8時檢測出異常的準確率，以查準率、查全率對檢測結果進行分析，查全率、查準率計算公式分別為[23-24]

(7)

(8)

式中w——被準確檢測出的異常平面數

W——異常平面數

V——被檢測出的平面數

v——總的被準確檢測出的異常平面數

分析結果顯示綏化倉、花都倉、綿陽倉的查全率分別為100%、100%、91.8%，查準率分別為94.7%、97.3%、97.1%。由分析結果可知，閾值為0.8時綏化倉、花都倉異常的查全率為100%，算法將所有異常檢測出。閾值為0.8時算法的查準率并未達到100%，表明存在誤檢測的異常平面，查詢檢測結果發現被誤檢測的異常平面為靠近糧堆表面，這是因為該層接近糧面，其溫度受倉溫變化影響較大，從而影響該層糧溫前后時間的相關性。綜上所述，設置XOY平面自相關系數閾值為0.8較為合適。

2.1.2XOY平面糧溫互相關性

選取糧堆中平行于XOY平面的測溫平面進行互相關性分析。計算第1層與第2層、第2層與第3層、第3層與第4層糧溫的互相關性，但因各層糧溫的互相關系數波動較大，不利于閾值設定，因此分析當天各平面互相關系數變化率(平面的互相關系數與前一天對應平面互相關系數差值與間隔時間(1 d)的比值)。XOY平面糧溫互相關系數變化率如圖4所示。

將圖4中4層糧溫互相關系數變化率在縱坐標軸方向進行投影，得到互相關系數變化率的分布直方圖。由直方圖統計分析發現，設置閾值區間為[-0.15 d-1,0.15 d-1]時，三糧倉檢測出異常的個數分別為34、35、31，三糧倉的實際異常個數分別為32、35、27，對比檢測異常與實際異常分析發現，此時算法的查全率分別為100%、97.1%、100%，查準率分別為94.1%、97.1%、87.1%。又經統計分析發現，閾值區間為[-0.15 d-1,0.15 d-1]時，查全率最高，查準率較高，因此將XOY平面互相關系數變化率閾值區間設為[-0.15 d-1,0.15 d-1]。

2.1.3縱截面糧溫互相關性

由上文分析可知，第1層糧溫受倉溫影響較大，不利于測溫平面的互相關性分析，因此選擇分析縱截面(XOZ平面與YOZ平面)糧溫的相關性。本文選擇分析XOZ測溫平面糧溫的互相關性。綏化倉XOZ平面共有6層，分析第1與縱中截面(縱中截面糧溫為第3層與第4層糧溫均值)，以及第6與縱中截面的互相關系數?；ǘ紓}XOZ平面共有7層，分析第1與第4截面，以及第4與第7截面的互相關系數。綿陽倉XOZ平面共有6層，分析第1與縱中截面，以及第6與縱中截面的互相關系數?？v中截面的糧溫為相鄰兩截面糧溫均值。糧倉XOZ平面糧溫的互相關系數變化率如圖5所示。

將圖5中3個糧倉XOZ平面互相關系數變化率在縱坐標軸方向進行投影，得到互相關系數的分布直方圖。由直方圖統計分析發現，設置閾值區間為[-0.15 d-1,0.15 d-1]時，三糧倉檢測出異常的個數分別為15、14、16，三糧倉的實際異常個數分別為13、12、13，對比檢測異常與實際異常分析發現，此時算法的查全率為100%，查準率分別為86.7%、85.7%、81.3%。后經統計發現，閾值區間為[-0.15 d-1,0.15 d-1]時，查全率最高，查準率較高，因此將XOZ平面互相關系數變化率閾值區間設為[-0.15 d-1,0.15 d-1]。

圖4 糧倉XOY平面糧溫互相關系數變化率Fig.4 Graphs of correlation coefficient changing value about grain temperature in XOY planes

圖5 糧倉XOZ平面糧溫互相關系數變化率Fig.5 Graphs of correlation coefficient changing value about grain temperature in XOZ planes

2.2 測溫線上糧溫的相關性

根據2.1.3節縱截面互相關系數變化率檢測的糧溫異常日期，從中選擇綏化倉2014-12-10，花都倉2010-06-16，綿陽倉2015-10-26的糧溫進行測溫線的自相關性與互相關性分析。分析各截面平行于X軸的測溫線(測溫點形成的線)自相關性，結果如圖6a～6c所示。

綏化倉平行于X軸的測溫線共24條，花都倉共28條，綿陽倉共24條，由圖6a～6c可以看出，綏化倉相關系數明顯較低的測溫線有兩條，花都倉有一條，綿陽倉有兩條，設置閾值為0.8時可將系數較低的測溫線檢測為異常，因此設置閾值為0.8。經查詢原始糧溫數據發現這幾條線上少部分測溫點前后兩天溫度變化超過4℃，屬于糧溫發熱異常[25]，而其他測溫線兩天糧溫變化正常，由此可見測溫線異常檢測結果與實際異常一致。為驗證測溫線自相關系數閾值為0.8的準確性，根據2.1.1節中的檢測結果，在3個糧倉中分別選擇5個異常日期，對日期內所有平行于X軸的測溫線進行自相關系數分析，分布直方圖如圖7所示。查詢原始糧溫數據發現糧溫異常的測溫線相關系數均小于0.8，共約98%無異常測溫線的相關系數分布于區間[0.8,1]之內(糧溫異常為相同位置相鄰兩天溫度變化超過4℃或當天糧溫局部高于其他位置4℃[25])。因此，測溫線自相關系數閾值0.8較為準確。

圖6 糧倉異常日期測溫線糧溫的自相關系數與互相關系數Fig.6 Autocorrelation coefficient and intercorrelation coefficient about grain temperature in temperature measuring lines of abnormal date

圖7 糧倉異常日期測溫線糧溫的自相關系數分布直方圖Fig.7 Distribution histogram of autocorrelation coefficient about grain temperature in temperature measuring lines of abnormal date

同時分析平行于X軸的相鄰兩條測溫線之間的互相關性，結果如圖6d～6f所示。由圖6d～6f可以看出，3個糧倉測溫線間的互相關系數分布無規律，無法反映糧溫的變化規律，因此測溫線的互相關性不作為異常檢測依據。

2.3 測溫點糧溫的相關性

根據2.2節中檢測的異常測溫線，選擇其上的測溫點進行自相關與互相關性分析。單個測溫點選擇以7 d為周期的糧溫數據進行自相關分析，因此綏化倉選擇12-03至12-09以及12-10至12-16之間的糧溫數據，花都倉選擇06-09至06-15以及06-16至06-22之間的糧溫數據，綿陽倉選擇10-19至10-25以及10-26至11-01之間的糧溫數據，測溫點自相關分析結果如圖8a～8c所示。

圖8 糧倉測溫點的自相關系數與互相關系數Fig.8 Autocorrelation coefficient and intercorrelation coefficient about grain temperature of temperature measurement points

查詢原始糧溫數據發現，綏化倉該測溫線所有測溫點糧溫均異常，花都倉測溫線起始測溫點糧溫異常，綿陽倉除起始位置外其他測溫點糧溫異常。對比異常測溫點與無異常測溫點相關系數，將自相關閾值設為0.8時，由圖8a～8c可以看出，檢測結果與實際異常測溫點相符合。為進一步驗證閾值0.8的準確性，計算上述時間段3個糧倉全部測溫點的自相關系數，如圖8d～8f所示。由圖8d～8f可以看出自相關系數存在負值，查詢綏化倉自相關系數為負值的測溫點，因12-10進行通風使糧溫變化超過4℃，屬于糧溫異常，同時查詢自相關系數低于0.8的測溫點糧溫數據，發現其均受到通風影響，因此測溫點糧溫異常檢測結果與實際異常一致；綿陽倉自相關系數為負值的測溫點屬于異常的測溫線，查詢原始糧溫數據發現，檢測結果與實際異常測溫點一致；同樣地，經查詢原始糧溫數據發現，檢測結果與實際異常測溫點一致。綜上所述，測溫點的自相關系數閾值可設為0.8。

分析單個測溫點糧溫與同一條線上相鄰測溫點糧溫的互相關性，同樣以7 d為周期，糧溫數據范圍同上?；ハ嚓P分析結果如圖8g～8i所示，可以看出，測溫點互相關系數無法反映糧溫變化，這是因為線上測溫點距離較遠，短期內相互之間的影響較弱，因此測溫點的互相關系數不作為異常檢測的依據。

按照上述測溫線與點相關性分析方法，檢測出異常日期內所有異常測溫點。統計當天異常測溫點的總數為e，若e>0.8E(其中E為糧倉總的測溫點數)，則判定為通風、環流熏蒸或糧堆數量異常等；若0.1E≤e≤0.8E，則判定糧堆發熱；若e<0.1E，且異常點離散分布，則判定為糧堆正常。

3 儲糧監管試驗

3.1 糧堆數量變化對溫度的影響分析

糧食儲藏進入穩定狀態后，糧堆內部溫度變化與環境溫度變化相比較為滯后，糧堆內部以及表層溫度均與環境溫度存在一定差異，這種差異在穩定儲藏的過程中基本存在，而且由2.1節分析可知穩定的儲藏過程中糧堆各測溫平面在相鄰時間點上的相關性較高。若某時刻糧堆局部糧食被替換，使得局部糧食溫度與外界溫度差異較小，與更換前的糧溫差異較大，則破壞了穩定儲藏過程中測溫平面在相鄰時間點較高的相關性。為證明糧食更換會對糧溫造成影響，選擇綏化倉第1層2014-11-01至2014-11-30之間的糧溫，模擬第1層糧堆被更換一半。綏化倉地處黑龍江省，此時處于冬季，冬季糧堆內部溫度偏高，更換進入的糧食溫度偏低，因此選擇將2014-12-15至2014-12-30第1層的糧溫更換到11-15至11-30之間。檢測得到糧溫自相關系數如圖9所示，由圖9可以看出，糧食更換時糧溫相關性較低，由此可見更換糧食會影響糧堆溫度。

圖9 糧食數量異常前后自相關系數變化曲線Fig.9 Changing curve of correlation coefficient before and after abnormal reserves

3.2 儲糧監管試驗與分析

為驗證基于糧溫時空相關性檢測方法進行儲糧監管的可行性，進行了儲糧監管試驗。在吉林省長春市朝陽區大嶺試驗基地建立小型糧倉，尺寸為3.8 m×3.8 m×6 m，倉內糧面高4.5 m，于2017-05-11入倉完畢，至2017-08-15儲藏結束，儲藏糧食為玉米，品種為先玉33，入倉時含水率為14%。倉內布置25條測溫電纜，橫、縱方向各5條，測溫電纜間距0.8 m，距倉壁0.3 m，每條電纜上布置4個測溫點，第1層距糧面0.3 m，第4層距倉底0.5 m

(靠近倉底)，中間兩層測溫點均勻布置。試驗倉如圖10所示。同時設計了基于糧溫時空相關性的儲糧數量監管系統，界面如圖11所示。大嶺糧倉面面自相關系數如圖12所示。

圖10 試驗倉Fig.10 Granary

圖11 儲糧數量監管系統Fig.11 Reserves supervision system for grain storage

圖12 大嶺糧倉XOY平面自相關系數曲線Fig.12 Graphs of autocorrelation coefficient about XOY plane in Daling granary

如圖12所示，第1層因靠近糧面，受外溫影響較大，自相關系數波動較大；第2、3、4層在儲藏開始糧堆進入穩定狀態前自相關系數較低，穩定后自相關系數基本大于閾值0.8，在07-01左右相關系數降低，而后07-11左右相關系數再次降低。選擇分析07-01和07-11的測溫線糧溫的自相關性。07-01自相關系數如圖13所示。

圖13 大嶺糧倉2017-07-01測溫線自相關系數Fig.13 Autocorrelation coefficient of temperature line in Daling granary on July 1st，2017

由圖13可以看出，07-01共有11條線的自相關系數低于閾值，因此存在11條糧溫線異常。采用測溫點自相關性分析檢測出的異常的糧溫線，統計糧溫異常點的個數e為35個，0.1E≤e≤0.8E(E=100)，因此判定07-01為糧堆發熱異常。采用同樣的方法分析07-11的糧溫異常點個數為37，因此判定為糧堆發熱異常。經查詢糧倉操作記錄手冊發現，糧倉于06-30通風，外界的熱風進入到糧倉內部，導致糧倉內部溫度升高，因此06-30也可稱為糧堆發熱。糧倉于07-11進行2 h的通風，同樣影響了糧堆溫度信息的相關性。由3.1節分析可知，糧食數量變化影響糧堆溫度，糧溫溫度變化對相關性影響的檢測方法同本節，因此基于糧溫時空相關性的儲糧數量監管方法能夠檢測出糧堆溫度信息突變的異常(糧倉通風、環流熏蒸、糧堆異動)，可滿足儲糧監管的技術需求。

4 結論

(1)提出了基于糧溫時空相關性的糧倉儲量監管方法。首先，分析糧堆測溫截面糧溫的自相關性與互相關性，確定異常發生的日期與平面；然后，利用測溫線的自相關特性確定發生異常的測溫線；最后，根據檢測到的異常測溫線，使用測溫點的自相關特性確定異常的測溫點，實現糧溫異常的定點檢測。

(2)設定了基于糧溫時空相關性檢測方法的糧堆異常判定閾值。在3個溫控儲糧區域(低溫區、中溫區、高溫區)內各選擇一個糧倉進行相關性分析，根據分析結果設定閾值，測溫平面自相關系數閾值為0.8，互相關系數變化率閾值為[-0.15 d-1,0.15 d-1]；測溫線的自相關系數閾值為0.8；測溫點的自相關系數閾值為0.8；分析結果表明，短周期內測溫線與測溫點的互相關系數無法作為異常判定依據。

(3)設計了儲糧數量監管系統，并進行了儲糧監管試驗，試驗結果表明，基于糧溫時空相關性的儲糧數量監管方法能夠檢測出糧堆溫度信息在時間或空間上突變的異常，可滿足儲糧監管的技術需求。