基于Kriging的加密自動氣象站要素場插值與改進(jìn)

王興　張琳焓　王麗娟　張瀟瀟

摘 要：為了解決加密自動氣象站資料的時空分辨率高，對資料中各個要素場插值計算時耗長，難以應(yīng)用到對實時性要求高的業(yè)務(wù)系統(tǒng)中的問題，提出一種改進(jìn)的Kriging算法，針對加密自動站的空間分布和氣象學(xué)特點，預(yù)先對自動站分區(qū)，并對要素場插值的算法流程進(jìn)行優(yōu)化，提高其計算速度，以滿足相關(guān)業(yè)務(wù)對氣象數(shù)據(jù)實時性的要求。實驗結(jié)果證明，本方法插值結(jié)果更接近實際觀測值，且計算速度提升了約20倍。

關(guān)鍵詞：克里金；插值；加密自動氣象站；模式初始場

中圖分類號：TP391.4 文獻(xiàn)標(biāo)識碼：A

1 引 言（Introduction）

我國是世界上受氣象災(zāi)害影響較大的地區(qū)之一，隨著社會、經(jīng)濟(jì)的發(fā)展，氣象災(zāi)害的影響越來越廣，造成的損失越來越大，所受到的關(guān)注也越來越多。據(jù)不完全統(tǒng)計，近年來我國每年因短時強降水、雷雨大風(fēng)等極端天氣所造成的經(jīng)濟(jì)損失高達(dá)數(shù)千億元，受影響人口超過4億人次[1]。因此，國家對氣象服務(wù)，特別是對災(zāi)害性天氣的預(yù)報、預(yù)警以及防災(zāi)減災(zāi)服務(wù)的要求越來越高。為了提高氣象監(jiān)測、預(yù)測的能力，國家投入了大量資金推動和促進(jìn)大氣監(jiān)測自動化的進(jìn)程[2]，自2000年開始，至今已投入業(yè)務(wù)化運行的加密自動氣象站（以下簡稱自動站或站點）數(shù)量超過4萬部，且仍在擴充建設(shè)。

自動站是一種適用于中尺度區(qū)域地面氣象觀測的專業(yè)監(jiān)測設(shè)備，可將溫度、濕度、氣壓、風(fēng)向、風(fēng)速等多種氣象要素的實況監(jiān)測數(shù)據(jù)實時發(fā)送到氣象業(yè)務(wù)部門。在我國中東部及沿海地區(qū)，約每5—10公里就有一套自動站系統(tǒng)，監(jiān)測數(shù)據(jù)發(fā)布的頻次一般為5分鐘，密集時可達(dá)到每分鐘一次。除了提供高時空分辨率的氣象實況觀測資料外，自動站數(shù)據(jù)還可應(yīng)用于中尺度數(shù)值天氣預(yù)報模式和災(zāi)害性天氣預(yù)報預(yù)警等業(yè)務(wù)系統(tǒng)中[3，4]，用以提供更加精細(xì)、準(zhǔn)確、客觀的模式初始場資料，進(jìn)而提高數(shù)值模式預(yù)報的水平。

由于自動站在空間分布上是若干個離散的點，而模式預(yù)報等氣象業(yè)務(wù)系統(tǒng)所需的初始場資料主要是等間距的格點場，因此，自動站觀測數(shù)據(jù)需要經(jīng)過插值處理。常用的插值算法有雙線性插值[5]、反距離權(quán)重插值[6]、樣條插值，以及Kriging（克里金）插值[7-10]等，其中反距離權(quán)重插值因其算法簡單、易于實現(xiàn)，且考慮了距離對權(quán)重的影響，應(yīng)用較為廣泛[11-13]。然而，各個氣象要素在空間上的分布并非均勻、平滑的，其受到地形、高程以及大氣運動的影響非常大。Kriging算法也稱為空間局部估計法，是在變異函數(shù)理論和結(jié)構(gòu)分析的基礎(chǔ)上，在有限區(qū)域內(nèi)對區(qū)域化變量的取值進(jìn)行無偏、最優(yōu)估計的一種方法[14]。該算法的插值結(jié)果往往更接近于實際值，從而能夠更加真實地反映各個氣象要素的空間分布情況。然而，Kriging算法的運算非常復(fù)雜，計算耗時長，對于更新頻次高、觀測要素多的自動站數(shù)據(jù)，難以滿足現(xiàn)代氣象業(yè)務(wù)系統(tǒng)的實時性要求。

本文針對自動站數(shù)據(jù)在空間分布上的特點，設(shè)計了一種方法來降低Kriging插值的運算量和計算的時耗，以滿足自動站數(shù)據(jù)插值相關(guān)應(yīng)用實時業(yè)務(wù)化的要求。

2 算法流程設(shè)計（The process design of algorithms）

2.1 自動站的分塊處理

在插值算法中，每一個待插值點都需要搜尋若干個與其位置相鄰的實測點，通行的方法就是對所有實測點進(jìn)行逐點遍歷。以自動站數(shù)據(jù)為例，假設(shè)有40 000個有效自動站數(shù)據(jù)，目標(biāo)格點數(shù)據(jù)的空間分辨率為10km×10km，則能夠覆蓋我國所有地區(qū)的矩形所含格點總數(shù)約為：

G1=經(jīng)向跨度×緯向跨度×單位經(jīng)向公里數(shù)×單位緯向公里數(shù)/目標(biāo)格點數(shù)據(jù)經(jīng)向分辨率/目標(biāo)格點數(shù)據(jù)緯向分辨率≈|35-73|×|53-3|×110×110/10/10=375 100。

如果對這G1個格點逐一搜尋40 000個自動站中最鄰近的n個用以插值計算，那么其運算規(guī)模達(dá)到1.5×1010?？紤]到每次搜尋相鄰的自動站時，絕大部分站點因相距待插值點甚遠(yuǎn)，沒有遍歷的必要，本文則使用如下方法改進(jìn)。以1×1的經(jīng)緯網(wǎng)格為單位，將遍布于全國各地的自動站分別映射到相應(yīng)地理位置的經(jīng)緯網(wǎng)格中。即創(chuàng)建一個二維數(shù)組G[X，Y]，對于任意一個待插值格點g[i，j]，g[i，j]∈G[x，y]，其中X為我國領(lǐng)土東西徑向跨度的值取整加1，Y為南北緯向跨度的值取整加1，x∈[0，X]，y∈[0，Y]。每計算一個格點時，只查找以其為中心，上下左右最相鄰的9個單位經(jīng)緯網(wǎng)格中所有自動站的記錄。如此將每個待插值格點遍歷自動站的數(shù)量大幅降低。

2.2 變異函數(shù)的計算

變異函數(shù)是指區(qū)域化變量Z（x）和Z（x+h）的差值的平方的數(shù)學(xué)期望，一般標(biāo)記為，其中，h為這兩個變量之間的距離。為了便于算法的工程化，假定與自變量x無關(guān)，即變異函數(shù)的值不隨Z（x）的位置而變化，僅與h有關(guān)，亦即函數(shù)滿足二階平穩(wěn)假設(shè)[15]，則變異函數(shù)的計算可簡化為如下式：

（1）

式（1）中，E表示數(shù)學(xué)期望，n為樣本數(shù)據(jù)的對數(shù)，xi表示樣本點的值，習(xí)慣上將稱為變異函數(shù)，相應(yīng)地稱為半變異函數(shù)，它們反映的是Z（x）與Z（x+h）隨距離變異的程度，本節(jié)以下內(nèi)容即是對公式1計算過程的描述。

由于自動站在全國的分布并不規(guī)則，總體上呈東南部多西北部少，同時，各氣象要素在空間分布上具有一定的氣象學(xué)規(guī)律，如發(fā)生天氣過程或冷暖鋒過境時，其所影響的地區(qū)的溫度、氣壓等要素在空間上可顯現(xiàn)一定的梯度變化特征和變化趨勢方向性，針對這些特征，本文在計算變異函數(shù)時，分為以下幾步進(jìn)行：

（1）按方位角將自動站分組。將每次參與插值計算的9個單位經(jīng)緯網(wǎng)格中所有自動站，以當(dāng)前待插值格點為圓心，正北方向為0，逐45劃分為8個扇形區(qū)域，則每個自動站必在且僅在一個扇形區(qū)域中。

（2）按距離將自動站分組。將上一步按扇形區(qū)域分組后的自動站再按各自與圓心間的距離分組。單位距離按實際情況可調(diào)整。距離的計算方法為：

（2）

其中，g[i，j]表示待插值格點，Z[x，y]表示自動站，lat和lon分別表示格點或自動站的緯度和經(jīng)度。如圖1所示，每個圓點代表1個自動站，每個自動站都落在某一扇區(qū)中，每個小方格覆蓋的地理范圍為一個經(jīng)緯度網(wǎng)格。其中，，d為單位距離的整數(shù)倍。

圖1 自動站按角度和距離分組示意圖

Fig.1 A schematic diagram of AWS grouped

by angle and distance

（3）計算的主要流程如下：

a.確定變異函數(shù)的方向數(shù)，即劃分的扇形區(qū)域的數(shù)量有N個。

b.確定距離步長為L，步數(shù)為M。

c.假定參與一次插值的自動站數(shù)量不超過K個。

d.定義三維數(shù)組Value[M，N，K]，用于存放自動站記錄，由于每條記錄包含若干個氣象要素信息，因此Value的每一個元素項又是一個一維數(shù)組。

e.定義數(shù)組Sum[M，N]，作為計數(shù)器，Value[M，N，K]每新增一條記錄，相應(yīng)的Sum[M，N]加1。

f.選定一個待插值格點g[i，j]，遍歷所在地理位置相對應(yīng)的9個單位經(jīng)緯網(wǎng)格中的自動站，直到所有自動站信息都記錄到Value[M，N，K]和Sum[M，N]中。

g.開始計算變異函數(shù)的值，計算方法為Value數(shù)組元素中各個氣象要素各自的累加值除以相應(yīng)Sum元素的值，即求取逐個扇區(qū)中，各個自動站要素的平均值。

2.3 變異函數(shù)模型和擬合

由2.2節(jié)可得到半變異函數(shù)的分布圖，如實驗部分圖2所示。在此基礎(chǔ)上，進(jìn)行半變異函數(shù)的曲線擬合。該步驟是一個建模的過程，類似于回歸分析，目的是將分布圖擬合成連續(xù)的曲線。常用的變異函數(shù)模型有球面模型、指數(shù)模型、高斯模型以及線性模型等[15，16]，由于變異函數(shù)模型的選擇很大程度上需要經(jīng)過大量試驗，反復(fù)選擇調(diào)整，且只能保證某一特定數(shù)據(jù)條件下的相對最優(yōu)解，本文球面模型為例進(jìn)行描述。

球面模型的公式為：

（3）

其中，表示變異函數(shù)，c0為塊金效應(yīng)常數(shù)，c為偏基臺值，c0+c為基臺值，a為變程。因為當(dāng)h=0和h>a時，為常數(shù)，因此只要分析的情況。

不妨記，b0=c0，b1=3c/2a，b2=-c/2a3，x1=h，x2=h3，，則可將：

變形為線性方程的形式：

（4）

設(shè)有一組數(shù)據(jù)h和，亦即一組x1，x2和y，將其代入式（4），得到線性方程組：

（5）

為方便求解，將該非齊次線性方程組轉(zhuǎn)換為AX=b的矩陣形式，其中Xi1=hi，Xi2=hi3，，i=1，2，…，n。

由此，將式（5）的求解轉(zhuǎn)化為線性規(guī)劃的求解問題，有多種有效的解法[17]，此過程不在本文中贅述。

通過上述步驟的計算，即可得到變異函數(shù)的各項參數(shù)，參數(shù)的確立為下一步Kriging插值計算提供了關(guān)鍵性支撐。有關(guān)上述變異函數(shù)擬合效果的檢驗，在下文的實驗部分進(jìn)行描述。

2.4 Kriging插值的計算

Kriging插值的基本公式為：

（6）

從形式上看，式（6）與反距離權(quán)重插值非常相似，每個待插值點的結(jié)果都是若干個相鄰實測點的加權(quán)平均值。但Kriging方法基于包含自相關(guān)的統(tǒng)計模型，其權(quán)重不僅取決于待插值點與實測點間的距離，還取決于基于實測點的整體空間排列，即實測點在空間分布上的影響。因此，Kriging方法的式（6）中，取決于實測點與待插值點的距離，以及待插值點周圍的實測點之間空間關(guān)系的擬合模型。

Kriging插值計算的流程如下：

a.創(chuàng)建一個二維數(shù)組G[X，Y]，對于任意一個待插值格點g[i，j]，g[i，j]∈G[X，Y]。

b.依次選擇一個g[i，j]，根據(jù)2.1節(jié)所述的分塊搜索策略，確定參與本次插值計算的若干自動站，作為實測點。

c.根據(jù)2.2和2.3節(jié)所述變異函數(shù)模型及函數(shù)參數(shù)計算方法，求解方程組，得到一組最優(yōu)的權(quán)重系數(shù)。

d.將得到的代入式（6），求得待插值點g[i，j]的值。需要說明的是，此處的g[i，j]為一個數(shù)值，而實際應(yīng)用中，由于自動站監(jiān)測數(shù)據(jù)是由多個氣象要素組成，因此，嚴(yán)格意義上說，g[i，j]的求解是對多個氣象要素分別進(jìn)行。

e.重復(fù)步驟b—d，直到所有g(shù)[i，j]的值全部求出。

f.輸出G[X，Y]即為預(yù)期的格點場形式的數(shù)據(jù)。此數(shù)據(jù)可用于天氣數(shù)值預(yù)報模式等業(yè)務(wù)系統(tǒng)中作為初始場資料。

3 實驗結(jié)果與分析（The experimental results and

analysis）

實驗數(shù)據(jù)選用2015年9月12日05時整點加密自動氣象站數(shù)據(jù)，經(jīng)過初步質(zhì)量檢驗，剔除缺測和異常數(shù)值后，該數(shù)據(jù)的有效站點數(shù)為28 343個。

3.1 變異函數(shù)擬合的檢驗

隨機選擇靠近南京地區(qū)的一個網(wǎng)格點，其對應(yīng)地理坐標(biāo)為東經(jīng)118.78，北緯32.05，以該點為圓心（即待插值點），西北方向至正北方向的45夾角為選定的方向，以6km為單位距離，計算此區(qū)域內(nèi)各扇區(qū)中的自動站實測溫度的半變異函數(shù)，結(jié)果如圖2所示。

圖2 實驗半變異函數(shù)圖

Fig.2 A schematic diagram of experimental

semivariogram function

圖2中橫坐標(biāo)為單位距離數(shù)，縱坐標(biāo)為半變異函數(shù)值，從該圖函數(shù)值的分布特征可以看出，使用對數(shù)函數(shù)模型可以得到相對較好的擬合效果，該模型的一般公式為：

（7）

根據(jù)式（7），由當(dāng)前實驗數(shù)據(jù)推算得到，c=0.423，c0=-0.189，擬合出的曲線如圖3所示。

圖3 實驗半變異函數(shù)擬合結(jié)果

Fig.3 A fitting result of experimental

semivariogram function

3.2 計算速度的檢驗

在對本文算法進(jìn)行編碼時，還考慮了提取循環(huán)語句中反復(fù)相同計算的部分，盡可能降低算法的時間復(fù)雜度。例如，將每個格點找鄰近自動站的代碼，改進(jìn)為先計算出各個1×1的單位經(jīng)緯網(wǎng)格周邊所包括的鄰近自動站，將此結(jié)果記錄下來，而當(dāng)每個格點找站時，直接讀取其所在單位經(jīng)緯網(wǎng)格預(yù)先遍歷好的鄰近自動站，從而大幅減少無意義的重復(fù)遍歷和計算。本實驗所使用服務(wù)器的硬件環(huán)境為：Intel Xeon E5-2697 V2 CPU，64GB DDR3內(nèi)存，SAS接口15 000轉(zhuǎn)/分硬盤。通過在程序中加入記時器，得到每秒處理的格點數(shù)量如圖4所示。

圖4 計算速度對比

Fig.4 Computing speed comparison

圖4中，橫坐標(biāo)為計算流逝的時間，單位為秒，縱坐標(biāo)為每一秒時間內(nèi)插值計算出的格點數(shù)。可以看出，優(yōu)化前每秒大約只能完成150個格點的計算，以此速度計算完所有G1個格點至少需要42分鐘，而采用本文優(yōu)化方法后，每秒可完成約3 200個格點的計算，所有格點的插值計算歷時約115秒。圖中還可以看出，在最開始的4秒“優(yōu)化后”算法完成的插值格點數(shù)量為0，這是因為此階段程序在進(jìn)行自動站分塊、單位經(jīng)度網(wǎng)格遍歷自動站等預(yù)處理工作。由此得出優(yōu)化前后計算速度的對比結(jié)果，即優(yōu)化后的Kriging插值速度提升了約20倍。

3.3 插值結(jié)果的檢驗

為了檢驗插值結(jié)果的正確性，本文使用交叉檢驗法[18]，隨機選取若干自動站實測溫度數(shù)據(jù)，分別與這些自動站最鄰近的格點插值溫度進(jìn)行比對，結(jié)果如圖5所示。

圖5 站點實測值與格點插值結(jié)果比對

Fig.5 The AWS measured values with the

interpolation results

圖5中，深色線條是自動站實測溫度，淺色線條是格點插值出的溫度。橫坐標(biāo)標(biāo)注的是各個自動站的站名，縱坐標(biāo)是溫度值，單位為攝氏度。為更加直觀地反映Kinging與其他插值方法的優(yōu)劣，使用反距離權(quán)重算法對樣本數(shù)據(jù)進(jìn)行插值，再對上述相同站點作交叉檢驗，這兩種插值算法的輸出結(jié)果與自動站實測溫度值的絕對誤差如圖6所示。

圖6 兩種插值結(jié)果與實測值的絕對誤差比對

Fig.6 Absolute errors of two kinds of interpolation

results and the measured values

從圖6可以看出，Kriging插值的效果普遍好于反距離權(quán)重法，Kriging插值結(jié)果與實測溫度值的絕對誤差一般在0—0.45℃，該波動主要有兩方面原因。一是由于各個格點與比對的自動站之間距離不同，通常情況下，兩者距離越近，插值結(jié)果越接近實測值，反之亦然。二是由于各個自動站的海拔高度不同，一般認(rèn)為在大氣非逆溫的情況下，近地面每升高100米，溫度下降0.65℃，如在青藏高原地區(qū)，相鄰的自動站之間可能有數(shù)百米甚至1km的高度差，在不考慮溫度—高度修正的情況下，這對插值的正確性造成了較大的影響。

4 結(jié)論（Conclusion）

本文利用Kriging插值算法在空間相關(guān)性方面的優(yōu)勢，針對自動站氣象要素在空間分布上的規(guī)律特征，提出了一種改進(jìn)的Kriging插值方法，闡述其計算的算法流程。通過實驗測試表明，改進(jìn)后的方法在計算速度上得到約20倍的提升，插值結(jié)果與實況監(jiān)測數(shù)據(jù)對照，其準(zhǔn)確性優(yōu)于反距離權(quán)重插值方法。特別是該方法在速度方面的大幅提升，為自動站數(shù)據(jù)作為精細(xì)化實時初始場資料的一部分應(yīng)用于數(shù)值模式天氣預(yù)報等業(yè)務(wù)系統(tǒng)成為了可能。

參考文獻(xiàn)（References）

[1] 劉彤，閆天池.我國的主要氣象災(zāi)害及其經(jīng)濟(jì)損失[J].自然災(zāi)

害學(xué)報，2011，20（2）：90-91.

[2] 吳洪.臨近預(yù)報系統(tǒng)SWAN培訓(xùn)手冊[M].中國氣象局氣象干

部培訓(xùn)學(xué)院，北京，2012.

[3] 劉高平，葉金印.基于分布式架構(gòu)的縣級氣象業(yè)務(wù)系統(tǒng)設(shè)計與

實現(xiàn)[J].軟件工程師，2015（04）：36-38.

[4] 華韻子，林紅.基于RIA的WebGIS在天氣預(yù)警系統(tǒng)中的應(yīng)用

研究[J].軟件工程師，2015（09）：37-39.

[5] 鄧林華，等.基于雙線性插值的圖像實時消旋系統(tǒng)的設(shè)計與實

現(xiàn)[J].計算機與現(xiàn)代化，2010（06）：64-66.

[6] 段平，等.自適應(yīng)的IDW插值方法及其在氣溫場中的應(yīng)用[J].

地理研究，2014（08）：1417-1426.

[7] Krige，D.G.Some Basic Considerations in the Application of

Geostatistics to the Valuation of Ore in South African Gold

Mines[J].Journal of the South African Institute of Mining and

Metallurgy，1976.

[8] 李迪，等.基于Kriging與IDW融合算法的氧化帶DTM擬合方

法[J].金屬礦山，2013（10）：93-96；100.

[9] 范玉潔，等.降雨資料Kriging與IDW插值對比分析—以漓江

流域為例[J].水文，2014（06）：61-66.

[10] 徐愛萍，胡力，舒紅.空間克里金插值的時空擴展與實現(xiàn)[J].計

算機應(yīng)用，2011（01）：273-276.

[11] 王興，等.基于OpenCL的雷達(dá)外推算法改進(jìn)與優(yōu)化[J].計算

機與現(xiàn)代化，2014（08）：81-86；90.

[12] Ying-jie XIA，Li KUANG，Xiu-mei.Accelerating geospatial

analysis on GPUs using CUDA[J].Journal of Zhejiang

University-Science（Computers & Electronics），2011（12）：990-999.

[13] 周豐，郝澤嘉，郭懷成.香港東部近海水質(zhì)時空分布模式[J].環(huán)

境科學(xué)學(xué)報，2007（09）：1517-1524.

[14] 翟進(jìn)乾.克里金（kriging）插值方法在煤層分布檢測中的應(yīng)用

研究[D].太原理工大學(xué)，2008.

[15] 占文俊，等.一種InSAR大氣相位建模與估計方法[J].地球物

理學(xué)報，2015（07）：2320-2329.

[16] Q.ZHU，H.S.LIN.Comparing Ordinary Kriging and Regression

Kriging for Soil Properties in Contrasting Landscapes[J].

Pedosphere，2010（05）：594-606.

[17] 馮桂蓮.有限元數(shù)值解法在MATLAB中的實現(xiàn)及可視化[J].

軟件工程師，2015，18（01）：61-64.

[18] Yongping Zhao，Kangkang Wang.Fast cross validation for

regularized extreme learning machine[J].Journal of Systems

Engineering and Electronics，2014（05）：895-900.

作者簡介：

王 興（1983-），男，博士生，講師.研究領(lǐng)域：氣象信息技術(shù)

與安全.

張琳焓（1993-），女，本科生.研究領(lǐng)域：大氣科學(xué).

王麗娟（1983-），女，碩士，講師.研究領(lǐng)域：短期天氣預(yù)報.

張瀟瀟（1987-），女，碩士，講師.研究領(lǐng)域：氣象信息技術(shù)與

教育.