基于插值方法的溫室溫度場可視化仿真分析*

肖雪朋，程曼，袁洪波，王起帆

(河北農業大學機電工程學院，河北保定，071000)

0 引言

溫室是設施農業的重要組成部分，其內部小氣候環境因子有溫度、濕度、CO2濃度，光輻射強度等，其中溫度是一個主導因素。溫室內保持適宜的溫度對作物健康生長非常重要，但是溫室內溫度空間分布不均衡，不同位置處的溫度分布差異性大，復雜的溫度分布模式不但導致其對熱能的利用率也相對較低，而且還會使不同區域的作物生長情況產生差異，甚至降低綜合產量。對溫室溫度的空間分布進行可視化仿真，可以直觀的描述溫度的空間分布狀態，對于提高溫室作物管理效率、環境調控和溫室的結構優化的具有重要的現實意義[1]。

近些年，國內外學者進行了大量溫室溫度變化的研究，并建立了豐富的溫室環境模型。總體來說，對于溫室內溫度變化的研究可以分為機理分析法和試驗測量法[2]兩種方法。機理分析法是一種基于理論的方法，機理分析法一般根據質量守恒方程和能量守恒方程，利用計算流體力學(Computational Fluid Dynamics，CFD)技術[3]進行分析和計算；主要用于研究溫室在通風[4-6]、加溫[7-9]、降溫[10-12]等過程中溫度場的變化情況，以及在此基礎上的溫室結構優化改進[13-14]如溫室墻壁厚度及材料選擇[15-18]、溫室跨度選擇[19-20]、排氣孔結構和開口尺寸選擇[21-23]等。機理分析法主要通過數值運算來進行溫室小氣候的模擬和仿真，計算過程繁瑣且需要確定的參數眾多，很多參數很難通過傳感器進行測量，只能通過估算的方式進行近似計算，所以其計算精度直接受到這些參數取值的影響[24]。此外，溫室在長期的使用過程中，有些參數如透光率、圍護結構導熱系數等參數會發生一定的改變[25]，而這些變化難以量化，這些因素都會影響到機理法的計算結果。試驗測量法主要通過各種測量儀器或傳感器進行溫室環境信息的獲取，在此基礎上對溫室的墻體[26]、地面和空氣的熱傳遞[27]，溫室的溫度分布[28]的進行研究。試驗法得到的結果與其獲取信息的傳感器的精度和數量直接相關，為了得到較為精確的模型，需要搭建相對復雜的試驗環境，對于不同的條件，可能還需要多次重復試驗，成本投入大、試驗周期長。將試驗測量和機理分析方法結合，可以簡化溫室溫度場分析的過程，簡單快速根據溫室數據構建模型研究溫室溫度環境的變化規律。

本文利用傳感器采用數據，基于不同的插值方法對溫室溫度場分布進行了計算，并進行了溫度場可視化仿真分析，比較了不同插值方法的優劣及其適用條件，為溫室溫度場分布規律研究提供了一種有益的方法。

1 材料與方法

1.1 試驗地點與時間

試驗地點為河北省保定市河北農業大學東校區，屬于暖溫帶大陸性季風氣候區。試驗用溫室為一面坡型鋼架結構玻璃溫室，坐北朝南，南偏西3°，東西方向長10 m，跨度為2.5 m，前屋面高2.1 m，后屋面高3.3 m。覆蓋物為4mm浮法玻璃，透光率大于83%，地面水泥硬化，栽培模式為槽式基質栽培。溫室主要能量來源為太陽輻射，試驗時間為2019年7—9月，根據室外小型氣象站氣象數據選取典型晴天和陰天環境進行瞬時，短期，中長期的溫度場的繪制。其中，8月11—12號為典型的陰天天氣，溫度為21 ℃～29 ℃，9月1—2號為晴天天氣，溫度為15 ℃～31 ℃，8月8—16號為短期的天氣數據，包含晴天及降溫的天氣類型，8月30號—9月10號為長期的天氣數據，包含多種天氣情況。

1.2 試驗環境構建

溫室內溫度采集使用溫濕度傳感器SHT20(濕度測量范圍：0～100%RH，測量精度：±3%RH，濕度分辨率：0.1%RH；溫度測量范圍：-40 ℃～125 ℃，測量精度：±0.3 ℃，溫度分辨率0.1 ℃)采集數據，測試之前對傳感器進行校正。以溫室東西方向1/2處溫室剖面為主要測量平面，在垂直方向共設置了0 m、0.35 m、1.2 m、2 m四個高度層面，水平方向設置0 m、0.4 m、0.8 m、1.2 m、1.6 m、2 m、2.5 m七個層面。每層均勻分布7個傳感器，溫室頂部均勻分布5個傳感器，共33個溫度傳感器，同時測量1.8 m處太陽輻射強度。玻璃溫室示意圖如圖1所示。溫室內部太陽輻射由JTBQ-2總輻射傳感器(光譜范圍：280～3 000 nm，測量范圍：0～2 000 W/m2，準確度±2%，靈敏度7～14 μV/(W/m2))采集。外界溫度與光照由氣象站采集記錄。

圖1 傳感器部署示意圖

1.3 試驗環境構建

溫度采集與處理功能主要由STM32單片機實現，單片機通過引腳控制傳感器采集數據，然后將信號處理后通過串口輸出環境參數數據。所有測量數據的采集間隔為10 min，數據上傳至電腦儲存。電腦將數據儲存為txt文件，可以實時儲存數據及查閱歷史數據。

1.4 數據分析方法

數據處理采用Matlab2018a，運算平臺計算機的配置為：CPU為AMD R-2600主頻3.8 GHz，內存為DDR4頻率3 000 Hz容量16 G。根據采樣點收集的溫度數據，使用插值算法對數據進行擴充，實現溫度場的擬合和圖像繪制。

插值是在一組已知數據點的范圍內添加新數據點，而且函數運算結果穿過樣本數據點的技術。使用插值算法可以實現填充缺失的數據、對現有數據進行平滑處理以及進行預測等功能。根據樣本數據的結構，插值分為兩類：插入網格數據(meshgrid)，數據部分為軸對齊網格格式的樣本數據，數據格式要求為有序數據。內插散點數據(griddata)，數據部分為散點格式的樣本數據，數據中的點沒有結構或其相對位置間的順序要求。由于網格上的數據插值運算時，數據不能缺失，要求數據為標準的矩陣。也不容易按實際比例顯示空間溫度分布。當數據類型標準時采用網格的數據點和散點數據點插值結果相同。散點插值方便查詢坐標上的數據靈活性更強。綜合考慮采用散點插值的方式。

插值運算過程分為幾個步驟：(1)創建散點數據集。(2)創建Delaunay三角剖分，并查詢插值點位置(查詢點是griddata執行插值的位置)。(3)選擇插值方式計算。通常情況下，整個插值過程通過一次函數調用完成，插值結果曲面始終穿過x、y和z定義的數據點即樣本值v。

散點插值中數據樣本點的樣本值v與坐標中x、y和z空間位置對應，樣本點唯一。以溫室的地面東北角為空間原點，建立傳感器的空間坐標，確定傳感器的空間位置。

在進行散點插值前對數據預處理進行劃分，構建三角剖分數據結構體使數據相對合理的分配以求包圍查詢點。三角剖分是指幾何中將平面對象細分為三角形，并且可以通過類似的擴展將高維度的幾何圖像細分為二維的三角形。劃分結果滿足所有的面均為三角面且沒有相交的邊、三角面內不含數據點(端點除外)、所有三角面組成數據的凸包。

根據插值數據的需求，在劃分的三角形中找到查詢點位置，計算值取決于后續步驟的插值算法。散點插值提供5種運算方式，如表1所示。

表1 散點插值方式區別

1.5 插值運算結果

在MATLAB中調用插值函數將空氣溫度數據3×5差值為6×10，對比實際測量值與插值數據區別。MATLAB編程：vq=griddata(x, y, z, v, xq, yq, zq, method)，method可以根據插值方式選擇‘Linear’、‘Nearest’、‘Natural’、‘Cubic’或‘V4’。

表2 原始數據

Linear線性插值在查詢點插入的值是基于三角網結構體頂點樣本鄰點數值的線性插值。已知兩數據點坐標用直線連接，以求在直線上其他的點進行插值計算。表3為Linear插值數據結果。

表3 Linear插值數據

Nearest最近鄰點插值在查詢點插入的值是距三角網結構體頂點樣本最近的值。已知插值點的坐標與兩數據點間的空間位置，選取與插值點距離近的數據點數值為插值結果。表4為Nearest插值數據結果。

表4 Nearest插值數據

Natural自然鄰點插值在查詢點插入的值是基于三角網結構體臨近點對查詢點的貢獻率計算的插值結果。自然鄰點插值時先對所有樣本點創建泰森多邊形，當對未知點進行插值時,就會修改這些泰森多邊形并對未知點生成一個新的泰森多邊形。表5為Natural插值數據結果。

表5 Natural 插值數據

Cubic三次插值在查詢點插入的值基于三角網結構體頂點樣本鄰點數值的三次卷積插值結果。根據已知數據構造出三次多項式，用構造多項式的極小點逼近尋求已知點函數的極小點的方法，求構造多項式導數為零的根，作為已知點函數的極小點的近似，重復應用這一方法進行迭代計算，直到得出滿足事先給出的精度要求為止，然后進行插值。表6 為Cubic插值數據結果。

表6 Cubic 插值數據

V4插值方式數據劃分非基于三角剖分，采用的是MATLAB自定義的4個樣本點進行樣條函數插值方式。表7為V4插值數據結果。

表7 V4插值數據

Cubic和V4方法生成平滑曲面，而Linear和Nearest分別具有零階導數和一階導數不連續的特性。Nearest插值方式不符合空氣溫度漸變的規律，最后選擇Cubic，Linear，Natural，V4這四種方式對比。

試驗過程：(1)將儲存的環境數據由TXT轉換成Excel文檔。(2)在Matlab編程，導入傳感器空間分布數據,環境數據。(3)設置多核并行運算，插值運算類型，運算精度等參數。(4)輸出溫度場平面圖，運算時間，及數據保存。

2 結果與分析

2.1 插值數據量確定

將溫度數據經過圖像化處理形成溫度分布圖片，可以直觀、迅速、準確地呈現出溫室內部溫度的分布規律，為進行溫度調節研究提供方便。數據量的大小決定著畫面的清晰度和預測數據的準確性，同時影響著運算速度。首先需要確定適宜的插值后數據量，插值數據量可以由式(1)計算得到。

Number=(T_max-T_min)/Scale

(1)

式中：Number——插值數據量；

T_max——進行插值計算時傳感器測量的最高溫度；

T_min——進行插值計算時傳感器測量的最低溫度；

Scale——傳感器分辨率。

選取2019年9月1日溫室的環境數據，同一時刻測量平面最高溫度和最低溫度最大差值20 ℃，根據傳感器分辨率0.1 ℃，由于試驗溫室主要測量面近似方形，最高溫度出現在南側附近和最低溫度出現在北側墻體，橫向、縱向插值產生的200×200數據生成圖片。然后分別生成7×7，25×25，50×50，100×100，400×400的數據量生成圖片，將圖片保存成相同大小。在Matlab中將圖片灰度化化，逐像素點比較相似。結果如圖2不同算法插值數據量與數據可視化后相似度關系曲線，當數據為200×200時四種方式相似度均超過95%。圖3為不同插值方式200×200與400×400灰度圖，從圖片中已經看不到差別，可以認為200×200的數據量滿足圖形繪制要求。

圖2 不同算法插值數據量與數據可視化后相似度關系曲線

2.2 相同天氣下不同插值方式的溫度數據的比較

為了比較不同插值方式的溫度數據的準確性，采用插值與實際采集值對比的方式。9月1—2日為典型的晴天天氣(2019-09-01日最高溫度31 ℃最低溫度15 ℃晴南風微風；2019-09-02日最高溫度31 ℃最低溫度15 ℃西南風2級)，8月11—12日為典型的連陰天天氣(2019-08-11日最高溫度26 ℃最低溫度21 ℃，東北風4級；2019-08-12日最高溫度29 ℃最低溫度21 ℃東北風3級)將數據中各點依次刪除。利用剩余數據對該點進行預測，計算值與測量值差值反映運算方式的準確性，由于刪去1個數據點時對圖片的整體影響過小，采用坐標點取值的方式，取值與采集值做差，得到的數值取絕對值，最后計算每個點不同時間段的平均值，間隔半小時計算一次。數據處理結果如表8所示，晴天數據中最大值61.9 ℃，數據中最小值24 ℃，平均值38.8 ℃，陰天數據中最大值37.0 ℃，數據中最小值25.3 ℃，平均值30.2 ℃。結果顯示：空氣溫度尤其是空間中心位置的插值數據誤差較小，地面、墻體、覆蓋物周圍的溫度數據插值數據誤差較大，可能由于壁面溫度與氣溫相差較多，導致經過插值運算誤差較大。Cubic插值方式誤差最小，V4插值方式誤差最大，Natural和Linear相似。地面附近數據誤差較大。

(a) Cubic插值方式200×200與400×400灰度圖

表8 不同插值算法數據誤差

2.3 不同天氣數據下不同插值方式的溫度數據誤差平均值的比較

誤差平均值如表9所示，晴天時由于陽光充足溫室效應明顯，溫度變化較劇烈。不同部位的溫室壁面接受到的輻射不同，空氣與壁面溫差較大，導致誤差的平均值增大。陰天溫度變化平緩，溫室空氣、壁面溫度在小范圍內變化，誤差平均值小。晴天兩天，一周、兩周、三種環境預測數據誤差平均值相似。總體四種環境預測數據誤差平均值中：Cubic、Natural、Liner誤差平均值明顯小于V4誤差平均值。數據中最小值25.3 ℃，平均值30.2 ℃。結果顯示：空氣溫度尤其是空間中心位置的插值數據誤差較小，地面、墻體、覆蓋物周圍的溫度數據插值數據誤差較大，可能由于壁面溫度與氣溫相差較多，導致經過插值運算誤差較大。Cubic插值方式誤差最小，V4插值方式誤差最大，Natural和Linear相似。地面附近數據誤差較大。

表9 誤差平均值

2.4 不同插值方式運算時間的比較

不同插值方式運算對計算機的資源占用存在差別。Matlab中并行運算可以減少運算時間，但是選擇合適的運算方式也十分重要。采用相同的程序框架，讀取環境數據、傳感器位置數據、對數據加工形成矩陣、設置圖片的輸出樣式、只改變插值方式，比較運算時間。數據量采用200×200大小，運算三次取平均值，運算時間結果如表10所示。結果顯示夜晚溫度差異較小、運算溫度范圍小、計算量小時間短，中午溫室內溫度分布不均，中午溫度差異較大、運算溫度范圍大、計算量大運算時間較夜晚運算時間長。總體運算時間Liner運算時間最短，V4運算時間最長，Cubic和Natural處于兩者之間。

2.5 不同插值方式可視化比較

圖4為溫室不同插值方式溫度場分布。Cubic插值方式結果等溫線平滑，溫度變化平緩。Liner插值方式結果等溫線粗糙，等溫線存在突變情況。Natural插值方式結果處于兩者之間, V4生成的溫度場分布圖與其他三種方法的結果有著顯著的不同，壁面與空氣之間沒有變化趨勢，與實際情況存在較大差距，表明該方法不適用于溫度場插值化的可視表達。

圖5為2019年9月1日玻璃溫室不同插值方式溫度最低時刻溫度場分布。夜晚，沒有加溫時，熱源是地面、墻體的蓄熱。熱損失主要有三部分：通過地面沿橫向傳到室外地面；通過縫隙傳到室外空氣；溫室內的空氣以對流輻射和凝結等方式，向覆蓋材料的內壁面傳導，通過熱傳導方式再擴散到覆蓋材料的外表面，之后再以輻射方式散失到周圍。形成地面、墻體附近溫度高，覆蓋物玻璃附近溫度低溫，室內部大部分空氣溫度均勻的情況。Liner和Natural插值結果相似，靠近墻體部分出現了溫度變化梯度，Cubic插值結果在溫室上部出現了溫度分布溫度變化梯度，符合熱空氣上升的自然規律，V4插值結果溫室范圍內，頂部出現高溫區域不符氣溫變化規律。

(a) Cubic插值溫度場

(a) Cubic插值溫度場

圖6為2019年9月1日玻璃溫室不同插值方式早上日出時刻溫度場分布。早晨太陽升起，溫室效應使溫室內漸漸升溫，升溫從溫室下部地面附近開始。溫室內部的空氣溫度基本相同，內部的溫度分布較為穩定。對比夜晚的溫度數據，壁面與空氣均有明顯的溫度升高。Liner和Natural插值結果相似，地面附近溫度較高，室內空氣大部分溫度均衡，Cubic插值結果在南側出現了部分高溫區域，空氣升溫從南部開始, V4插值結果溫室范圍內，頂部、南部出現高溫區域。

圖7為2019年9月1日玻璃溫室不同插值方式中午溫度最高時刻溫度場分布。

(a) Cubic插值溫度場

(a) Cubic插值溫度場

從圖7可以明顯看出溫室南側靠近玻璃區域為高溫聚集區，而北側墻體的位置溫度較低，說明在中午，南側位置由于優先接收太陽輻射而且輻射較強，空氣溫度明顯高于墻體；溫室上部溫度較低的原因可能是由于溫室上部的建筑物遮擋，影響了光照條件，溫度在水平方向上中間高、北側較低，在豎直方向上部高、地面較低，Cubic、Liner和Natural三種方式插值結果相似，V4上部同樣出現高溫與其他三種方式存在差別。

圖8為2019年9月1日玻璃溫室不同插值方式下午日落時刻溫度場分布。傍晚時，太陽下落太陽輻射減少。溫室內地面、墻體溫度大于覆蓋物的溫度，室內空氣處于降溫的狀態，溫度較為均勻。Cubic插值結果在下部出現了部分高溫區域，空氣接近地面和墻體時溫度較高。Liner和Natural插值結果相似，Liner插值結果出現溫度突變，V4上部墻體同樣出現高溫，中下部出現低溫。

(a) Cubic插值溫度場

3 結論

借助MATLAB平臺高效的數值計算及完備的圖形處理功能，對獲取的大量溫度數據進行批預處理。將數字數據經過圖像化處理轉化為溫度分布圖，是數據可視化方面重要的方式。不同插值運算結果的溫度場研究總結如下。

1) 插入數據量以傳感器精度和溫度差為依據，插值結果可以反應溫室內的溫度場細節，圖片邊緣較為順滑，運算量適宜。

2) 不同的天氣情況對插值的誤差影響較大，與溫度的波動范圍有關。Cubic、Natural、Liner插值運算的誤差結果近似，平均誤差均在1.5 ℃內，且均小于V4。

3) Cubic、Natural、Liner三種插值方式中，Cubic運算時間最長，Natural運算時間其次，Liner運算時間最短，平均每次運算2.816 s，需要大量溫度場圖運算時Liner插值方式較為適宜。

4) Cubic插值運算輸出圖片溫度變化平滑，且具有溫度的水平梯度，接近溫室溫度的變化規律效果好。圖片能夠較好地反映出溫室的溫度場。Cubic插值方式適用于展示研究某一時刻的溫度場。