基于協(xié)克里金的雨量雷達(dá)融合置信處理方法

陳根華，莫正威

(南昌工程學(xué)院 信息工程學(xué)院，江西 南昌 330099)

降雨量是降雨型災(zāi)害預(yù)警與降雨徑流模型研究的關(guān)鍵信息。目前，較廣泛運用的降雨測量儀器包括雨量雷達(dá)、雨量計和氣象衛(wèi)星等，利用采集的雨量相關(guān)數(shù)據(jù)，結(jié)合數(shù)據(jù)分析能力形成數(shù)字化天氣預(yù)報模型[1]。本文分析雨量雷達(dá)和雨量計兩類儀器在降雨測量中的融合利用。

雨量計提供較為準(zhǔn)確的降雨測量數(shù)據(jù)[2]。但只能提供局部點數(shù)據(jù)測量，區(qū)域雨量的2-D形成需要使用插值技術(shù)。受限于陡峭地形、寬廣水域及布設(shè)成本，高密度雨量計站網(wǎng)難以建立，造成降雨測量的時間與空間分辨率較低，極易導(dǎo)致“測不到”問題[3]，降低了區(qū)域降雨測量性能，在強(qiáng)對流天氣下尤為明顯，觀測時易漏掉強(qiáng)降雨中心。而雨量雷達(dá)依靠其高時空分辨率，可實現(xiàn)區(qū)域雨量的無縫測量，提供區(qū)域降雨空間變異性的準(zhǔn)確描述，在降雨估測中被廣泛應(yīng)用[4]。雨量雷達(dá)采用電磁波反射理論，通過建立反射率因子與降雨量之間的函數(shù)關(guān)系，實現(xiàn)區(qū)域降雨估計，但由于地物回波、超折射回波、零度亮層[5]以及雷達(dá)元器件等影響，在定量降水估計(Quantitative Precipitation Estimation，QPE)方面的準(zhǔn)確度低于雨量計[6]，形成“測不準(zhǔn)”現(xiàn)象[4]。常用的QPE方法有反距離加權(quán)法[7]、雙調(diào)和樣條插值法[8]、基于三角剖分的插值法等[9]，但以雨量計為單一數(shù)據(jù)來源的方法，因雨量計測量范圍小的缺點，難以詳細(xì)反映區(qū)域內(nèi)降雨的變異性，導(dǎo)致區(qū)域估計降雨情況與實際情況產(chǎn)生偏差，甚至漏掉強(qiáng)降雨中心等。

本文結(jié)合數(shù)據(jù)融合技術(shù)，提出了雨量雷達(dá)數(shù)據(jù)融合置信處理算法。該算法先對雨量雷達(dá)獲取的反射率因子(Radar Reflectivity Factor)進(jìn)行Z-score標(biāo)準(zhǔn)化處理，再以標(biāo)準(zhǔn)化處理后的反射率因子數(shù)據(jù)為輔變量，雨量計數(shù)據(jù)為主變量，利用空間變異理論[10]中的協(xié)克里金法[11]，實現(xiàn)雨量雷達(dá)與雨量計數(shù)據(jù)的融合置信處理。實現(xiàn)了異質(zhì)數(shù)據(jù)的融合，結(jié)合兩種儀器的良好品質(zhì)，較好地解決了雨量計“測不到”和雨量雷達(dá)“測不準(zhǔn)”問題，改善了常規(guī)雨量計站網(wǎng)的區(qū)域降雨觀測能力，提高了雨量雷達(dá)定量降水估計的準(zhǔn)確度。本文同時定義了雨量雷達(dá)降雨估計的相對置信度，并分析了融合處理后雨量雷達(dá)降雨估計的相對置信度分布情況。仿真試驗結(jié)果驗證了雨量雷達(dá)數(shù)據(jù)的融合置信處理算法的正確性與有效性，提高了雨量雷達(dá)QPE的準(zhǔn)確度及測量的穩(wěn)健性。

1 數(shù)據(jù)模型

雨量雷達(dá)在區(qū)域降水探測中，具有探測范圍廣、時空分辨率高等優(yōu)勢，但探測環(huán)境的復(fù)雜多變常影響雨量雷達(dá)的探測精度，導(dǎo)致俗稱的“測不準(zhǔn)”現(xiàn)象[4]。氣象部門常利用高精度的雨量計作為真實降水量的度量儀器，但測量范圍受成本、布設(shè)場地等多種因素的制約無法實現(xiàn)高密度測量，出現(xiàn)所謂的“測不到”問題，如圖1所示。

圖1 雨量雷達(dá)工作示意圖

雨量雷達(dá)通過接收水滴對入射波的反射回波，可計算出散射截面σb，即

σb≈(π5/λ4)|Kω|2D6，

(1)

其中λ是入射波波長；D是水滴直徑，且D≤λ/16；Kω=(m2-1)/(m2+2)，m=n-jnκ，m是水的復(fù)折射率，n是折射率，κ是衰減指數(shù)。

利用雷達(dá)氣象方程

(2)

反射率因子Z為

(3)

其中N(D)表示單位體積內(nèi)直徑范圍在D和D+dD之間水滴粒子的粒子數(shù)；Z的單位為mm6/m3。

由于雷達(dá)反射率因子Z和降雨量R是雨滴尺寸分布(DSD)的兩種不同矩(Moment)形式，因此，Z-R關(guān)系可表示為

Z=aRb

(4)

其中系數(shù)a的變化范圍為16～1200，系數(shù)b的變化范圍為1～2.8[5]。雷達(dá)降水估計的影響因素復(fù)雜，包括異常雜波、零度亮層、衰減、風(fēng)場、溫度、濕度等，導(dǎo)致Z-R關(guān)系在定量降雨估計方面精度不足。因此，本文提出一種融合置信處理算法，結(jié)合雨量雷達(dá)和雨量計兩種儀器的優(yōu)勢，提高雨量雷達(dá)定量降雨估計的準(zhǔn)確度。

工程上，雨量雷達(dá)數(shù)據(jù)是以最小分辨單元為基本單位輸出雨量數(shù)據(jù)，包含距離、方位角、反射率、多普勒譜寬等多種探測信息。本文針對反射率因子進(jìn)行融合處理，因此以雷達(dá)為參考中心，以最小分辨單元為基本單位提取反射率因子和位置數(shù)據(jù)，包括距離r、方位角φ和反射率因子Z，可將雷達(dá)數(shù)據(jù)表示成矢量形式，即D=[r，φ，Z]T，則可得雨量雷達(dá)數(shù)據(jù)集R∈3×N為

R=[D1，D2，…，DN]，

(5)

其中N為雷達(dá)分辨單元個數(shù)。

同理，雨量計數(shù)據(jù)包含距離r、方位角φ和降雨量g，表示成矢量形式，即S=[r，φ，g]T，則雨量計數(shù)據(jù)集G∈3×N為

G=[S1，S2，…，SM]，

(6)

其中M為雨量計站點個數(shù)。

顯然，雨量雷達(dá)分辨單元數(shù)N遠(yuǎn)大于雨量計站點數(shù)M，即N?M。在氣象、水文、地質(zhì)災(zāi)害監(jiān)測等應(yīng)用中，如何利用少量雨量計提高雨量雷達(dá)降雨估計精度具有重要的工程意義。本文擬采用數(shù)據(jù)融合技術(shù)，充分利用雨量雷達(dá)和雨量計的良好品質(zhì)，實現(xiàn)雨雷達(dá)和雨量計在數(shù)據(jù)層[6]的融合，以提高雨量雷達(dá)定量降雨估計的精度。

因此，本文提出融合置信處理算子F (·)，對雨量雷達(dá)和雨量計數(shù)據(jù)進(jìn)行融合置信處理，以提高雨量雷達(dá)降雨估計精度，得到區(qū)域雨量雷達(dá)融合后的降雨數(shù)據(jù)Rf，即

Rf=F(R，G).

(7)

2 雨量雷達(dá)融合置信處理算法

在工程應(yīng)用中，定量降雨估計常以雨量計降雨數(shù)據(jù)為準(zhǔn)，利用二維插值技術(shù)[7-9]估計區(qū)域雨量。以雨量計站點數(shù)據(jù)插值的方法難以反映整體區(qū)域內(nèi)的完整變異性，且在雨量計未覆蓋區(qū)域估計精度下降明顯。雨量雷達(dá)具有探測范圍廣、時空分辨率高等特點，但是受復(fù)雜環(huán)境影響，導(dǎo)致其精度不足。因此，本文擬提出融合置信處理算法，提高雨量雷達(dá)降雨估計精度。本算法在M?N時，亦可實現(xiàn)線性無偏和均方誤差最小[11]，同時對于雨量計布設(shè)密度不均問題，本算法全面分析區(qū)域內(nèi)數(shù)據(jù)的空間相關(guān)性，使得融合后的數(shù)據(jù)更可靠，雨量雷達(dá)測量更穩(wěn)定。

本文的融合置信處理算法，以協(xié)克里金法為基礎(chǔ)[11]，綜合變量的空間連續(xù)和變量間的相關(guān)性，以精度更高的雨量計數(shù)據(jù)為主變量，以雷達(dá)反射率因子為輔變量，采用不同的空間變異函數(shù)，運用協(xié)克里金法實現(xiàn)二者的數(shù)據(jù)融合，以提高雨量雷達(dá)降雨估計的精度。

2.1 標(biāo)準(zhǔn)化處理

由于雨量雷達(dá)反射率因子數(shù)量級相較降雨量數(shù)據(jù)差距過大，直接融合難以體現(xiàn)數(shù)據(jù)優(yōu)勢，本文先利用Z-score標(biāo)準(zhǔn)化方法[13]，實現(xiàn)變異性轉(zhuǎn)換及標(biāo)準(zhǔn)化。Z-score標(biāo)準(zhǔn)化的公式如下：

(8)

(9)

(10)

經(jīng)標(biāo)準(zhǔn)化后得雷達(dá)反射率因子RS為

(11)

其中Z(·)為標(biāo)準(zhǔn)化算子。

2.2 基于協(xié)克里金(CoKriging)算法的融合估計

協(xié)克里金估計算法是一種多元空間估計方法，利用存在區(qū)域協(xié)同關(guān)系(Coregionalization)的主、輔變量間的空間相關(guān)性，建立交叉協(xié)方差函數(shù)(Cross Covariance)和交叉變異函數(shù)(Cross Variogram)，融合分析多變量的空間相關(guān)性和統(tǒng)計相關(guān)性，實現(xiàn)區(qū)域變量的均方誤差最小和線性無偏估計[11]，而雨量雷達(dá)和雨量計在降雨方面的測量，是非完全采樣問題，屬于協(xié)克里金估計中的特例。

因此，本文采用協(xié)克里金估計算法[11]融合雷達(dá)雨量反射率因子及雨量計數(shù)據(jù)，以較高精度的雨量計數(shù)據(jù)G和標(biāo)準(zhǔn)化后的雷達(dá)反射率因子RS分別為主、輔變量，由協(xié)克里金算法得到校正后雨量雷達(dá)數(shù)據(jù)，即

(12)

為求解式(12)中的權(quán)系數(shù)，下面簡要介紹空間變異函數(shù)等基本概念。由空間變異理論[10]可知，主、輔變量的變異函數(shù)γR、γG及交叉變異函數(shù)γRG分別為

(13)

(14)

(15)

式中N′(γ)表示半徑Δr內(nèi)數(shù)據(jù)對的個數(shù)。

由估計子的均方誤差定義可知[14]，協(xié)克里金算法融合后的均方誤差為

(16)

將式(16)展開成協(xié)方差形式，即

(17)

式中CR(·)和CG(·)分別表示雨量雷達(dá)和雨量計數(shù)據(jù)的協(xié)方差函數(shù)；CRG(·)表示交叉協(xié)方差函數(shù)[11]。

由于協(xié)克里金估計算法是無偏且均方誤差最小，因此，求解式(12)中主、輔變量的權(quán)系數(shù)等效為求解式(18)的約束優(yōu)化問題，即

(18)

其中λ=[λ1，λ2，…，λN]T，ξ=[ξ1，ξ2，…，ξM]T。由優(yōu)化理論可知，式(18)可轉(zhuǎn)化成為無約束優(yōu)化問題。由Lagrange乘子法可得：

(19)

式中μ1和μ2為Lagrange系數(shù)。

由無約束優(yōu)化問題的求解方法[14]可知求解目標(biāo)函數(shù)L的零梯度即可，即

?L=0，

(20)

因此，由式(19)和式(20)得到N+M+2階協(xié)克里金估計線性方程組，即

(21)

(22)

由線性代數(shù)理論[14]可將協(xié)克里金估計方程組表示成矩陣形式，即

AX=b，

(23)

其中A為變異矩陣，X為廣義權(quán)系數(shù)矢量，b為廣義交叉變異矢量，且

(24)

X=[λ1，…，λN，ξ1，…，ξN，μ1，μ2]T，

(25)

b=[γGR01，…，γGR0N，γG01，…，γG0M，1，0]T.

(26)

由最小二乘法(LS)[14]可得廣義權(quán)系數(shù)矢量X為

X=(ATA)-1ATb，

(27)

將權(quán)系數(shù)矢量λ和ξ代入式(12)可得融合后雨量數(shù)據(jù)，即

(28)

(29)

2.3 相對置信度(Relative Confidence Degree，RCD)

針對雷達(dá)反射率數(shù)據(jù)動態(tài)變化范圍大的問題，傳統(tǒng)的絕對誤差和相對誤差難以反映雨量雷達(dá)的測雨性能，在雨量雷達(dá)全局誤差分析方面，歸一化的評價方法無法實現(xiàn)有效分析，因此，結(jié)合概率中的置信區(qū)間概念，本文提出了一種新的誤差評判標(biāo)準(zhǔn)，即相對置信度，以實現(xiàn)誤差的全局分析。

本文將相對置信度Cα定義為在給定相對誤差門限α下，所有雨量雷達(dá)測量相對誤差小于α的分辨單元數(shù)Nα與雷達(dá)測量分辨單元總數(shù)N之比，即

(30)

其中α定義為相對誤差，即

(31)

由相對置信度的定義可知，相對置信度可分析不同相對誤差門限下的雨量雷達(dá)校準(zhǔn)誤差的空間分布情況，從而反映校準(zhǔn)后雨量雷達(dá)降雨數(shù)據(jù)的可信度。

綜上所述，本文提出的雨量雷達(dá)融合置信處理算法流程圖如圖2所示。

圖2 雨量雷達(dá)融合置信處理算法流程圖

3 實測與仿真實驗

設(shè)雨量雷達(dá)分辨單元數(shù)為N=71×71，隨機(jī)選取雨量站位置，雨量雷達(dá)數(shù)據(jù)的隨機(jī)誤差滿足高斯分布，仿真實驗中每個數(shù)據(jù)點為100次Monte Carlo實驗。

試驗1分析雨量雷達(dá)和雨量計實測數(shù)據(jù)。圖3(a)為2016年3月12日8時到16時PR-11A雨量雷達(dá)在江西省贛州市信豐縣古陂河流域的雷達(dá)監(jiān)測雨量分布圖，圖3(b)表示雨量雷達(dá)相對雨量計的相對誤差，其最大相對誤差高達(dá)70%，平均誤差約為55%。因此，雨量雷達(dá)的準(zhǔn)確度仍較低，須提高雨量雷達(dá)數(shù)據(jù)有效性。

圖3 PR-11A雨量雷達(dá)實測數(shù)據(jù)分析

圖4 融合置信處理算法分析示意圖

表1 不同降雨量參數(shù)設(shè)置

表2 不同降雨類型Z-R關(guān)系系數(shù)設(shè)置

試驗4分析雨量計站點數(shù)M對融合置信處理算法的影響。仿真條件同試驗3，利用相對置信度分析不同M時融合處理算法的性能。由圖6可知，融合置信處理算法的相對置信度均優(yōu)于其他算法，且對M不敏感，驗證了融合置信處理算法對雨量計站點數(shù)M的穩(wěn)健性。

圖5 不同變異性下相對置信度分析

圖6 不同M時相對置信度分析

4 結(jié)束語

針對雨量雷達(dá)測雨精度不足、雨量計測量范圍有限的問題，本文結(jié)合空間變異理論和數(shù)據(jù)融合技術(shù)，提出了一種中小流域雨量雷達(dá)融合置信處理算法。該算法先對雨量雷達(dá)反射率因子數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，實現(xiàn)變異性轉(zhuǎn)換及標(biāo)準(zhǔn)化，再以標(biāo)準(zhǔn)化后雷達(dá)反射率因子作為輔變量，雨量計數(shù)據(jù)為主變量，利用協(xié)克里金法實現(xiàn)數(shù)據(jù)層融合，進(jìn)行定量降水估計。實測數(shù)據(jù)與仿真結(jié)果表明，融合置信處理算法實現(xiàn)了雨量雷達(dá)與雨量計數(shù)據(jù)層融合，提高了雨量雷達(dá)的定量降雨測量的相對置信度，驗證了本文融合置信處理算法的適用性與穩(wěn)健性。本文為工程上雨量雷達(dá)的定量降水估計提供了理論指導(dǎo)，具有重要的工程意義。未來也將研究天空地多源降雨數(shù)據(jù)的融合技術(shù)，并結(jié)合深度學(xué)習(xí)等方法實現(xiàn)降雨估計及降雨趨勢預(yù)測等。