曹端廣，張子民*，王 海，常魯群

(1.山東建筑大學(xué)測繪地理信息學(xué)院，山東 濟(jì)南 250101；2.濟(jì)南市勘察測繪研究院，山東 濟(jì)南 250101)

0 引言

氣溫不僅是自然地域系統(tǒng)界限劃分的關(guān)鍵指標(biāo)，還是陸地環(huán)境過程模擬的重要參數(shù)[1,2]，廣泛應(yīng)用于生態(tài)環(huán)境評價、自然災(zāi)害檢測、農(nóng)業(yè)生產(chǎn)和氣候變化等領(lǐng)域[3-6]，故對氣溫進(jìn)行空間化對于研究全球氣候變暖對生態(tài)系統(tǒng)的影響至關(guān)重要[7，8]。已有研究?；跉庀笳军c實測氣溫數(shù)據(jù)利用空間插值方法(反距離加權(quán)法、克里金法、自然鄰近法、樣條函數(shù)法、趨勢面法等[9])預(yù)測未知點數(shù)據(jù)[10]，例如：Stahl等比較了12種基于回歸和加權(quán)平均的氣溫插值方法，發(fā)現(xiàn)選擇控制站點與計算權(quán)重的方法比計算高程的方法對預(yù)測精度的影響更小[11]；Hutchinson等在樣條函數(shù)插值法的基礎(chǔ)上考慮了海拔對氣候要素的影響，提出了薄盤光滑樣條插值方法(Thin Plate Smoothing Spline，TPS)，并開發(fā)了相應(yīng)的插值軟件ANUSPLIN[12]；Nalder等通過結(jié)合多元線性回歸法和地理加權(quán)回歸法提出了梯度下降反比法(GIDW)，提高了插值結(jié)果的精度[13]；徐翔等通過比較普通克里金法、反距離加權(quán)法以及將高程作為協(xié)變量的薄盤光滑樣條插值法，發(fā)現(xiàn)薄盤光滑樣條插值法適于復(fù)雜山地環(huán)境下氣候要素的空間插值[14]；劉正佳等研究得出，薄盤光滑樣條插值法適用于全國大部分地區(qū)在區(qū)域尺度上以DEM為協(xié)變量的插值[15]。

綜上可知，當(dāng)前氣溫空間插值方法多將高程、經(jīng)緯度等作為影響因素，未考慮氣象數(shù)據(jù)中風(fēng)向和風(fēng)速對氣溫空間插值的影響。我國冬季盛行的西北風(fēng)起源于西伯利亞高壓，其將西北冷空氣吹向東南方，使我國大部分地區(qū)氣溫降低，而夏季則相反[16]。陸福志等通過建立秦嶺—大巴山氣溫與降水格點數(shù)據(jù)集，發(fā)現(xiàn)秦嶺阻擋冬季季風(fēng)南下，從而影響冬季氣溫空間分布[9]，可見季風(fēng)對氣溫影響很大。與高程數(shù)據(jù)不同，風(fēng)向和風(fēng)速數(shù)據(jù)既有大小也有方向，因此，不能將風(fēng)向和風(fēng)速數(shù)據(jù)直接作為薄盤光滑樣條插值的協(xié)變量，需對其進(jìn)行量化。謝瑤瑤等對中國1951-2008年的氣溫站點數(shù)據(jù)進(jìn)行了嚴(yán)格的統(tǒng)計學(xué)分析，證明氣溫服從正態(tài)分布，符合高斯擴(kuò)散模型的假設(shè)條件[17]；Li等利用風(fēng)向數(shù)據(jù)與最短路徑算法提出風(fēng)場最短路徑，并將其代替反距離加權(quán)插值(IDW)中的歐氏距離進(jìn)行空氣污染物空間插值[18]。

本研究在已有氣溫與污染物插值研究的基礎(chǔ)上，提出一種顧及風(fēng)向和風(fēng)速的氣溫空間插值方法，并分別以累計移動成本最短路徑(SPOCMC)、高程(DEM)和SPOCMC-DEM作為協(xié)變量進(jìn)行薄盤光滑樣條插值，通過分析不同方法插值結(jié)果的誤差驗證本文方法的準(zhǔn)確性。

1 研究方法

本文研究方法具體實現(xiàn)步驟(圖1)為：1)根據(jù)氣象站點觀測的風(fēng)向、風(fēng)速數(shù)據(jù)進(jìn)行插值，生成連續(xù)的風(fēng)場表面；2)基于該風(fēng)場數(shù)據(jù)利用高斯擴(kuò)散模型構(gòu)建擴(kuò)散成本表面；3)利用Dijkstra最短路徑算法計算觀測點與待求點的累計移動成本最短路徑(SPOCMC)；4)將SPOCMC作為協(xié)變量進(jìn)行薄盤光滑樣條插值以實現(xiàn)氣溫插值。

圖1 本文方法流程Fig.1 Flowchart of the proposed method in this paper

1.1 雙線性內(nèi)插生成風(fēng)場表面

由于風(fēng)場是矢量數(shù)據(jù)，對風(fēng)向和風(fēng)速進(jìn)行插值生成連續(xù)的風(fēng)場表面，需將矢量風(fēng)分解為東西、南北兩個方向的笛卡爾分量，然后采用雙線性插值法，即在兩個方向分別進(jìn)行一次線性插值[19]。如圖2所示，假設(shè)已知函數(shù)f在Q11=(x1,y1)、Q12=(x1,y2)、Q21=(x2,y1)、Q22=(x2,y2)4個點處的值，為求得未知函數(shù)f在點P=(x,y)處的值，首先在x方向進(jìn)行線性插值，得到式(1)、式(2)，然后在y方向進(jìn)行線性插值，得到式(3),即可得到f(x,y)(式(4))。

圖2 雙線性插值示意Fig.2 Schematic diagram of bilinear interpolation

(1)

(2)

(3)

(4)

1.2 基于高斯擴(kuò)散模型構(gòu)建成本表面

為計算累計移動成本最短路徑，需要求出待求點到已知點的移動成本。本文利用高斯擴(kuò)散模型(式(5))模擬氣溫在風(fēng)場作用下的移動距離，以此作為兩點之間的移動成本，并假設(shè)氣溫受風(fēng)力作用在水平和垂直方向符合正態(tài)分布，風(fēng)力為均勻連續(xù)的作用力[18]。

(5)

式中：C0(x,y,z,u)表示氣溫；x和y分別表示下風(fēng)向距離和待求點與風(fēng)向中心線的水平距離；z表示污染物在擴(kuò)散中的釋放高度，本研究將z設(shè)為固定值；u表示水平風(fēng)速；σy和σz分別表示水平和垂直方向的擴(kuò)散標(biāo)準(zhǔn)差。

由于高斯擴(kuò)散模型用于模擬污染物從特定來源的擴(kuò)散，而熱量不能點源化，且不像污染物有穩(wěn)定的擴(kuò)散源，故該模型不能直接應(yīng)用于本研究，但能說明熱量在風(fēng)場作用下移動成本的起源。將高斯擴(kuò)散模型中的笛卡爾坐標(biāo)(x,y)轉(zhuǎn)換為極坐標(biāo)后，可以看出濃度會隨距離γ和方位角θ的增大而減小，所以將高斯擴(kuò)散模型化簡[18]可得：

CostAB=[F(DA,DM)+F(DB,DM)]×LAB

(6)

F(DA,DM)=[|cos(|DA-DM|)|×VA]-sgn[cos(|DA-DM|)]

(7)

F(DB,DM)=[|cos(|DB-DM|)|×VB]-sgn[cos(|DB-DM|)]

(8)

式中：VA和VB分別代表A、B網(wǎng)格單元的風(fēng)速；cos(|DA-DM|)為A網(wǎng)格單元的風(fēng)速在AB方向上的分量，其值的正負(fù)分別代表對下風(fēng)向和上風(fēng)向的作用：作用于下風(fēng)向時，風(fēng)速越大，站點溫度對下風(fēng)向溫度貢獻(xiàn)越大，反之則貢獻(xiàn)越小；而作用于上風(fēng)向時，風(fēng)速越大，站點溫度對上風(fēng)向溫度貢獻(xiàn)越小，反之則貢獻(xiàn)越大。

1.3 計算累計移動成本最短路徑(SPOCMC)

在移動成本表面的基礎(chǔ)上，采用最短路徑Dijkstra算法[20]計算待求點到已知點的累計移動成本最短路徑(SPOCMC)。算法的輸入條件包括起點與終點在格網(wǎng)中的位置和所有格網(wǎng)單元的有向加權(quán)鄰接矩陣(其權(quán)值為式(5)計算出的移動成本)，算法的輸出結(jié)果為累計移動成本最短路徑，包含最短路徑和最短路徑長度。計算過程如圖3所示，將風(fēng)場表面劃分為規(guī)則格網(wǎng)，網(wǎng)格線上的數(shù)值代表相鄰網(wǎng)格點之間的移動成本，通過式(5)將風(fēng)向和風(fēng)速的矢量數(shù)據(jù)分解到網(wǎng)格的每條邊上進(jìn)行相鄰網(wǎng)格間移動成本的求解，圖中加粗路徑為累計移動成本的最短路徑，其方向與分解到規(guī)則格網(wǎng)后的風(fēng)向一致。

圖3 累計移動成本最短路徑計算Fig.3 Calculation of the shortest path of cumulative moving costs

1.4 薄盤光滑樣條插值法的氣溫插值

薄盤光滑樣條插值法實現(xiàn)了模型的平滑度和精確度的最佳融合，極大提高了插值精度。將累計移動成本最短路徑(SPOCMC)作為協(xié)變量進(jìn)行薄盤光滑樣條插值，以實現(xiàn)顧及風(fēng)向和風(fēng)速的氣溫空間插值。薄盤光滑樣條的理論統(tǒng)計模型為：

zi=f(xi)+bTyi+ei

(9)

式中：zi為位于空間i(i=1,…,N)點的因變量；xi為d維樣條獨立變量矢量；f為要估算的關(guān)于xi的未知光滑函數(shù)；yi為p維獨立協(xié)變量矢量；b為yi的p維向量系數(shù)；ei為期望值為0、方差為ωiσ2(ωi為作為權(quán)重的已知局部相對變異函數(shù)，σ2為誤差方差，在各數(shù)據(jù)點上均為常數(shù)，且一般未知)的自變量隨機誤差。

函數(shù)f與系數(shù)b通過最小二乘法進(jìn)行估計：

(10)

式中：Jm(f)為函數(shù)f(xi)的m階偏導(dǎo)數(shù)，稱為粗糙測次數(shù)或樣條次數(shù)；ρ為正的光滑參數(shù)，用于平衡曲面的粗糙度和數(shù)據(jù)的保真性。

2 方法驗證

2.1 研究區(qū)與數(shù)據(jù)

考慮到氣象觀測點分布以及地形的影響，選擇風(fēng)向和風(fēng)速較為穩(wěn)定、地勢較為平坦的山東省作為試驗區(qū)域。山東省位于中國東部沿海、黃河下游，省內(nèi)中部山地突起，西南、西北低洼平坦，東部緩丘起伏，形成以山地丘陵為骨架、平原盆地交錯分布的地形大勢。實驗數(shù)據(jù)(圖4)包括：1)氣象站點數(shù)據(jù)，從中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)獲取山東省109個氣象站點逐小時觀測溫度、風(fēng)向和風(fēng)速數(shù)據(jù)；2)DEM數(shù)據(jù)，從地理空間數(shù)據(jù)云(http://www.gscloud.cn/)上下載山東省90 m空間分辨率的SRTM DEM UTM 數(shù)據(jù)。

圖4 山東省氣象站點分布Fig.4 Distribution of meteorological stations in Shandong Province

2.2 不同協(xié)變量的薄盤光滑樣條插值

為驗證本文方法可有效提高氣溫空間插值的準(zhǔn)確性，分別以累計移動成本最短路徑(SPOCMC)、高程(DEM)和SPOCMC-DEM作為協(xié)變量進(jìn)行薄盤光滑樣條插值，對比不同插值方法的精度。選取山東省109個觀測站點2019年每月15日9：00的氣溫、風(fēng)向和風(fēng)速數(shù)據(jù)作為實驗數(shù)據(jù)，采用上述3種方法進(jìn)行氣溫空間插值，圖5為7月15日9：00的氣溫空間插值效果圖。

從圖5a可以看出，以高程(DEM)為協(xié)變量進(jìn)行薄盤光滑樣條插值時，研究區(qū)中西部出現(xiàn)氣溫低值點，緣于泰山位于該處，高程較大，氣溫偏低；本文SPOCMC法顧及了風(fēng)向和風(fēng)速影響，插值結(jié)果(圖5b)顯示，受東南風(fēng)影響，泰山低溫區(qū)域較圖5a向西北方向有所偏移，東部沿海低溫區(qū)域也向內(nèi)陸擴(kuò)散；從SPOCMC-DEM法的插值效果圖更能明顯看出受風(fēng)力影響的氣溫擴(kuò)散趨勢，在泰山低溫點處呈現(xiàn)出近似條帶狀的擴(kuò)散路徑。

圖5 空間插值效果Fig.5 Spatial interpolation effect

2.3 誤差分析方法

采用交叉驗證方法[21,22]對109個觀測站點12個月的數(shù)據(jù)進(jìn)行插值，運用表征估計值誤差范圍的平均絕對誤差(MAE，式(11))、反映數(shù)據(jù)估值靈敏度和極值效應(yīng)的均方根誤差(RMSE，式(12))評估SPOCMC法、DEM法和SPOCMC-DEM法的氣溫插值精度并進(jìn)行對比分析(表1)。

(11)

(12)

由表1可知，SPOCMC-DEM法的MAE和RMSE均值分別為0.517、0.779，比SPOCMC法(0.583、1.016)分別降低了11.32%、23.32%，比DEM法(0.809、1.231)分別降低了36.09%、36.72%，表明本文累計移動成本最短路徑(SPOCMC)法的氣溫插值精度要高于DEM法，當(dāng)同時使用SPOCMC和DEM作為協(xié)變量進(jìn)行薄盤光滑樣條插值時效果更優(yōu)。

表1 3種插值方法誤差統(tǒng)計Table 1 Error statistics of three interpolation methods ℃

由圖6可直觀地看出3種插值方法的精度對比結(jié)果，總體上SPOCMC-DEM法的MAE和RMSE最小，其次為SPOCMC法，DEM法最大。3種方法的MAE、RMSE在6-8月略有增大，原因是研究區(qū)夏季氣溫空間差異明顯，導(dǎo)致插值結(jié)果誤差增大;春、秋季氣溫空間差異較小，致使3-4月、9-11月的MAE、RMSE值降低，插值結(jié)果精度提升。

圖6 誤差分析折線Fig.6 Broken lines of error analysis

2.4 方法適用性分析

為驗證本文方法在不同地區(qū)的適用性，選取江蘇省和貴州省兩個具有地形代表性的省份進(jìn)行實驗，其中江蘇位于我國東部沿海地區(qū)，地貌主要為平原，貴州地處我國西南內(nèi)陸腹地，地貌類型復(fù)雜，包括山地、丘陵、高原和盆地。分別選取兩省份2019年3月、6月、9月和12月15日9：00的氣溫觀測數(shù)據(jù)，采用3種方法進(jìn)行氣溫插值并統(tǒng)計誤差(表2)。

表2 3種插值方法在兩個研究區(qū)域的誤差統(tǒng)計Table 2 Error statistics of three interpolation methods in two study areas ℃

從表2可以看出：江蘇的實驗結(jié)果中，SPOCMC法的MAE和RMSE值明顯小于DEM法,這是由于江蘇地處平原，區(qū)內(nèi)高程差異較小，所以高程對氣溫插值造成的影響較小；貴州的實驗結(jié)果中，SPOCMC法的MAE和RMSE值比DEM法大，原因是貴州的地形多為丘陵和山地，高程對氣溫插值影響較大；而使用SPOCMC-DEM法在兩省進(jìn)行插值的MAE和RMSE值均小于SPOCMC法和DEM法，充分證明了添加SPOCMC協(xié)變量可提高氣溫插值的精度。

2.5 算法復(fù)雜度分析

算法的計算耗時與樣本數(shù)量直接相關(guān)，本文分別采用SPOCMC法、DEM法和SPOCMC-DEM法對江蘇、貴州和山東3省的數(shù)據(jù)進(jìn)行氣溫插值，同時計算3種方法在不同樣本數(shù)量情況下的執(zhí)行時間(表3)?？梢钥闯?，3種方法的計算耗時隨樣本數(shù)量的增加而增大，SPOCMC法對于不同樣本數(shù)量的計算時間均明顯大于DEM法，SPOCMC-DEM法的計算耗時最多。

表3 3種插值方法的性能統(tǒng)計Table 3 Performance statistics of three interpolation methods s

設(shè)n為樣本數(shù)量，m為待插值點的數(shù)量。傳統(tǒng)的以DEM為協(xié)變量進(jìn)行薄盤光滑樣條插值需要迭代m次，時間復(fù)雜度為Ο(mn)，空間復(fù)雜度為Ο(n2)。本文算法在傳統(tǒng)方法的基礎(chǔ)上需要額外計算累計移動成本最短路徑，第一步采用雙線性內(nèi)插法生成風(fēng)場表面需要進(jìn)行兩次循環(huán)次數(shù)為n的循環(huán)，時間復(fù)雜度為Ο(n2)，空間復(fù)雜度為Ο(1)；第二步利用高斯擴(kuò)散模型的改進(jìn)公式計算移動成本的時間復(fù)雜度為Ο(n)，空間復(fù)雜度為Ο(1)；最后采用最短路徑Dijkstra算法計算待求點到已知點的累計移動成本最短路徑，其對應(yīng)的時間復(fù)雜度為Ο(n2)，空間復(fù)雜度為Ο(E)(E代表路徑數(shù))。綜上，本文算法的總時間復(fù)雜度為T(m,n)=Ο(mn)+Ο(n2)+Ο(n)+Ο(n2)=Ο(max{mn,n,n2})，空間復(fù)雜度為S(n,E)=Ο(n2)+Ο(1)+Ο(1)+Ο(E)=Ο(max{n2,E})，與傳統(tǒng)方法相比，時間復(fù)雜度和空間復(fù)雜度均有所增加。由于實際氣溫插值工作中樣本點數(shù)量有限，且算法增加的時空復(fù)雜度對于實際工作并不會產(chǎn)生較大影響，所以本文算法仍具有較高的實用性和有效性。

3 結(jié)論

本文提出一種顧及風(fēng)向和風(fēng)速的氣溫空間插值方法，并以SPOCMC、DEM和SPOCMC-DEM分別作為協(xié)變量對山東省109個站點12個月份的氣溫數(shù)據(jù)進(jìn)行薄盤光滑樣條插值比較，結(jié)果表明：以SPOCMC-DEM作為協(xié)變量進(jìn)行插值的MAE、RMSE均值分別為0.517、0.779，略小于以SPOCMC作為協(xié)變量插值結(jié)果的均值(0.583、1.016)，明顯小于以DEM作為協(xié)變量插值結(jié)果的均值(0.809、1.231)，說明該方法相比通常以DEM為協(xié)變量的薄盤光滑樣條氣溫插值方法精度明顯提高，當(dāng)同時以SPOCMC和DEM為協(xié)變量進(jìn)行插值時，效果更優(yōu)。江蘇省和貴州省的插值實驗結(jié)果證明了該方法具有普適性。由于本研究的風(fēng)場模型僅考慮了風(fēng)在水平方向上的作用力，若將垂直方向上的作用力也納入模型中，可能會進(jìn)一步提高插值的精度。目前該算法計算量較大，有待進(jìn)一步優(yōu)化，以降低算法的時空復(fù)雜度。