劉文毫，衛(wèi)建國,2，張春梅，楊 豫，劉兆宇

(1.北方民族大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，寧夏 銀川 750021；2.中國氣象局旱區(qū)特色農(nóng)業(yè)氣象災(zāi)害監(jiān)測(cè)預(yù)警與風(fēng)險(xiǎn)管理重點(diǎn)實(shí)驗(yàn)室，寧夏 銀川 750002；3.寧夏大學(xué) 農(nóng)學(xué)院，寧夏 銀川 750021)

0 引 言

為對(duì)中寧太陽梁試驗(yàn)區(qū)農(nóng)作物霜凍情況進(jìn)行研究，需要對(duì)研究區(qū)域的溫度分布狀況進(jìn)行推算。目前對(duì)特定區(qū)域的溫度分布狀況進(jìn)行推算研究的成果較多，Doug M. Smith等[1]通過建立的模型對(duì)全球的溫度進(jìn)行了推算，Y.Radhika等[2]采用支持向量機(jī)技術(shù)對(duì)區(qū)域的大氣溫度進(jìn)行了預(yù)測(cè)。國內(nèi)學(xué)者在溫度推算方面也有很深入的研究。例如，唐圣鈞等[3]對(duì)貴州省南北兩個(gè)區(qū)域，結(jié)合多元回歸等方法，分區(qū)建立了基于緯度、海拔和坡向以及坡度因子的氣溫降水模型，在研究區(qū)域進(jìn)行氣候要素的小網(wǎng)格推算，結(jié)果表明緯度因子對(duì)溫度和降水的貢獻(xiàn)最大。衛(wèi)建國[4]根據(jù)溫度分布和熱量理論，在250 m的地理網(wǎng)格資料中，采用多點(diǎn)推算方法，綜合改進(jìn)距離加權(quán)平均法，對(duì)運(yùn)算中的多層推算結(jié)果進(jìn)行權(quán)值訂正，結(jié)果表明經(jīng)過多點(diǎn)推算的溫度圖像平滑而富有變化，從根本上消除了單點(diǎn)推算的臺(tái)階現(xiàn)象。劉靜等[5]利用線性內(nèi)插確定界限溫度的初終日期，利用逐旬氣候資料進(jìn)行小網(wǎng)格訂正，從而得出了紅寺堡等地區(qū)的溫度分布。但是受限于自然條件及其他方面的限制，這些對(duì)溫度推算的研究都是采用傳統(tǒng)實(shí)測(cè)數(shù)據(jù)對(duì)地區(qū)的溫度進(jìn)行分布，而沒有采用地形數(shù)據(jù)結(jié)合經(jīng)緯度信息等資料對(duì)區(qū)域的溫度進(jìn)行推算，在實(shí)現(xiàn)方式上大多是使用傳統(tǒng)的C/C++編程技術(shù)進(jìn)行實(shí)現(xiàn)，沒有采用最新的Python開發(fā)技術(shù)進(jìn)行實(shí)現(xiàn)，因此處理效率相對(duì)較低，無法應(yīng)對(duì)基于大批量數(shù)據(jù)的溫度推算工作。基于此現(xiàn)狀，文中基于GDAL開源庫，以無人機(jī)航測(cè)5 m分辨率DEM數(shù)據(jù)及實(shí)測(cè)的經(jīng)緯度數(shù)據(jù)為數(shù)據(jù)源，實(shí)現(xiàn)了基于小區(qū)域溫度推算模型的溫度分布推算。

1 試驗(yàn)區(qū)概況

太陽梁位于衛(wèi)寧平原之上的中寧縣渠口農(nóng)場(chǎng)，南與中寧縣棗園鄉(xiāng)接壤，以白馬梁為界；北與青銅峽與中寧交界線的胡子溝為界；西倚內(nèi)蒙古阿拉善左旗頭道湖，東靠包蘭鐵路。試驗(yàn)區(qū)地形復(fù)雜多變，由土石山地和緩坡丘陵組成，絕大部分為緩坡丘陵，局地為微坡地。

試驗(yàn)區(qū)處于中溫帶干旱區(qū)，大陸性季風(fēng)氣候特征明顯，光照充足、熱量豐富、蒸發(fā)強(qiáng)烈、晝夜溫差大，氣候干旱少雨，春暖快、夏熱短、秋涼早、冬寒長。主要?dú)庀鬄?zāi)害有干旱、沙塵暴、霜凍、冰雹、暴雨、干熱風(fēng)等。

由于研究區(qū)域氣候特征獨(dú)特，農(nóng)業(yè)災(zāi)害類型多樣，因此該試驗(yàn)區(qū)適合農(nóng)業(yè)氣象災(zāi)害的研究。又因試驗(yàn)區(qū)沒有詳細(xì)的溫度分布數(shù)據(jù)，因此，需要建立一套合理的溫度推算模型，推算該試驗(yàn)區(qū)溫度分布狀況。

太陽梁桃園試驗(yàn)區(qū)布局如圖1所示。

2 數(shù)據(jù)來源及溫度推算模型

項(xiàng)目采用太陽梁試驗(yàn)區(qū)無人機(jī)航測(cè)的5 m分辨率細(xì)網(wǎng)格DEM數(shù)據(jù)作為數(shù)據(jù)源。通常情況下，某地區(qū)的平均溫度分布情況與該地區(qū)大氣環(huán)境、海拔高度、經(jīng)緯度及小氣候環(huán)境有關(guān)[6]。不同季節(jié)太陽輻射[7]和大氣環(huán)流背景是導(dǎo)致氣溫季節(jié)變化的主要原因，同時(shí)，坡度坡向也會(huì)影響該地區(qū)的輻射分布情況[8]。該項(xiàng)目采用前人[3]對(duì)貴州全區(qū)氣候要素的推算模型為基礎(chǔ)，并結(jié)合太陽梁試驗(yàn)區(qū)的具體氣候特征，設(shè)定每項(xiàng)氣候要素的權(quán)重系數(shù)，使其能夠更好地實(shí)現(xiàn)對(duì)太陽梁試驗(yàn)區(qū)的溫度推算。具體推算模型如下：

圖1 太陽梁桃園試驗(yàn)區(qū)布局 (左側(cè)及底部陰影區(qū)為研究區(qū)域)

設(shè)基準(zhǔn)站的平均氣溫為S，則預(yù)測(cè)點(diǎn)的溫度值可表示為：

S=f(λ,φ,h,χ,ω)+T

(1)

其中，S為某預(yù)測(cè)點(diǎn)氣候要素值；λ為預(yù)測(cè)點(diǎn)與基準(zhǔn)點(diǎn)的經(jīng)度差值；φ為預(yù)測(cè)點(diǎn)與基準(zhǔn)點(diǎn)的緯度差值；h為預(yù)測(cè)點(diǎn)與基準(zhǔn)點(diǎn)的海拔差值；χ為預(yù)測(cè)點(diǎn)與基準(zhǔn)點(diǎn)的坡向差值；ω為預(yù)測(cè)點(diǎn)與基準(zhǔn)點(diǎn)的坡度差值。

基于氣象科研人員對(duì)模型中的每個(gè)因子進(jìn)行分析評(píng)估得到各個(gè)因子的系數(shù)，同時(shí)結(jié)合研究區(qū)域2018年3月的基準(zhǔn)點(diǎn)實(shí)測(cè)的平均溫度為14.8℃，確定最終的溫度推算模型如下：

T=[1λ+1φ+0.03h+(-0.01)χ+0.01ω]*

0.05+14.8

(2)

3 關(guān)鍵技術(shù)

3.1 GDAL地理空間數(shù)據(jù)抽象庫

GDAL[9]全稱為地理空間數(shù)據(jù)抽象庫(geospatial data abstraction library)，是眾多操作地理空間柵格數(shù)據(jù)的開源庫之一。它利用抽象數(shù)據(jù)模型來表達(dá)所支持的各種文件格式，并使用一系列命令行工具來進(jìn)行數(shù)據(jù)的轉(zhuǎn)換和處理。GDAL包含兩個(gè)庫：用于操作地理空間柵格數(shù)據(jù)[10]的GDAL和用于操作地理空間矢量數(shù)據(jù)的OGR，且其支持的數(shù)據(jù)格式也相當(dāng)廣泛，包括GeoTiff(.tiff)、Erdas Image(.img)、(.grd)等155種柵格數(shù)據(jù)，以及(.shp)、(.gml)、(.GeoPackage)等95種矢量數(shù)據(jù)。GDAL庫還提供了一系列算法的接口[11]，如矢量柵格化、柵格矢量化、遙感數(shù)據(jù)的空間糾正等，并對(duì)這些算法提供了可以運(yùn)行的文件，方便用戶的使用。利用GDAL和OGR這兩個(gè)開源庫，可實(shí)現(xiàn)對(duì)柵格和矢量數(shù)據(jù)的操作及處理。該項(xiàng)目采用開源GDAL庫讀取了.grd格式的原始DEM文件并轉(zhuǎn)換成通用GeoTiff文件，基于GeoTiff文件提取溫度推算所需坡度、坡向等因子并寫入到柵格文件中。通過調(diào)用GDAL庫操作柵格數(shù)據(jù)的具體流程如圖2所示。

圖2 GDAL庫操作柵格數(shù)據(jù)流程

3.2 Numpy庫及Matplotlib庫

Numpy庫是Python進(jìn)行科學(xué)計(jì)算的基礎(chǔ)模塊，它是一個(gè)提供多維數(shù)據(jù)對(duì)象的Python庫，包含了多種衍生對(duì)象以及一系列為快速計(jì)算數(shù)組而生的例程[12]。Numpy庫最核心的部分是ndarray對(duì)象，通過ndarray對(duì)象可以實(shí)現(xiàn)對(duì)柵格數(shù)據(jù)的有效存儲(chǔ)。同時(shí)，Numpy庫具有強(qiáng)大的矩陣運(yùn)算能力，可以實(shí)現(xiàn)矩陣求逆、求特征值、求特征向量等操作。Numpy科學(xué)運(yùn)算庫為溫度推算工作的實(shí)現(xiàn)提供了支撐。

Matplotlib庫[13]是Python編程語言及Numpy擴(kuò)展包的可視化操作界面，其提供了面向?qū)ο蟮腁PI用于通用GUI工具包。Matplotlib庫類似于MATLAB中的Pyplot模塊。基于Matplotlib庫可將溫度推算中的每個(gè)因子以圖像的形式直觀顯示，并自動(dòng)生成圖例，生成的圖例不僅方便檢查數(shù)據(jù)的準(zhǔn)確性，而且可以幫助靈活地調(diào)整權(quán)重系數(shù)，從而實(shí)現(xiàn)對(duì)溫度分布狀況更加準(zhǔn)確的推算。

4 項(xiàng)目實(shí)現(xiàn)流程

項(xiàng)目首先采用無人機(jī)航拍的細(xì)網(wǎng)格DEM數(shù)據(jù)生成海拔格點(diǎn)數(shù)據(jù)，生成該地區(qū)的坡度格點(diǎn)數(shù)據(jù)、坡向格點(diǎn)數(shù)據(jù)。依據(jù)實(shí)測(cè)的研究區(qū)域經(jīng)緯度范圍，生成包含該區(qū)域每個(gè)格點(diǎn)的經(jīng)度緯度值的格點(diǎn)數(shù)據(jù)。由于生成的格點(diǎn)數(shù)據(jù)結(jié)果為標(biāo)準(zhǔn)矩形矩陣，因此需要使用DEM數(shù)據(jù)生成二值矩陣將格點(diǎn)數(shù)據(jù)裁剪為研究區(qū)形狀。獲取了溫度推算所需5個(gè)因子的格點(diǎn)數(shù)據(jù)以后，基于Python的ConfigParser模塊，讀取了每個(gè)因子的權(quán)重系數(shù)并基于Numpy庫實(shí)現(xiàn)模型的運(yùn)算操作。溫度推算技術(shù)路線圖如圖3所示。

圖3 溫度推算技術(shù)路線圖

4.1 溫度推算因子的提取

(1)海拔因子。

由于DEM數(shù)據(jù)[14-15]已經(jīng)包含研究區(qū)域的海拔數(shù)據(jù)，因此只需對(duì)DEM數(shù)據(jù)進(jìn)行格式轉(zhuǎn)換，將原始.grd格式文件轉(zhuǎn)換為模型推算所需GeoTiff格式。具體步驟如下：通過GDAL.open()函數(shù)打開.grd格式的DEM文件，讀取為Numpy庫支持的ndarray格式，然后通過GDAL庫提供的WriteArray()函數(shù)將數(shù)據(jù)寫入為GeoTiff格式。

(2)坡度坡向因子。

由于GDAL庫包含生成坡度和坡向的函數(shù)，因此該項(xiàng)目通過調(diào)用GDAL庫的函數(shù)生成了包含坡度坡向兩個(gè)因子的文件。坡度坡向因子生成步驟如圖4所示。

圖4 坡度坡向因子生成步驟

生成坡度坡向因子需調(diào)用subprocess模塊，subprocess模塊可以通過命令行的方式執(zhí)行GDAL命令，這是執(zhí)行生成坡度坡向因子命令的前提。引用完成subprocess后開始編寫生成坡度坡向因子的命令，需調(diào)用GDAL庫中g(shù)daldem模塊，然后指定如下參數(shù)：指定需要生成的因子(slope、aspect)；指定輸出圖像的壓縮方式，該項(xiàng)目采用LZW壓縮方式；指定其他參數(shù)，如-s指定計(jì)算比例，-alg指定計(jì)算算法。命令編寫完成后，便可通過subprocess.call()函數(shù)將命令執(zhí)行，分別生成包含研究區(qū)域坡度坡向因子的柵格文件。

(3)經(jīng)度緯度因子。

由于無人機(jī)航測(cè)的DEM數(shù)據(jù)中包含研究區(qū)的經(jīng)緯度范圍，因此使用GDAL庫中的GetGeoTransform()函數(shù)即可獲取研究區(qū)域的經(jīng)緯度范圍，但獲取到的坐標(biāo)是投影坐標(biāo)，還需要編寫函數(shù)將投影坐標(biāo)轉(zhuǎn)換為經(jīng)緯度坐標(biāo)，函數(shù)主要用到的是TransformPoint()方法。最終得到研究區(qū)域的經(jīng)緯度范圍：東經(jīng)105°46′～105°50′；北緯37°38′～37°43′?；诮?jīng)度緯度的區(qū)間，使用Numpy的linspace()函數(shù)分別生成該地區(qū)的經(jīng)度緯度分布矩陣。

4.2 溫度推算模型的實(shí)現(xiàn)

溫度推算模型的實(shí)現(xiàn)需要通過GDAL庫將各個(gè)因子數(shù)據(jù)的GeoTiff文件進(jìn)行讀取，轉(zhuǎn)化為可以進(jìn)行運(yùn)算的矩陣。主要操作如下：

(1)疊加因子權(quán)重系數(shù)：由于每個(gè)因子對(duì)該區(qū)域的溫度影響作用不同，因此需要對(duì)每個(gè)因子疊加不同權(quán)重系數(shù)；

(2)疊加總系數(shù)：每個(gè)因子疊加完權(quán)重系數(shù)后，推算出的溫度區(qū)間可能與實(shí)際溫度分布區(qū)間不符，因此需要設(shè)置總系數(shù)對(duì)溫度分布區(qū)間進(jìn)行調(diào)整；

(3)疊加基準(zhǔn)站溫度T：前兩項(xiàng)運(yùn)算得出的是目標(biāo)點(diǎn)與基準(zhǔn)站溫度的差值，因此需要與基準(zhǔn)站溫度進(jìn)行疊加，從而實(shí)現(xiàn)對(duì)整個(gè)區(qū)域的溫度進(jìn)行推算。

以上三種對(duì)矩陣的操作皆可通過Numpy庫進(jìn)行實(shí)現(xiàn)，這里需要注意的是，在執(zhí)行完每一步操作以后，都要使用Matplotlib庫中的imshow方法將運(yùn)算完成的矩陣以圖形化展示，這樣可以保證在運(yùn)算出現(xiàn)異常值時(shí)及時(shí)發(fā)現(xiàn)，同時(shí)可以對(duì)每一步運(yùn)算結(jié)果進(jìn)行驗(yàn)證，提升溫度推算工作的魯棒性。

基于用戶體驗(yàn)需求，工具在使用時(shí)需要隱藏具體的實(shí)現(xiàn)細(xì)節(jié)，同時(shí)又需要將權(quán)重系數(shù)配置接口開放給用戶。因此引入了Python的ConfigParser模塊，通過此模塊可以實(shí)現(xiàn)讀取配置文件中的權(quán)重系數(shù)實(shí)現(xiàn)溫度的推算，因此用戶僅需在配置文件中配置好每個(gè)因子的權(quán)重系數(shù)即可生成研究區(qū)域的溫度分布圖。

4.3 溫度分布圖的生成

基于Python技術(shù)對(duì)提取到的5個(gè)因子通過溫度推算模型進(jìn)行運(yùn)算以后，最終生成包含每個(gè)網(wǎng)格點(diǎn)溫度值的GeoTiff文件。雖然基于python的Matplotlib模塊可以將溫度推算中的每個(gè)因子進(jìn)行圖形化展示，方便檢查數(shù)據(jù)的準(zhǔn)確性，但是在實(shí)際應(yīng)用中發(fā)現(xiàn)Matplotlib模塊在繪制最終正果溫度分布圖時(shí)存在很多不足的地方，如無法設(shè)置圖例樣式，無法添加指南針比例尺等。因此決定采用ArcMap制圖工具[16-17]制作最終的正果溫度分布圖。ArcMap是一款強(qiáng)大的制圖工具，可以方便地添加圖例、比例尺等組件，同時(shí)可以靈活地根據(jù)數(shù)值范圍設(shè)置不同的配色。

5 實(shí)驗(yàn)結(jié)果

最終生成的正果溫度分布如圖5所示，由圖中不同標(biāo)注的區(qū)域可以看到，研究區(qū)域的西側(cè)區(qū)域溫度較低，而研究區(qū)域北部區(qū)域及東南部區(qū)域溫度偏高，其中溫度最高的區(qū)域?yàn)闁|北部區(qū)域，溫度達(dá)到了15.3 ℃以上。

圖5 正果溫度分布

為了對(duì)項(xiàng)目建立的溫度推算模型進(jìn)行驗(yàn)證，分別在太陽梁試驗(yàn)區(qū)選取三個(gè)觀測(cè)點(diǎn)，對(duì)觀測(cè)點(diǎn)的溫度進(jìn)行逐時(shí)監(jiān)測(cè)，統(tǒng)計(jì)了2018年3月1日至31日的每小時(shí)溫度。具體觀測(cè)位置如圖6所示。

圖6 監(jiān)測(cè)點(diǎn)位置分布

監(jiān)測(cè)點(diǎn)A位于試驗(yàn)區(qū)的北部區(qū)域，監(jiān)測(cè)點(diǎn)B位于試驗(yàn)區(qū)的西南部區(qū)域，監(jiān)測(cè)點(diǎn)C位于試驗(yàn)區(qū)的東南部區(qū)域。通過統(tǒng)計(jì)得出了26日至30日的平均溫度：觀測(cè)點(diǎn)A的平均溫度為15.1℃，觀測(cè)點(diǎn)B的平均溫度為15.13℃，觀測(cè)點(diǎn)C的平均溫度為14.9℃。經(jīng)過與模型推算結(jié)果對(duì)比，得出推算結(jié)果與地區(qū)實(shí)際溫度分布狀況比較吻合，達(dá)到了預(yù)期的目標(biāo)。同時(shí)，溫度推算設(shè)定的五個(gè)因子中，海拔對(duì)溫度的推算影響最大，其次是坡度和坡向?？梢娫谛^(qū)域網(wǎng)格溫度推算中，經(jīng)度和緯度的影響權(quán)重較小，基本上可以忽略不計(jì)，這一點(diǎn)與大區(qū)域的網(wǎng)格點(diǎn)溫度推算不同。

6 結(jié)束語

基于對(duì)溫度推算模型及開源柵格處理GDAL庫的研究，設(shè)計(jì)并實(shí)現(xiàn)了一種基于DEM數(shù)據(jù)的小網(wǎng)格溫度推算模型，并基于開源柵格處理庫GDAL庫及科學(xué)運(yùn)算庫Numpy庫進(jìn)行實(shí)現(xiàn)。首先介紹了研究區(qū)概況及模型的設(shè)計(jì)，其次介紹了實(shí)現(xiàn)模型時(shí)用到的關(guān)鍵技術(shù)，并詳細(xì)介紹了溫度推算模型的實(shí)現(xiàn)步驟，最后對(duì)溫度推算結(jié)果進(jìn)行了分析與驗(yàn)證?；陂_源GDAL庫實(shí)現(xiàn)的溫度推算模型相較于傳統(tǒng)的實(shí)現(xiàn)方式具有可擴(kuò)展性強(qiáng)且運(yùn)算速度快等特點(diǎn)。該成果同樣可推廣應(yīng)用于其他模型的推算工作，如基于該模型中的坡度、坡向、海拔因子，可以應(yīng)用于小區(qū)域風(fēng)場(chǎng)分布情況的推算，這也是該項(xiàng)目接下來的研究方向。