楊士進,張明祥,彭 波,袁三明,李慧璇,廖 瑤

(1.貴州省氣象服務中心,貴州 貴陽 550002;2.貴州省生態(tài)氣象和衛(wèi)星遙感中心,貴州 貴陽 550002)

0 引言

低溫冷害是貴州常見的氣象災害之一[1]。低溫能使茶樹等作物受損,品質降低[2-3],每年造成的損失數(shù)以億計,嚴重制約著貴州農(nóng)業(yè)產(chǎn)業(yè)發(fā)展。因此,開展低溫監(jiān)測對作物區(qū)劃、災害影響評估及保險理賠等都具有重要意義[4-5]。低溫是判斷降溫過程強度的指標[6],常規(guī)的低溫監(jiān)測是根據(jù)氣象站在離地面1.5～2 m處觀測的數(shù)據(jù),依賴于臺站的分布,對于無臺站區(qū)域的數(shù)據(jù)則主要靠插值推算。目前常用的有樣條、反距離、克里金、多項式等插值算法,不同插值算法的適應性不同,不同區(qū)域、不同季節(jié)得到的最優(yōu)插值算法不同[7-8],通常距離氣象站點距離越遠,插值的誤差就越大。由于貴州地形復雜,目前站網(wǎng)間隔大小不一,有的站間隔達11 km[9],不可避免會對插值結果產(chǎn)生影響。

衛(wèi)星遙感具有覆蓋廣、時效性高等特點,成為地面觀測的有力補充。常用的極軌衛(wèi)星傳感器能夠提供高達1 km分辨率的熱紅外遙感數(shù)據(jù),在無云覆蓋的區(qū)域能夠反演獲得無縫的地表溫度(LST)數(shù)據(jù)。地表溫度數(shù)據(jù)可用于估算氣溫,常用的遙感氣溫估算方法在最低氣溫、最高氣溫和平均氣溫估算方面都有應用。氣溫是由地表溫度間接影響的,因此可以由LST參數(shù)化推算而來[10],這是氣溫可以由地表溫度進行估算的理論基礎。常用的推算方案有簡單統(tǒng)計法、高級統(tǒng)計法、溫度植被指數(shù)法(Temperature-Vegetation Index,TVX)、能量平衡方法和大氣溫度廓線外推法等[11]。簡單統(tǒng)計法利用地表溫度參數(shù)進行氣溫估算,此法抓住地表溫度這一影響氣溫的最主要矛盾,即可得到較好效果,其平均標準偏差達1.3～2.0 ℃[12-13],但也有研究表明在作物生長季平均標準偏差達到3.8 ℃[14]。高級統(tǒng)計法利用影響溫度分布的地表類型、經(jīng)緯度、高程等參數(shù),構建多元回歸模型進行氣溫估算,估算的旬平均最高氣溫、月平均氣溫誤差在3 ℃以內(nèi)[15-16]。研究表明,夜間地表溫度和日最低氣溫的一致性最高[14],利用夜間地表溫度和白天地表溫度推算氣溫發(fā)現(xiàn)夜間地表溫度的效果更好[17]。最初的溫度植被指數(shù)法[18]利用地表溫度和植被指數(shù)之間的負相關關系進行氣溫推算,只適用于濃密植被地區(qū),對于有云、水體或植被稀少地區(qū)和復雜山地時誤差較大,而通過去除窗口內(nèi)云和水體可以擴大TVX算法的適用性[19]。能量平衡方法利用能量平衡方程對氣溫進行計算,涉及的參數(shù)計算眾多,挪威中部、西部和東部的中低高山帶與地面數(shù)據(jù)的決定系數(shù)>0.9[20]。大氣溫度廓線外推法利用遙感反演的大氣溫度廓線產(chǎn)品和海拔高程數(shù)據(jù)開展氣溫推算,采用NOAA或MODIS大氣溫度垂直廓線產(chǎn)品推算近地面氣溫的研究較多,誤差多在3 ℃以內(nèi)[21-22]。從以上研究來看,高級統(tǒng)計法物理及統(tǒng)計意義明確,模型復雜度相對較低,但也能獲得精度較優(yōu)的結果,具有較高的實用價值。但高級統(tǒng)計法建立的模型依賴特定地區(qū),推廣到不同時空條件下誤差可能較大,在地形復雜的貴州地區(qū)缺乏效果驗證。

本文以貴州省為例,通過分析MODIS地表溫度、高程等參數(shù)和日最低氣溫的相關關系,建立復雜地形下的日最低氣溫推算模型,并對模型精度進行驗證,以期為低溫指數(shù)保險服務開展提供氣象站點外的數(shù)據(jù)支持。

1 研究數(shù)據(jù)

1.1 數(shù)據(jù)

本文選取2016年2月10日1742個、2017年4月1日1581個、2023年1月31日1864個貴州省自動氣象站小時最低氣溫數(shù)據(jù)作為日最低氣溫推算的氣象站點數(shù)據(jù)(因每個日期具有有效數(shù)據(jù)的自動站數(shù)量不同,故每個日期選取的自動站數(shù)量有差異),站點類型包括國家站和區(qū)域站,分別代表冬季和春季數(shù)據(jù),站點分布見圖1。

圖1 氣象站點分布示意圖Fig.1 Schematic diagram of meteorological station distribution

遙感數(shù)據(jù)選取2016年2月10日、2017年4月1日和2023年1月31日夜間過境的AQUA衛(wèi)星上的MODIS傳感器的L2級日數(shù)據(jù)作為LST數(shù)據(jù)源,空間分辨率為1 km,同時,還選取2016年2月10日、2017年4月1日和2023年1月31日L2級的MODIS水汽數(shù)據(jù)參與分析建模。

一般來說,高程對氣溫的分布具有重大影響,因而選取了ASTER DEM V2數(shù)據(jù)作為高程的數(shù)據(jù)源。

1.2 數(shù)據(jù)預處理

針對MODIS LST數(shù)據(jù),采用NASA官方發(fā)布的MRT軟件進行投影轉換為等經(jīng)緯度投影,將DN值轉換為攝氏溫度,并根據(jù)質量控制字段(QC_Night),只保留LST質量字段值為0、8、17、65、73、81的數(shù)據(jù),這些質量字段代表沒有受到云影響且反演誤差<3 ℃的地表溫度。同時,對于2017年4月1日前一晚的LST數(shù)據(jù),發(fā)現(xiàn)西部的LST數(shù)據(jù)異常,主要表現(xiàn)為大多LST<0 ℃,且比較破碎不連續(xù),考慮是受到云污染導致LST值偏低,所以將其剔除。

對MODIS水汽數(shù)據(jù),除了采用HEG進行投影轉換為等經(jīng)緯度投影外,還將DN值轉換為真實的水汽含量(單位:cm)。

氣象數(shù)據(jù)的質量控制采取異常值剔除法,去掉小時最高、最低、平均溫均為0 ℃的數(shù)據(jù)及最高溫>50 ℃的數(shù)據(jù),通過小時最低溫數(shù)據(jù)計算得到日最低氣溫數(shù)據(jù)。

將高程數(shù)據(jù)重采樣到1 km的空間分辨率,并計算全省的坡度、坡向和地形開闊度數(shù)據(jù),用于后續(xù)建模分析。

將所有數(shù)據(jù)的坐標系統(tǒng)進行統(tǒng)一,并根據(jù)所有氣象站點的坐標,提取相應位置所有遙感柵格數(shù)據(jù)的值,去除每種數(shù)據(jù)的無效值,2016年2月10日、2017年4月1日和2023年1月31日剩余的有效數(shù)據(jù)分別為1369個、780個和1864個。

利用高程數(shù)據(jù),添加經(jīng)緯度字段,通過GIS軟件構建全省1 km分辨率的經(jīng)度、緯度柵格數(shù)據(jù)。

2 日最低氣溫推算模型構建

本文采用多元線性回歸方法開展貴州省最低氣溫推算。首先根據(jù)單因素相關分析,利用一元線性相關分析模型,分別得到所有要素和日最低氣溫的相關系數(shù),篩選出線性相關性顯著的要素。其中2016年2月10日線性相關性高的要素從高至低依次為地表溫度、緯度、高程、經(jīng)度(圖2);2017年4月1日線性相關性高的要素從高至低依次為地表溫度、高程、坡度、經(jīng)度(圖3);2023年1月31日相關性高的要素從高至低依次為地表溫度、坡度(圖4)。均通過0.01的顯著性檢驗。

圖2 2016年2月10日線性顯著通過0.01水平的參數(shù)Fig.2 The linearity significantly passed the 0.01 level of the reference substance on February 10, 2016

圖3 2017年4月1日線性顯著性通過0.01水平的參數(shù)Fig.3 The linearity significantly passed the 0.01 level of the reference substance on April 1, 2017

圖4 2023年1月31日線性顯著性通過0.01水平的參數(shù)Fig.4 The linearity significantly passed the 0.01 level of the reference substance on January 31, 2023

隨機劃分2016年2月10日的171個數(shù)據(jù)作為驗證樣本,1198個數(shù)據(jù)作為建模樣本;2017年4月1日的99個數(shù)據(jù)作為驗證樣本,681個數(shù)據(jù)作為建模樣本;2023年1月31日的188個數(shù)據(jù)作為驗證樣本,1866個作為建模樣本。驗證樣本占比10%～13%,建模樣本占比87%～90%。根據(jù)每個日期得到的顯著參數(shù),構建多元線性回歸模型:

Y=a1x1+a2x2+…+anxn+b

(1)

式中,Y為日最低氣溫,x1、x2、…、xn分別為選擇的參數(shù),a1、a2、…、an、b為擬合系數(shù)。

采用平均絕對誤差(MAE)、均方誤差(MSE)和均方根誤差(RMSE)來評價所建立的推算模型的精度。

3 結果分析

采用多元線性回歸方法,利用不同相關性顯著的參數(shù),逐步增加建模參數(shù),建立冬季和春季的不同線性回歸模型,并用驗證樣本開展精度評價。

3.1 冬季

冬季的樣本有2016年2月10日和2023年1月31日2個時段的數(shù)據(jù),建立的不同推算模型及其精度評價見表1和表2。

表1 2016年2月10日不同線性模型推算最低氣溫驗證精度對比(單位: ℃)Tab.1 Comparison of verification accuracy of different linear models for calculating minimum temperature on February 10, 2016 (unit:℃)

表1對比結果表明:2016年2月10日的數(shù)據(jù),模型6使用地表溫度、緯度、經(jīng)度、高程和水汽5個顯著參數(shù)建立的線性模型精度最高,相對于只使用地表溫度的線性模型,精度有較大提升,其MAE提高了7.3%。

表2對比結果表明:2023年1月31日,使用地表溫度和坡度構建的線性模型精度比只使用地表溫度的線性模型精度更高一些,但是精度提升不到1%,故此數(shù)據(jù)只使用地表溫度進行日最低氣溫度推算。

利用2016年2月10日的數(shù)據(jù)建立多元線性回歸模型如下:

y=1.1501x1-0.1228x2+2.6685x3+0.4423x4+0.3330x5-48.7213

(2)

式中,x1、x2、x3、x4、x5分別為地表溫度、緯度、高程、經(jīng)度和水汽變量,變量名前面的數(shù)據(jù)為擬合的回歸系數(shù), -48.7213為常數(shù)項。利用2023年1月31日的數(shù)據(jù)建立一元線性回歸模型如下:

y=0.7095x1-1.0104

(3)

式中,x1為地表溫度,其中0.7095為擬合的回歸系數(shù),-1.0104為常數(shù)項。

3.2 春季

冬季使用了2017年4月1日的數(shù)據(jù),建立的推算模型及其精度評價見表3。

表3 2017年4月1日不同線性模型推算最低氣溫驗證精度對比(單位:℃)Tab.3 Comparison of verification accuracy of different linear models for calculating minimum temperature on April 1, 2017 (unit:℃)

表3對比結果表明:使用地表溫度、坡度、高程和經(jīng)度4個顯著參數(shù)建立的線性模型精度較高,相對于只使用地表溫度的線性模型,精度也有較高提升,其MAE提高了3.7%。

利用2017年4月1日的數(shù)據(jù)建立多元線性回歸模型如下:

y=0.8029x1-1.127x2+0.0749x3-0.2189x4+24.5729

(4)

式中,x1、x2、x3、x4分別為地表溫度、高程、坡度和經(jīng)度變量,變量名前面的數(shù)據(jù)為擬合的回歸系數(shù),24.5729為常數(shù)項。

3.3 誤差分析

利用模型(2)、(3)和(4),采用經(jīng)度、緯度、地表溫度、水汽含量和高程5種柵格數(shù)據(jù),通過ArcGIS的柵格運算功能,分別計算得到全省日最低氣溫的推算結果(圖5～7)。從圖5中驗證點的平均絕對誤差統(tǒng)計結果來看,2016年2月10日平均絕對誤差<1 ℃、<2 ℃、<3 ℃的點占比分別為45.0%、81.9%和95.9%,平均絕對誤差>3 ℃的點僅占比4.1%。從驗證點的平均絕對誤差空間分布來看,平均絕對誤差>2 ℃的主要分布在中部以南,即日最低氣溫度較高的區(qū)域。從圖6中驗證點的平均絕對誤差分布來看,2017年4月1日平均絕對誤差<1 ℃、<2 ℃、<3 ℃的點占比分別為65.7%、87.9%和99.0%,平均絕對誤差>3 ℃的點僅占比1%。通過冬季和春季的情況來看,平均絕對誤差<2 ℃的驗證點均占比超過80%,>3 ℃的點均較少,誤差較大的點分布在北部。從圖7中驗證點的平均絕對誤差分布來看,2023年1月31日平均絕對誤差<1 ℃、<2 ℃,<3 ℃的點占比分別為33.5%、63.3%、83%,平均絕對誤差>3 ℃的點占比17%,誤差分布沒有明顯的區(qū)域分布特征。

圖5 2016年2月10日的日最低氣溫多元線性模型推算結果及驗證點誤差空間分布Fig.5 Prediction results and verification point error space distribution of the multivariate linear model for daily minimum temperature on February 10, 2016

圖6 2017年4月1日的日最低氣溫多元線性模型推算結果及驗證點誤差空間分布Fig.6 Prediction results and verification point error space distribution of the multivariate linear model for daily minimum temperature on April 1,2017

4 結論

通過上文分析,在利用遙感資料進行日最低氣溫推算的過程中,使用經(jīng)度、緯度、地表溫度、水汽含量和高程5個參數(shù)建立日最低氣溫推算的多元線性回歸模型,其平均絕對誤差(MAE)均在1 ℃左右,其結果比只使用相關系數(shù)最高的地表溫度效果有所提高,表明采用夜間地表溫度參與的日最低氣溫推算在冬季和春季能獲得較高的精度。

根據(jù)線性分析的過程來看,有的參數(shù)與日最低氣溫的線性關系不顯著,但是非線性關系較顯著,因此可考慮構建復雜的非線性模型,預期可以進一步提高日最低氣溫的推算結果。

5 討論分析

為何夜間地表溫度和日最低氣溫的關系如此顯著?文獻[11]闡明了地表溫度和氣溫之間具有明確的物理關系,但是因為有的參數(shù)難以獲取導致無法直接反演,因此利用地表溫度和實測氣溫直接建立相關關系成為可行的途徑。日最低氣溫一般出現(xiàn)在凌晨日出前后,夜間地表溫度因為沒有太陽的輻射能量,主要是地表自身的熱輻射,因此比較穩(wěn)定,而日最低氣溫出現(xiàn)的時間也受太陽輻射的影響較小,所以夜間地表溫度與日最低氣溫的相關性很高。當然,日最低氣溫還受到地形和大氣狀況的影響,所以理論上增加地形和大氣狀況的參數(shù)參與建模能一定程度上提高日最低氣溫的推算效果。通過與其他研究的對比,本文得到的效果較好,如文獻[11]日最低氣溫的RMSE為2.8～4.1 ℃,這可能與本文使用的站點數(shù)量更多有關,也可能因為本文所分析的日期個例較少。

雖然遙感地表溫度推算日最低氣溫的方法具有一定可行性,但由于貴州晴空條件較少,特別是低溫冷害時天氣多陰雨,更難獲得熱紅外遙感地表溫度等地面數(shù)據(jù),限制了此方法的應用范圍,需繼續(xù)開發(fā)不受天氣條件限制的其他推算方法,比如微波反演的地表溫度,才能滿足全天候條件下的實際業(yè)務應用需要。