陳 華，盛 晟，夏潤亮，馬 瑞，朱躍龍，李 濤

(1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072;2.黃河水利委員會黃河水利科學(xué)研究院，河南 鄭州 450003;3.長江空間信息技術(shù)工程有限公司，湖北 武漢 430010；4.河海大學(xué)計算機(jī)與信息學(xué)院，江蘇 南京 210098)

降水?dāng)?shù)據(jù)在國家防汛、抗旱等自然災(zāi)害應(yīng)急決策管理中具有重要的作用，精確的雨量估算和空間分布確定能為自然災(zāi)害的應(yīng)急決策提供科學(xué)準(zhǔn)確的數(shù)據(jù)[1]，提升應(yīng)急管理水平，社會效益顯著。雨量站網(wǎng)是目前最主要的降水觀測和采集系統(tǒng)，其主要特點是便捷、實時和準(zhǔn)確；但由于站點是離散的，需要借助于空間插值方法來得到降水的空間連續(xù)分布。

現(xiàn)有的降水空間插值方法，不論是整體插值(如邊界內(nèi)差法、趨勢面分析法)還是局部插值(如克里金法[2-3]、泰森多邊形法[4]、反距離權(quán)重法[5])，都是基于插值點與雨量站的地理位置關(guān)系，利用插值的當(dāng)前時刻降水信息進(jìn)行插補(bǔ)。這些方法主要是基于“地理學(xué)第一定律”[6]的基本假設(shè)：空間位置上越靠近的點，具有相似特征值的可能性越大；而距離越遠(yuǎn)的點具有相似特征值的可能性越小。這些方法都僅僅考慮了空間上的關(guān)系，沒有考慮降水序列的過去臨近信息對當(dāng)前值的影響，無法將歷史降水信息應(yīng)用于當(dāng)前時刻降水判斷和未來趨勢預(yù)測。也有不少學(xué)者開展了降水時空插值方法研究，比如STARMA和具有時間擴(kuò)展模塊的克里金插值方法[7-9]?，F(xiàn)有的結(jié)果都表明，考慮時間維度的插值方法精度都比只考慮空間維度的插值方法精度高，但由于算法過于復(fù)雜，并沒有被廣泛推廣應(yīng)用。

為了充分利用降水的時間序列和空間信息以提高降水空間分布及雨量估算的精度，本文引進(jìn)了推薦系統(tǒng)中的矩陣分解方法，來實現(xiàn)降水?dāng)?shù)據(jù)的時空插值，并對降水空間插值結(jié)果進(jìn)行了評價。

1 插 值 方 法

1.1 Funk-SVD模型

矩陣分解是解決協(xié)作過濾問題最常用的方法之一，其將用戶的商品偏好程度看作用戶-商品評分矩陣，并使用已知用戶對商品的評分來預(yù)測用戶在商品選擇中的偏好并進(jìn)行推薦[10-11]。由于具有較高的預(yù)測準(zhǔn)確度，矩陣分解已成為環(huán)境科學(xué)、生物醫(yī)學(xué)等許多領(lǐng)域常用的插值方法[12-14]。González-Macías等[15]使用正矩陣分解法來識別和確定海底沉積物中人為重金屬的來源；Lee等[16]將非負(fù)矩陣分解法應(yīng)用于基因數(shù)據(jù)的表達(dá)，量化了不同劑量的pan-PPAR興奮劑在小鼠體內(nèi)引起的分子變化。Funk-SVD模型是Funk[17]提出的一種矩陣分解的變體模型，其基本思想是用戶和商品可以分別由從用戶-商品評分矩陣推斷出的潛在特征向量來描述，由用戶和項目特征之間的高度對應(yīng)性形成預(yù)測和推薦[18]，這種模型也稱作隱語義模型。

基于Funk-SVD模型的降水時空插值模型結(jié)構(gòu)如圖1所示?？臻g上不同點在多個時刻的降水?dāng)?shù)據(jù)可以看作是一個內(nèi)在相關(guān)的矩陣P=(Pij)，其中列代表時間，行代表空間。無降水記錄的點在矩陣P中對應(yīng)位置為空值，即為待插值點。為了得到該位置的值，可以根據(jù)已有的降水記錄對矩陣P進(jìn)行分解，通過分解后兩個矩陣的乘積來補(bǔ)全原矩陣P中空缺位置的值。

圖1 基于Funk-SVD模型的降水時空插值方法結(jié)構(gòu)示意圖Fig.1 Structure diagram of spatiotemporal interpolation method of rainfall based on Funk-SVD model

具體來說，降水?dāng)?shù)據(jù)矩陣的時間和空間之間存在隱式關(guān)系，無法直接求解，所以需要潛在因子來建立間接關(guān)系。假設(shè)存在影響降水的q個潛在因子，從圖1可以看出，該矩陣P可以分解為空間特征矩陣X和時間特征矩陣Y，分別包含m個q維行向量X1、X2、…、Xm和n個q維列向量Y1、Y2、…、Yn，其中Xi(i=1，2，…，m)代表i站點q個潛在因子的影響程度，其包含的特征數(shù)據(jù)越相近，潛在因子對降水的影響效果越相近；Yj(j=1，2，…，n)描述j時刻潛在因子的關(guān)聯(lián)關(guān)系，相關(guān)程度越高，對應(yīng)的特征值越大。因此，可以基于潛在因子的影響程度和關(guān)聯(lián)關(guān)系計算某時某地的降水量。也就是說，Pij可以由Xi和Yj相乘得到，其計算公式如下：

(1)

1.2 插值方法計算流程

基于Funk-SVD模型的降水時空插值方法具體計算流程如圖2所示。圖2涉及的主要有4個計算步驟：

圖2 基于Funk-SVD模型的降水時空插值方法計算流程Fig.2 Calculation flowchart of spatiotemporal interpolation method of rainfall based on Funk-SVD model

a.對于j時刻在流域上的某點I，可以確定一個由m個站點和n個時刻組成的降水?dāng)?shù)據(jù)矩陣。所涉及的m個站點是在計算所有可能排列組合后由空間均勻度L選出的最優(yōu)集合，空間均勻度是點集的空間關(guān)系度量，其值越大，點分布越均勻，計算公式為

(2)

式中：A——涵蓋所有站點的矩形網(wǎng)格面積；a——獨(dú)占圓面積的總和，某點的獨(dú)占圓的定義為以該點為圓心，與最近點距離的一半為半徑的圓?？紤]到矩陣分解的效率和精度，時刻數(shù)n的取值如下：

(3)

式中：t0——降水事件的開始時刻;N——前期影響時刻總數(shù)閾值。如果降水持續(xù)時間不長，n個時刻就從降水開始時刻t0開始取到插值時刻j。如果降水持續(xù)時間長，對于距開始時間超過N個時刻的插值時刻來說，n個時刻就取j時刻以及前N個時刻。

b.增加一行代表流域上一點I，得到一個新矩陣R。矩陣第m+1行前n-1個值為已知的j時刻前的降水，由傳統(tǒng)插值方法計算得到。第n個值為空，代表待插值點。

c.利用Funk-SVD模型將矩陣R分解成時間特征矩陣X和空間特征矩陣Y。Funk-SVD模型通過最小化所有已知點的預(yù)測值平方誤差來計算時間和空間的潛在特征。此外，為了防止出現(xiàn)過度擬合的情況，將正則化方法引入目標(biāo)函數(shù)：

(4)

(5)

(6)

式中：α——機(jī)器學(xué)習(xí)的學(xué)習(xí)率。

d.將優(yōu)化得到的時空特征矩陣相乘，可以得到一個無限逼近于原矩陣的重構(gòu)矩陣R′，且R′中每個位置都有值。原始矩陣R和重構(gòu)矩陣R′中各元素為一一對應(yīng)的關(guān)系，重構(gòu)矩陣R′第m+1行第n列的值即為點I在j時刻的降水量。

1.3 評價指標(biāo)

使用交叉驗證的方法進(jìn)行插值的精度驗證，即每次將一個站點看作未知點，利用其余所有站點進(jìn)行插值，將結(jié)果和觀測值進(jìn)行比較。采用均方根誤差(RMSE)、平均絕對誤差(MAE)、百分誤差(PERC)和雙樣本Kolmogorov-Smirnov檢驗的概率值(KS)[19-20]等4個評價指標(biāo)來驗證插值精度，采用反距離權(quán)重法和普通克里金法這兩種廣泛應(yīng)用的空間插值方法與Funk-SVD模型結(jié)合進(jìn)行插值。這兩種方法都是基于站點的觀測值和權(quán)重來得到待插值點的計算值，計算公式如下：

(7)

式中：Z——待插值點的計算值；Zi——站點的觀測值；wi——站點的權(quán)重，其中IDW和OK兩種方法計算站點權(quán)重的公式分別為

(8)

(9)

式中：di——站點到目標(biāo)點的距離；p——距離的冪；μ——拉格朗日乘子；γ——空間上兩點間的半變異函數(shù)；xi、xh——插值利用到的站點；x——待插值點。

2 實 例 驗 證

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

選擇黃河流域的小浪底到花園口區(qū)間為研究區(qū)域，該區(qū)域位于黃河下游(東經(jīng)109°42′～113°54′、北緯33°17′～37°30′)，是黃河流域主要暴雨來源區(qū)之一，暴雨來勢急、匯流快，短時間內(nèi)可能出現(xiàn)洪水陡漲局面，對黃河防汛安全威脅較大。區(qū)間內(nèi)有伊洛河和沁河兩條國家一級支流，是暴雨洪水多發(fā)地區(qū)。降水量年內(nèi)分布不均勻，集中在6—9月。研究區(qū)有16個雨量站(圖3)，選取2009—2012年不同氣象條件下8次降水事件的日降水?dāng)?shù)據(jù)來驗證本文提出的基于Funk-SVD模型的降水時空插值方法(以下簡稱本文方法)。表 1為8個降水事件的詳細(xì)信息，每場降水的空間分布如圖4所示。這些降水事件持續(xù)時間為6～10 d，分布在6—9月，面降水量跨度很大，從17～67 mm，具有較好的代表性。

圖3 黃河流域小浪底到花園口區(qū)間測站分布Fig.3 Distribution of gauge stations between Xiaolangdi and Huayuankou

圖4 8場降水事件的降水量空間分布(單位：mm)Fig.4 Spatial distribution of rainfall in eight selected events

表1 降水事件信息Table1 Information of selected rain events

2.2 分解矩陣參數(shù)的影響

采用Funk-SVD模型利用矩陣進(jìn)行插值時,需要m和N兩個參數(shù)來確定矩陣的規(guī)模，其中m是插值所需的站點數(shù)，和矩陣的行數(shù)有關(guān)；結(jié)合N和插值時刻可以確定矩陣的列數(shù)。為了得到最佳的插值效果，對這兩個參數(shù)選取了多種值進(jìn)行演算和比較，結(jié)果見表2?？梢钥闯觯?dāng)m=7時，L達(dá)到最大值；當(dāng)N=6時，RSME達(dá)到最小值，所以本研究中m取值為7，N取值為6。

表2 不同m、N值的計算結(jié)果Table 2 Calculation results of different m and N

2.3 整體精度評價

使用4個不同的評價指標(biāo)對各插值方法的計算結(jié)果進(jìn)行了總體評價，結(jié)果如表3所示。從表3可以看出，4種插值方法的精度均在合理范圍內(nèi)，表明它們可以很好地應(yīng)用于研究區(qū)域的降水插值計算。使用不同評價指標(biāo)的評價結(jié)果是一致的，即通過與Funk-SVD結(jié)合，IDW和OK精度均得到了明顯的提升。在RSME上精度提升了9%左右，在其余3個指標(biāo)上的精度提升均在15%以上。除去KS，其他評價指標(biāo)對IDW和OK的精度評價結(jié)果的差異很小，幾乎可以忽略。

表3 8場降水事件的綜合評價結(jié)果Table 3 Comprehensive evaluation results of 8 rainfall events

圖5為4種插值方法在每場降水事件交叉驗證中的所有站點的誤差均值，可以看出，對于每一場降水事件，通過與本文方法結(jié)合，IDW和OK的精度在各項指標(biāo)上均有所改善，這種改善在PERC和KS上更加明顯。RSME和MAE對不同事件的評價結(jié)果是類似的，這是因為它們的數(shù)值與降水量大小有非常大的關(guān)聯(lián)，所以降水量越大，計算得到的誤差就越大。而對個別事件，使用PERC和KS會給出不一樣的評價結(jié)果。比如5號降水事件，它的面總降水量是所有事件中最大的，所以絕對誤差也大，相應(yīng)RSME和MAE也最大；而PERC的數(shù)值大小與相對誤差的關(guān)系更大，受絕對誤差影響較小，因此PERC反而小。

圖5 降水事件插值精度評價結(jié)果Fig.5 Interpolation accuracy evaluation results of each rainfall event

2.4 代表性事件的站點誤差比較

為了更好地評價本文方法的精度，選擇5號大雨事件和7號小雨事件這兩場典型的降水事件進(jìn)行比較分析。圖6和圖7為使用不同插值方法得到的插值與誤差結(jié)果。從圖6和圖7的(b)(c)和(e)(f)可以發(fā)現(xiàn)，無論是大雨還是小雨，降水量越大的點，往往插值產(chǎn)生的誤差越大，并且大雨事件的插值誤差會比小雨事件的大；同時，周圍站點分布不均的點也會產(chǎn)生相對大的誤差。圖6和圖7的(d)(g)中紅色(改善)的點比藍(lán)色(未改善)的點多，說明無論大雨還是小雨，通過與Funk-SVD結(jié)合，可以明顯提升傳統(tǒng)方法插值的精度。

圖6 5號降水事件插值結(jié)果Fig.6 Interpolation results of No.5 rainfall event

3 結(jié) 語

本文基于站點的空間分布關(guān)系，利用當(dāng)前時刻降水信息和歷史降水信息，結(jié)合傳統(tǒng)插值方法，提出了一種基于Funk-SVD模型的降水插值方法。采用黃河流域小浪底到花園口區(qū)間的8場降水事件中16個雨量站的日降水?dāng)?shù)據(jù)，通過交叉驗證的方式對提出的方法進(jìn)行精度檢驗，發(fā)現(xiàn)與Funk-SVD模型結(jié)合可以提升IDW和OK插值的精度，證實了本文方法的可行性。本文方法提高了降水空間估算的精度，為降水時空插值的相關(guān)研究提供了新的思路。