基于局部薄盤光滑樣條函數(shù)的重慶地區(qū)氣溫空間插值

趙美艷，余 君，胡蕓蕓

(重慶市氣象信息與技術(shù)保障中心，重慶 401147)

早在19世紀(jì)，我國就有了現(xiàn)代方法的氣象觀測記錄，并逐步出現(xiàn)了較為完善的氣象數(shù)據(jù)[1]，但這些數(shù)據(jù)只是離散且不規(guī)則的氣象臺(tái)站數(shù)據(jù)，難以反映空間的連續(xù)變化特征。而站點(diǎn)外的數(shù)據(jù)一般由鄰近臺(tái)站的觀測值用一定的數(shù)學(xué)算法進(jìn)行推算求得，插值算法便是利用已有的采樣點(diǎn)數(shù)據(jù)對未采樣點(diǎn)進(jìn)行估算的一種數(shù)學(xué)方法，被廣泛應(yīng)用于對連續(xù)空間的數(shù)值計(jì)算[2-3]。

插值算法的選擇是數(shù)據(jù)類型和計(jì)算效率的一種平衡，其中任何一種方法都不是絕對的[4]，只有特定條件下的最優(yōu)[5]。近年來，用于氣象要素的空間插值算法有距離權(quán)重法(distance weighting)、克里格法(Kriging)、多項(xiàng)式插值法(interpolating polynomials)、Delaunay三角剖分線性插值、薄盤樣條法(spline methods)等[6-10]，但對于不同的變量所適用的插值方法不同[11]。在對多種插值方法進(jìn)行對比分析時(shí)發(fā)現(xiàn)，基于地質(zhì)統(tǒng)計(jì)技術(shù)的Kriging法和薄盤樣條法較為通用[12-13]。Collins[14 ]用多種插值算法對最高和最低氣溫進(jìn)行了插值效果對比分析，認(rèn)為在不同的時(shí)空尺度下，氣溫的插值誤差估計(jì)是不同的。馮錦明等[15]采用4種內(nèi)插方法對中國160個(gè)臺(tái)站降水觀測資料進(jìn)行空間插值結(jié)果分析，研究認(rèn)為，臺(tái)站分布的密集度對插值方法的選擇有一定影響。對于不同變量，其“最優(yōu)”內(nèi)插法是相對的，而不是絕對的[16-19]。樣條法能夠有效優(yōu)化數(shù)據(jù)逼真度和擬合曲面光滑度之間的平衡，具有不受空間尺度影響、不直接依賴空間平穩(wěn)的協(xié)方差等優(yōu)點(diǎn)；因此在綜合考慮誤差估計(jì)、數(shù)據(jù)結(jié)構(gòu)及計(jì)算簡便時(shí)，使用樣條法進(jìn)行氣候數(shù)據(jù)插值不失為一個(gè)好的選擇[14]。Hutchinson等[20]在利用經(jīng)度、緯度和海拔高度之間線性關(guān)系的基礎(chǔ)上，提出局部薄盤光滑樣條插值算法[21-22]，并根據(jù)氣候要素插值的特點(diǎn)，設(shè)計(jì)編寫了針對氣候數(shù)據(jù)進(jìn)行曲面擬合的專用軟件ANUSPLIN[23]。在ANUSPLIN軟件中允許引入多元協(xié)變量線性子模型，可以平穩(wěn)處理二維以上的樣條，并且能同時(shí)完成兩個(gè)以上表面的空間插值，所以對于時(shí)間序列的氣象數(shù)據(jù)插值尤為適用。

在對氣象要素進(jìn)行空間插值時(shí)，地形是影響誤差的一個(gè)重要因素。氣溫隨高度的上升而下降的現(xiàn)象具有普遍性且這種現(xiàn)象隨著時(shí)間和位置的不同而變化。重慶地處中國西南地區(qū)，地形以山地為主，且坡地面積較大，地形復(fù)雜，本研究將利用基于薄盤光滑樣條函數(shù)的曲面擬合程序ANUSPLIN，并依托數(shù)字高程模型(DEM)實(shí)現(xiàn)對重慶地區(qū)氣溫空間分布模型的建立。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

所用資料為重慶市氣象信息與技術(shù)保障中心提供的2017年12月31日21時(shí)至2018年12月31日20時(shí)重慶1 934個(gè)區(qū)域級自動(dòng)站逐小時(shí)氣溫資料，均經(jīng)過質(zhì)量控制[24]。

1.2 臺(tái)站的選取

為了確保試驗(yàn)數(shù)據(jù)的完整性和可用性，對1 934個(gè)自動(dòng)站進(jìn)行了篩選，選取原則和步驟如下。

(1)柵格挑選。將研究區(qū)域(28°N～32.2°N、105°E～110.2°E)按經(jīng)緯度每0.05°×0.05°為一個(gè)柵格進(jìn)行劃分。若一個(gè)柵格里僅有一個(gè)站，則選取該站；若此柵格里有2個(gè)以上的站點(diǎn)，則進(jìn)行下一步挑選。

(2)計(jì)算所有臺(tái)站的氣溫平均可用率和各臺(tái)站的氣溫可用率，對柵格內(nèi)的站點(diǎn)按可用率進(jìn)行排序，選取臺(tái)站可用率大于平均可用率的站點(diǎn)；若柵格中沒有大于平均可用率的站點(diǎn)，則挑選數(shù)據(jù)可用率最高的一個(gè)站點(diǎn)。

(3)計(jì)算所有臺(tái)站的氣溫平均標(biāo)準(zhǔn)差和各臺(tái)站的氣溫標(biāo)準(zhǔn)差，對柵格內(nèi)的站點(diǎn)按標(biāo)準(zhǔn)差進(jìn)行排序，選取標(biāo)準(zhǔn)差小于平均標(biāo)準(zhǔn)差的站點(diǎn)，若柵格中沒有小于平均標(biāo)準(zhǔn)差的站點(diǎn)，挑選標(biāo)準(zhǔn)差值最小的一個(gè)站點(diǎn)。

在滿足條件(2)或(3)的臺(tái)站中，本研究最終選取了數(shù)據(jù)可用率達(dá)99.9%以上且標(biāo)準(zhǔn)差值相對較小的1 000個(gè)站點(diǎn)進(jìn)行網(wǎng)格化試驗(yàn)。

1.3 數(shù)字高程模型

數(shù)字高程模型(digital elevation model，下簡稱“DEM”)，它是用一組有序數(shù)值陣列形式表示平面坐標(biāo)(x，y)及其海拔高度(z)的一種實(shí)體地面模型，主要描述區(qū)域地貌形態(tài)的空間分布，一般采用連續(xù)等間距的海拔高度點(diǎn)反映地形的變化。氣象要素插值的地形效應(yīng)和空間尺度通常是通過與DEM結(jié)合來實(shí)現(xiàn)的，因此，擁有合適的空間尺度的DEM是構(gòu)造氣象要素空間分布的基礎(chǔ)。本文采用1/20經(jīng)緯度(約5 km)作為插值要表達(dá)的空間尺度而建立與之相對應(yīng)的DEM。地形數(shù)據(jù)來自1∶5 000 000世界數(shù)字地圖。投影方式選用Albert投影。投影范圍為28°N～32.2°N，105°E～110.2°E(重慶范圍)。

1.4 曲面擬合

局部薄盤光滑樣條法在包含普通樣條自變量的基礎(chǔ)上，允許加入線性協(xié)變量子模型，所以它是薄盤光滑樣條原型的一個(gè)擴(kuò)展[25]。如它對氣溫插值時(shí)，可以引入海拔高度等。局部薄盤光滑樣條理論統(tǒng)計(jì)模型如下

zi=f(xi) +bTyi+ei(i=1,…,N)，

(1)

式中，zi是位于空間i點(diǎn)的因變量，xi是樣條獨(dú)立變量的d維向量，f是關(guān)于xi的平滑函數(shù)，yi是獨(dú)立協(xié)變量P維向量，ei是隨機(jī)誤差。當(dāng)式中缺少第二項(xiàng)，即模型無協(xié)變量時(shí)，該模型就變?yōu)橐粋€(gè)普通的薄盤光滑樣條模型。當(dāng)缺少第一項(xiàng)獨(dú)立自變量時(shí)，模型便變?yōu)橐粋€(gè)多元線性回歸模型。

最早的擬合程序通常需要至少兩個(gè)獨(dú)立樣條變量，(即f(xi) 中i為2維矩陣)，一般是經(jīng)度、緯度(以度為單位)。但是在擬合氣溫或降水量時(shí)，可增加第三個(gè)獨(dú)立變量，即海平面以上的高程(海拔高度)。在擬合多變量氣象表面時(shí)，只需知道樣點(diǎn)處的獨(dú)立變量的值，因此，氣象站點(diǎn)的坐標(biāo)和海拔信息必須準(zhǔn)確。坐標(biāo)或海拔信息錯(cuò)誤的點(diǎn)會(huì)在輸出的最大殘差日志里反映出來，即以降序排列的殘差文件中，排在首位的幾個(gè)極大殘差值對應(yīng)的站點(diǎn)，可用于檢驗(yàn)原始數(shù)據(jù)在位置和數(shù)值上的錯(cuò)誤。

ANUSPLIN在插值過程中逐步迭代產(chǎn)生一系列統(tǒng)計(jì)參數(shù)，用來判斷插值效果。如表征擬合曲面復(fù)雜程度的信號自由度Signal值需小于站點(diǎn)數(shù)的一半，且在以月為單位進(jìn)行曲面擬合時(shí)，Signal值應(yīng)有較平穩(wěn)的月間過渡；廣義交叉驗(yàn)證GCV(generalized cross validation)估算插值誤差是通過移去一個(gè)站點(diǎn)，用剩余站點(diǎn)進(jìn)行曲面擬合時(shí)得到該點(diǎn)的估算值，從而計(jì)算該點(diǎn)原始觀測值與估算值之間的誤差；GCV的平方根(RTGCV)是由輸入數(shù)據(jù)誤差和估算誤差組成，在模型選取時(shí)，應(yīng)確保RTGCV是最小的；期望真實(shí)均方誤差(RTMSE)是所有樣點(diǎn)的預(yù)計(jì)均方根誤差的估算，相當(dāng)于插值過程的真實(shí)誤差，同樣要選擇RTMSE值最小的模型；另外，RTGCV和RTMSE的差值越大,可間接說明模型的解釋率越高。

1.5 獨(dú)立性檢驗(yàn)

為驗(yàn)證ANUSPLIN所選方案的插值精度，本文除采用ANUSPLIN自帶的統(tǒng)計(jì)誤差進(jìn)行分析外，還將基于重慶范圍內(nèi)未參與插值的35個(gè)國家級自動(dòng)站，采用交叉驗(yàn)證和相關(guān)分析兩種方法對插值結(jié)果進(jìn)行精度檢驗(yàn)。平均絕對誤差(MAE)和均方根誤差(RMSE)可以作為衡量估算值與真實(shí)值誤差的兩個(gè)重要指標(biāo)，即MAE和RMSE值越小，表明插值效果越好。

(2)

(3)

式中，n為臺(tái)站數(shù)，Toi和Tei分別表示第i個(gè)臺(tái)站的觀測值與估算值，同時(shí)，還計(jì)算了相關(guān)系數(shù)來反應(yīng)臺(tái)站的估算值與觀測值之間的相關(guān)性。

2 結(jié)果與分析

2.1 插值表面分析

時(shí)間序列的氣象要素空間插值結(jié)果既要能保證插值表面的插值精度，又要保證所選插值模型的穩(wěn)定性，使其在時(shí)間和空間的連續(xù)上具有可比性。為尋找合適的氣溫插值方案，本研究共設(shè)計(jì)了6種模型(表1)，即以高程數(shù)據(jù)為自變量或協(xié)變量，改變樣條次數(shù)。

表1 薄盤光滑樣條函數(shù)模型

針對重慶市1 000個(gè)站點(diǎn)氣溫要素的空間插值，在參照模型判別標(biāo)準(zhǔn)的條件下，當(dāng)Signal值小于站點(diǎn)數(shù)的一半時(shí)，選取模型最穩(wěn)定，且GCV值最小的方案，經(jīng)過反復(fù)試驗(yàn)，最終確定以經(jīng)、緯度為函數(shù)自變量，海拔高度為協(xié)變量，樣條次數(shù)為2的三變量局部薄盤光滑樣條函數(shù)。

圖1給出了2018年8月1日10時(shí)重慶的氣溫插值。從圖1a可以看出，插值表面帶有明顯的地帶性差異。就整個(gè)重慶來看，中西部氣溫明顯高于東北及西南地區(qū)，其中重慶東北部有一白色區(qū)域，氣溫明顯低于其他地區(qū)，主要因?yàn)檫@里海拔高度較高(2 500 m)，平均氣溫值比周邊低8～10 ℃。值得注意的是，在重慶中東部地區(qū)，有幾條明顯的條帶狀氣溫低值區(qū)，如梁平、墊江、萬州及忠縣等地，而這些帶狀低值區(qū)正好對應(yīng)著明月山、精華山等山脈；因此從圖中可以明顯看出氣溫隨高度的梯度變化，這與常見的氣溫插值趨勢面不太一樣。從估算標(biāo)準(zhǔn)誤差(圖1b)可以看出，整個(gè)重慶的氣溫誤差均較小，誤差值基本在1.0 ℃以下，而重慶地區(qū)以外，誤差值逐漸增大。就重慶內(nèi)部而言，東北及東南部的高海拔地區(qū)誤差比其他地區(qū)偏大0.1 ℃左右，因?yàn)楦吆０螀^(qū)，站點(diǎn)相對較少，從而導(dǎo)致誤差稍大。

圖1 2018-08-01T10重慶氣溫插值(a)和估算標(biāo)準(zhǔn)誤差(b)(單位為℃)

由此可以看出，引用高程線性子模型的局部薄盤光滑樣條函數(shù)可以較好實(shí)現(xiàn)對氣溫的空間插值，且能實(shí)現(xiàn)對站點(diǎn)稀少的山脈地帶氣象要素的插值估算，而插值誤差因地形的差異會(huì)有不同表現(xiàn)，即站點(diǎn)稀少高海拔區(qū)相對于站點(diǎn)密集低海拔區(qū)，估算誤差較大。

2.2 統(tǒng)計(jì)結(jié)果分析

月平均氣溫的插值曲面統(tǒng)計(jì)分析結(jié)果見表2。從表中可以看出，信號自由度Signal值遠(yuǎn)遠(yuǎn)小于站點(diǎn)數(shù)的一半，由此可以說明試驗(yàn)所用站點(diǎn)數(shù)能夠滿足插值的需求。氣溫插值的期望真實(shí)均方誤差RTMSE值除7、8月份大于0.2 ℃外，其余多數(shù)月份均小于0.2 ℃，且2018年各月RTGCV值的大小分布也表現(xiàn)出了秋冬季較小，夏季較高的分布形式。夏季，重慶中西部地區(qū)(低海拔區(qū))高溫悶熱，氣溫高達(dá)40 ℃，而東部高海拔區(qū)的氣溫最高在30 ℃左右，氣溫的空間分布差異較大；冬季，重慶高、低海拔區(qū)的氣溫差異相對夏季來說則較小。由此可以看出，重慶復(fù)雜的地形(海拔差異較大)對氣溫空間差異的影響夏季較冬季明顯。

表2 2018年各月平均氣溫插值統(tǒng)計(jì)結(jié)果

由于模型中引入了第三變量，即海拔高度作為協(xié)變量，因此便存在一個(gè)隨高程變化的線性常數(shù)，ANUSPLIN在此提供了一個(gè)氣候變量隨海拔高度的變化率(lapse rate)。從圖2可以看出，氣溫隨海拔高度下降的幅度在夏季為0.6 ℃/100 m，春秋季較小為0.5 ℃/100 m左右，冬季最小，為0.4 ℃/100 m左右。由此可看出，不同的季節(jié)，氣溫隨海拔高度的變化率并不完全相同，這跟一些學(xué)者研究其他地方得出的結(jié)論相似[26]。

圖2 重慶氣溫隨海拔高度變化率的月際變化

2.3 獨(dú)立性檢驗(yàn)

為了驗(yàn)證模型所選插值方案對氣溫的插值精度，將重慶范圍內(nèi)未參與插值的35個(gè)國家級自動(dòng)站的氣溫觀測值與模型插值結(jié)果，求取平均絕對誤差(MAE)和均方根誤差(RMSE)(圖3)。整體上看，所用插值方案插值效果較好，月平均氣溫的MAE值為0.69 ℃，且冬季優(yōu)于夏季，其中1月最小(0.60 ℃)，9月最大(0.85 ℃)。RMSE值隨時(shí)間的分布與MAE相似，冬季相對較小。雖然獨(dú)立檢驗(yàn)的插值均方根誤差RMSE相對于模型本身計(jì)算的期望真實(shí)均方誤差RTMSE稍偏大(這或許跟模型考慮了地形因素有關(guān))，但二者隨時(shí)間的分布特征相似。另外，插值月平均氣溫值與臺(tái)站觀測值的相關(guān)系數(shù)達(dá)到0.995，相關(guān)性較高。由此可以看出，本研究所采用的插值方案，即以經(jīng)、緯度為函數(shù)自變量，海拔高度為協(xié)變量，樣條次數(shù)為2的三變量局部薄盤光滑樣條函數(shù)對重慶地區(qū)的氣溫插值較為適用。

圖3 重慶插值氣溫的平均絕對誤差(a)和均方根誤差(b)

3 結(jié)論與討論

(1)利用薄盤光滑樣條函數(shù)的曲面擬合程序ANUSPLIN和依托數(shù)字高程模型(DEM)，以經(jīng)、緯度為函數(shù)自變量，海拔高度為協(xié)變量，樣條次數(shù)為2的三變量局部薄盤光滑樣條函數(shù)作為插值方案，建立重慶地區(qū)氣溫要素的空間分布模型，實(shí)現(xiàn)了對重慶市1 000個(gè)站點(diǎn)氣溫的最優(yōu)空間插值。

(2)從氣溫插值結(jié)果可以發(fā)現(xiàn)，插值方案實(shí)現(xiàn)了對站點(diǎn)稀少的高海拔區(qū)氣溫要素較為精確的插值估算，且插值表面能夠明顯看出氣溫隨高度的梯度變化，再現(xiàn)了地形因素對氣溫空間差異的影響在夏季較冬季明顯的特征。由此可以看出，研究所采用的方案對重慶地區(qū)的氣溫插值是適用的。