基于矩陣分解的降水時(shí)空插值方法

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

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

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

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

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

1 插 值 方 法

1.1 Funk-SVD模型

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

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

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

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

(1)

1.2 插值方法計(jì)算流程

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

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

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

(2)

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

(3)

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

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

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

(4)

(5)

(6)

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

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

1.3 評(píng)價(jià)指標(biāo)

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

(7)

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

(8)

(9)

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

2 實(shí) 例 驗(yàn) 證

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

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

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

圖4 8場(chǎng)降水事件的降水量空間分布(單位：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)行插值時(shí),需要m和N兩個(gè)參數(shù)來確定矩陣的規(guī)模，其中m是插值所需的站點(diǎn)數(shù)，和矩陣的行數(shù)有關(guān)；結(jié)合N和插值時(shí)刻可以確定矩陣的列數(shù)。為了得到最佳的插值效果，對(duì)這兩個(gè)參數(shù)選取了多種值進(jìn)行演算和比較，結(jié)果見表2。可以看出，當(dāng)m=7時(shí)，L達(dá)到最大值；當(dāng)N=6時(shí)，RSME達(dá)到最小值，所以本研究中m取值為7，N取值為6。

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

2.3 整體精度評(píng)價(jià)

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

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

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

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

2.4 代表性事件的站點(diǎn)誤差比較

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

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

3 結(jié) 語

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