顧及風向和風速的氣溫空間插值方法

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

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

0 引言

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

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

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

1 研究方法

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

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

1.1 雙線性內插生成風場表面

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

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

(1)

(2)

(3)

(4)

1.2 基于高斯擴散模型構建成本表面

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

(5)

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

由于高斯擴散模型用于模擬污染物從特定來源的擴散，而熱量不能點源化，且不像污染物有穩定的擴散源，故該模型不能直接應用于本研究，但能說明熱量在風場作用下移動成本的起源。將高斯擴散模型中的笛卡爾坐標(x,y)轉換為極坐標后，可以看出濃度會隨距離γ和方位角θ的增大而減小，所以將高斯擴散模型化簡[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網格單元的風速；cos(|DA-DM|)為A網格單元的風速在AB方向上的分量，其值的正負分別代表對下風向和上風向的作用：作用于下風向時，風速越大，站點溫度對下風向溫度貢獻越大，反之則貢獻越小；而作用于上風向時，風速越大，站點溫度對上風向溫度貢獻越小，反之則貢獻越大。

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

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

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

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

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

zi=f(xi)+bTyi+ei

(9)

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

函數f與系數b通過最小二乘法進行估計：

(10)

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

2 方法驗證

2.1 研究區與數據

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

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

2.2 不同協變量的薄盤光滑樣條插值

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

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

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

2.3 誤差分析方法

采用交叉驗證方法[21,22]對109個觀測站點12個月的數據進行插值，運用表征估計值誤差范圍的平均絕對誤差(MAE，式(11))、反映數據估值靈敏度和極值效應的均方根誤差(RMSE，式(12))評估SPOCMC法、DEM法和SPOCMC-DEM法的氣溫插值精度并進行對比分析(表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法，當同時使用SPOCMC和DEM作為協變量進行薄盤光滑樣條插值時效果更優。

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

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

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

2.4 方法適用性分析

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

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

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

2.5 算法復雜度分析

算法的計算耗時與樣本數量直接相關，本文分別采用SPOCMC法、DEM法和SPOCMC-DEM法對江蘇、貴州和山東3省的數據進行氣溫插值，同時計算3種方法在不同樣本數量情況下的執行時間(表3)。可以看出，3種方法的計算耗時隨樣本數量的增加而增大，SPOCMC法對于不同樣本數量的計算時間均明顯大于DEM法，SPOCMC-DEM法的計算耗時最多。

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

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

3 結論

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