兩種面向數值同化的風廓線雷達資料質量控制方法比較分析

汪學淵 林銀杰 劉德強 林立崢

(1 廈門市氣象局 海峽氣象開放實驗室，福建 廈門 361012；2 福建省災害天氣重點實驗室，福州 350001;3 福建省大氣探測技術保障中心，福州 350001；4 福建省南平市氣象局，福建 南平 353000；5 福建省氣象臺，福州 350001)

引 言

風廓線雷達是利用大氣湍流對電磁波的散射作用進行探測的遙感設備，可以提供探測高度范圍內的大氣水平風速、風向、垂直氣流、大氣折射率結構常數等氣象要素的觀測，具有較高的時空分辨率，彌補了常規探空觀測時空密度不足的缺陷。目前，風廓線雷達數據在監測預警、預報和數值同化中應用較為廣泛，且取得了豐碩的成果。美國和日本的業務應用表明：風廓線雷達資料的同化對于數值模式0～12 h，尤其是3～6 h的預報具有正效果[1-2]；北京、廣東等地都初步開展了一些同化應用的個例試驗，結果表明：在同化了經過質量控制處理的風廓線資料后，區域模式的預報效果取得了顯著的改善，其中完善的質量控制流程則是資料得到有效同化應用的關鍵[3-5]。

近年來，中國氣象局氣象探測中心建立了完善的風廓線雷達資料質量控制和評估業務，分為臺站級和國家級質控體系，臺站級主要對功率譜資料進行質控，國家級主要對徑向數據質控，為風廓線雷達資料的同化應用奠定了基礎。采用變分方法進行資料同化時，觀測誤差和模式背景誤差都必須要滿足高斯分布的假設[4]。因此，在同化應用之前，必須識別和消除觀測數據中不可靠或包含不能滿足數據同化要求的離群值，確保觀測場與背景場的差值(觀測增量)近似與高斯分布相一致。

大氣中的各氣象要素基本上都是一維觀測向量，目前針對單一要素(如溫度、濕度等)的質量控制普遍采用了雙權重標準差(Biweight Standard Deviation, BSD)方法，它通過給定的閾值來剔除離群值，質控效果較好[6-8]。然而，對于水平風場(u/v)而言，BSD方法無法實現對二維觀測向量的同時質控。迭代加權最小協方差行列式(the Iterated Reweighted Minimum Covariance Determinant,IRMCD)[9]是在最小協方差行列式(MCD)[10-11]基礎上發展起來的方法。MCD是應用穩健統計中最早的仿射同變和高魯棒性多元離群點檢測規則之一。自從引入計算效率較快的fast-MCD算法以來[12]，MCD已被應用于醫學，金融，圖像分析和化學等領域。然而，由于傳統MCD方法在檢測離群值時存在一定量的誤判，Cerioli[9]在其基礎上引入了防“假陽性”機制以減少誤判，應用于多元變量離群點檢測。IRMCD可以對多維向量同時進行處理，ZHANG, et al[13]將IRMCD方法用于風廓線雷達水平風離群值檢測發現：IRMCD對于二維風廓線雷達水平風觀測資料的質控效果要好于BSD方法。研究從實際應用角度加深了對這兩種質控方法的認識。然而，由于IRMCD依賴于形狀分布參數，這些參數隨數據集的大小而變化，ZHANG, et al[13]沒有就這些參數對于質控效果的影響進行深入討論。此外也沒有給出晴雨條件下兩種方法質控效果的對比研究。

為了進一步全面深入考察兩種方法的差異性，本文將從統計指標、波形指標、概率密度分布、離群值分布多方面對IRMCD和BSD方法處理風廓線雷達資料離群值的能力和效果進行更深入的對比分析，揭示兩種方法的差異性和優異性。

1 資料和方法

1.1 資料

風廓線資料挑選了福建省運行比較可靠的9部CFL-06型號的雷達資料，分別是：建甌(58737)、建寧(58822)、羅源(58845)、連城(58912)、武平(58917)、德化(58935)、秀嶼(58938)、平和(59125)和翔安(59140)。由于本文的重點在于考察IRMCD方法與BSD方法在混合雷達站點資料處理離群值過程中的性能和效果，所以將生成的風場小時數據作為原始觀測數據。前期關于臺站級和國家級質量控制有關部門和學者已做了大量研究，并取得了積極的研究成果，不再贅述。

利用9部風廓線雷達2018年2月2—11日10 d的小時風場數據作為原始觀測數據，將觀測數據分為降水和非降水天氣，在這里降水和非降水的判定準則按照風廓線雷達垂直速度w≥2 m·s-1判定為降水，獲得了65 000個非降水觀測數據并在其中隨機抽取5 000、10 000、30 000、60 000個觀測數據；同時也獲得了12 750個降水觀測數據并在其中抽取5 000、12 750個觀測數據，以考察IRMCD方法和BSD方法處理不同天氣情況下不同觀測樣本量在統計指標和波形指標上是否有較大差異。

模式背景場數據選取了歐洲數值預報中心(ECWMF)哥白尼CS35數據庫中高空u/v分量的小時再分析數據，并對模式背景數據在垂直和水平方向進行了插值處理，以獲得與觀測數據相同高度的背景場u/v分量，因此，u/v分量觀測增量可以定義為:

ombu(i)=obsu(i)-mu(i)，

(1)

ombv(i)=obsv(i)-mv(i)，

(2)

其中：i=1,2,....n，n表示風觀測數據總量；u，v分別表示風在水平方向兩個分量。ombu(i)表示u分量的觀測增量；ombv(i)表示v分量的觀測增量；obsu(i)表示u分量的觀測值，由OOBS產品文件中的風速V和風向θ根據-V×sinθ計算公式獲得；obsv(i)表示v分量的觀測值，由OOBS產品文件中的風速V和風向θ根據-V×cosθ計算公式獲得；mu(i)表示u分量的模式背景值，mv(i)表示v分量的模式背景值。以下所有指標和參數的計算都是基于u/v分量的觀測增量進行運算，如果觀測增量判定為離群值，那么對應的原始觀測數據定義為離群值。

1.2 迭代權重的最小協方差矩陣方法(IRMCD)

假設n個樣本p個維度的數據集可以表示為：

Y=[y(1)......y(n)]T，

(3)

那么y(i)=(yi1......yip)T為第i個樣本點，矩陣Y的平均值μ和協方差矩陣∑，如果Y中存在離群值，那么μ和∑已經被離群值污染。本文應用穩健統計分析方法，通過檢測每個觀測值魯棒距離的平方與χp,1-α分布相差較大的距離定義為Y中的離群值，可以得到μ和∑的穩健估計值。其中1-α為χ分布的分位數，α一般取0.025。IRMCD是一種基于重加權MCD估計值而發展起來的穩健估計方法[14-15]。對于有限樣本離群值檢測的IRMCD方法的步驟如下：

(1)在樣本Y中，如果h(n/2≤h

(4)

協方差估計為：

，(5)

其中：C0為比例常數[9]。

(2)在Y中,y(i)的魯棒距離的平方可以定義為：

，(6)

它測量了觀測值到假定非離群值的中心位置的距離。樣本Y中所有觀測值的權重系數可以通過DIS的值確定：

(7)

(3)為了增強效率，對y(i)進行加權步驟：

(8)

[y(i)-μRMCD]T，

(9)

那么重新加權后魯棒距離的平方為：

(10)

(4)參考文獻[9]中，

(12)

那么數據集Y中沒有離群值。

按照上述步驟，使用預設的γ值，可以檢測多變量數據集Y中的離群值。

1.3 雙權重標準差方法(BSD)

雙權重離群值判別計算方法(簡稱雙權重標準法，又稱 Z-Score 法)如下:設有n個樣本(xi,i=1,2,...n)

(1)計算每個樣本量xi(i=1,2,..,n)的權重函數:

(13)

其中：C為“敏感參數”，取C=7.5，當|wi|>1.0時，設定wi為1，M為樣本量的中位數，MAD為絕對偏差中位數，即|xi-M|的中位數。

(14)

計算雙權重標準差(BSD)：

(15)

對每一個xi計算Z-Score值：

(16)

如果Zi>Zthresh，那么xi被認定為離群值[16]，Zthresh為設定好的閾值，一般取2～4。

2 結果分析

2.1 基于正態波形指標的最優參數判定準則和指標分析

這里引入了峰度和偏度兩個統計指標來形容觀測增量數據的波形是否符合正態分布情況，峰度(Kurtosis)是描述總體中所有取值分布形態陡緩程度的統計量，峰度為0表示該總體數據分布與正態分布的陡緩程度相同；偏度(Skewness)是統計數據分布偏斜方向和程度的度量，當偏度接近0則可認為分布對稱。兩個指標都是以接近0值為最優值，因此可以組合峰偏值KS指標，表示如下：

KS=|Ku|+|Kv|+|Su|+|Sv|，

(17)

其中：Ku，Kv表示u,v分量的峰度；Su,Sv表示u,v分量的偏度。

那么當IRMCD和BSD方法分別取不同的參數γ和Zthresh時，質控后的觀測增量的KS值應該具有最小值，KS取最小值所對應的參數γ和Zthresh值就是兩種方法的最優解，就是本文所需要的最優觀測增量數據。在以往的研究中，γ參數的典型取值為0.025[17-18]，表示在樣本集中期望2.5%比例的離群值，本文設定γ范圍為0.080～0.001[13]，每0.001的間隔考察KS值是否達到最小值，KS最小值所對應的γ值就是IRMCD處理此次觀測樣本增量的最優解；以同樣的方式對Zthresh的取值范圍設定在4.0～1.0,每0.01的間隔考察KS值是否達到最小值，KS最小值所對應的Zthresh值就是BSD處理此次觀測樣本增量的最優解。從總樣本中隨機抽取了無降水樣本60 000個和降水樣本12 000個，分別繪制了KS值隨γ參數和Z閾值變化曲線(圖1、2)，無降水樣本用藍色表示，降水用紅色表示，γ參數以0.001的間隔在0.080～0.001取值對應一個KS值，從圖1中可以看出，KS值的變化曲線呈現不規則拋物線形狀，有且僅有一個最低點，所對應γ參數就是IRMCD方法所需的最優解，當然對于不同的數據集KS最小值以及γ參數都會有所不同；同樣，Zthresh以0.1的間隔在4.0～1.0取值對應一個KS值，從圖2中可以看出，KS值的變化曲線同樣呈現不規則拋物線形狀，總能找到KS最小值，所對應Zthresh就是BSD方法所需的最優解。這說明所制定的通過峰偏值KS指標判定數據達到最優正態分布的合理性。

圖1 KS值隨γ參數變化曲線

圖2 KS值隨Z閾值變化曲線

從樣本數據中隨機抽取5 000、10 000、30 000、60 000個非降水觀測數據和5 000和12 750個降水觀測數據,分別利用IRMCD和BSD兩種方法通過調整γ和Zthresh使KS值達到最小值，各個參數值如表1所示，其中Ku表示原始觀測u分量增量數據峰度指標，Ku′表示經過IRMCD或BSD方法質控后的u分量增量數據峰度指標，以此類推。從峰度和偏度指標來看，在非降水樣本中u分量的峰度Ku值保持在7.2左右，經過質控后Ku′下降到0.01左右，v分量的峰度Kv值保持在25左右，經過質控后Kv′下降到0.15左右；u分量的偏度Su值保持在-1.6左右，經過質控后Su′下降到0.1左右，v分量的偏度Sv值保持在-4.3左右，經過質控后Sv′下降到0.2左右。從波形指標上看，兩種方法都起到很好的質控效果，在降水天氣下峰度和偏度指標有著類似的趨勢。但是從KS指標和離群值的數量來看，IRMCD始終比BSD方法的質控效果更好。圖3展示了KS指標在不同樣本下的變化曲線，IRMCD方法始終在0.4左右，而BSD方法始終在0.5左右，兩者之間相差0.1，說明IRMCD方法質控后的數據更符合高斯或正態分布；從離群值的數量上來看，IRMCD方法始終比BSD方法判斷的離群值要多，由表2可見，兩種方法能夠判別離群值占總樣本的比例在11%～13%之間，但前者比后者要多0.6%,Avgu和Sdu分別代表u分量的絕對平均值和標準差，以此類推，經過兩種方法的處理后，相對于原始數據都有極大的改進，質控后的Sdu基本保持在2.1～2.3，總體上IRMCD在絕對平均值和標準差指標都優于BSD方法。說明IRMCD方法無論在波形指標、統計指標和離群值數量上都優異于BSD方法，而且兩種方法在樣本的數量多少以及是否降水天氣都不影響各自離群值判斷能力。

表2 IRMCD和BSD不同采樣統計指標表

圖3 兩種方法的KS指標對比曲線

2.2 兩種方法的概率密度和散點分布差異

一般來說，IRMCD和BSD方法在判定離群值的本質上是等價的：給定一個穩健的均值和標準差，數據集向量Y中的離群值通過它們與穩健擬合存在較大距離來識別。以非降水天氣下60 000樣本為例，圖4、5分別為u/v分量的觀測增量在不同方法處理后的概率密度和分位數—分位數(Q-Q)圖，其中U-質控前表示u分量原始觀測增量；U-IRMCD表示u分量觀測增量經過IRMCD質控后的觀測增量；U-BSD表示u分量觀測增量經過BSD質控后的觀測增量，以此類推。這能反映觀測增量數據的分布情況，U-質控前和V-質控前的概率密度分布類似于高斯分布，但不是嚴格的高斯分布，可以看出陡峭的峰值和左右兩側分布的不對稱存在異常值。更準確地說，在相對應Q-Q散射的兩端存在較大差異，與其相對應的u/v觀測增量的峰度值分別為7.35/25.09以及偏度值分別為-1.62/-4.29都說明原始觀測增量數據分布嚴重偏離正態分布。從U-IRMCD和V-IRMCD的概率密度分布和Q-Q散點可以看出質控后的概率密度分布更接近于標準正態分布，Q-Q散點幾乎以直線收斂，表明幾乎所有離群點已被剔除，從相對應u/v觀測增量的峰度值分別為0.0/-0.13以及偏度值分別為-0.07/-0.17，從數值上也說明質控后的數據逼近標準正態分布。同樣的U-BSD和V-BSD的概率密度分布和Q-Q散點以及相對應的峰度值分別為-0.01/-0.21以及偏度值分別為-0.09/-0.20能得到相同的結論，說明兩種方法在剔除離群值后都具有較好的正態分布，但是從峰度值、偏度值、峰偏值和標準差的指標對比來看，明顯IRMCD方法的指標優于BSD方法，從概率密度直方圖的底部兩側還是能看出IRMCD比BSD來得更加平緩；Q-Q散點兩側IRMCD比BSD更加靠近中線位置。值得注意的是，表1的兩種方法的v分量偏度值始終保持在0.2左右，仍然需要最后的偏倚校正[13]。

表1 IRMCD和BSD不同采樣數量波形指標表

圖4 u增量概率密度直方圖和相對應的Q-Q分布

圖5 v增量概率密度直方和相對應的Q-Q分布

圖6、7為u/v分量離群和非離群值散點分布，將進一步理清兩種方法的差異之處。其中“+”表示非離群值，“.”表示離群值，并以不同的顏色代表觀測值所在的高度，為了更加清晰地表示離群值和非離群值，在4～7 km的非離群值用藍色表示，4～7 km的離群值青藍色表示，可以看出0～2 km的離群值以綠色實心圓分布，表明v分量的觀測值大于模式值，4～7 km的離群值以青藍色實心圓分布，表明v分量的觀測值小于模式值為主，在所有的離群值中4～7 km占據了一半以上，這是因為2月的溫度與濕度低造成風廓線雷達的有效探測高度在6 km以下，在有效探測高度以上信噪比越來越弱，生成的風場可靠性降低，造成大量的離群值，同時也可以看到7 km以上存在很少的離群值，因為2月探測高度很少能達到7 km以上。從整體上來看，很明顯，IRMCD和BSD兩者最大的不同在非離群值聚集的形狀上，BSD的非離群值更趨向于“方形”，而IRMCD的非離群值更趨向于“橢圓形”，這是由各自的算法所決定，BSD方法只能處理單向量，根據觀測點偏離標準差的倍數來決定是否為離群值，而IRMCD方法能同時處理二維向量，通過二維向量距離最小協方差矩陣中心的距離是否滿足特定分布來判定是否為離群值，這也是IRMCD方法的優勢所在。

圖6 u/v增量BSD離群和非離群值散點

圖7 u/v增量IRMCD離群和非離群值散點

為了更進一步地理清兩種方法在判定離群值的不同之處，將兩種方法進行對比(圖8)，在非降水情況下兩種方法都判定為離群值用紅色表示，都判定為非離群值用藍色表示，僅僅IRMCD方法為離群值但BSD方法為非離群值用綠色表示，僅僅BSD方法為離群值但IRMCD方法為非離群值用黑色表示，可以看出，紅色點離群值所占比例為11.55%，IRMCD方法判定的離群值所占比例為12.41%，BSD方法判定的離群值所占比例為12.24%，因此大部分離群值兩種方法都能識別，不同的是僅IRMCD方法的非離群值分布更趨向于0值軸附近，在圖8中用綠色部分表示，僅BSD方法判定非離群值分布更趨向于“方形”對角線附近，在圖8中用黑色部分表示，明顯看出黑色點在4個角處且必然存在著離群值，但是BSD方法并沒有識別出來，造成對非離群值的污染，而IRMCD方法識別的非離群值顯得更加的平滑，雖然IRMCD方法也有存在錯誤識別離群值的可能性，但是相對于離群值來說小得多，幾乎可以忽略不計。在降水情況下，如圖9所示，展示了如上所述相近的分布，僅僅BSD識別出的離群值聚集在“方形”的對角線附近，而僅僅IRMCD識別的離群值聚集在0值軸附近。

圖8 u/v增量無降水BSD和IRMCD散點

圖9 u/v增量降水BSD和IRMCD散點

2.3 IRMCD方法質控前后風場變化

為了更好地展示原始觀測風場和質控后數據的變化，圖10、11分別用風羽圖展示了雷達站點(58944)的風廓線，2018年2月8日08時(北京時，下同)至9日14時共計30 h的原始風場和IRMCD質控后的小時水平風廓線，對比發現，原始數據最大探測高度在7 200 m，質控后探測高度在5 000 m，圖11風場廓線顯示明顯比圖10干凈、整潔、有規律，可見離群值主要分布在高空(4.5～7.5 km)和低空(0～0.5 km)，原因是風廓線雷達在4.5 km以上接收到的回波信號很弱，幾乎淹沒在噪聲信號中，造成功率譜信號識別錯誤，就會生成錯誤的水平風；同時由于風廓線雷達低空接收到的回波信號容易受地物雜波的干擾，這些在零頻位置很強地物信號完全將大氣湍流回波信號淹沒，因此生成的水平風風速很小，方向雜亂沒有規律。從圖11中可以看出，IRMCD方法剔除離群值的能力優異，這里不再展示BSD方法處理后的廓線，因為處理后幾乎與圖11一樣，在這么小的樣本情況下幾乎只有2～3個點的區別，這也能從前面表1的指標也能看出。

圖10 2018年2月8—9日風廓線原始小時水平風廓線

圖11 2018年2月8—9日IRMCD質控后小時水平風廓線

因此，這兩種方法在3個方面有所不同：

(1)在雙權重標準差檢查中，Y必須是單變量數據集。當應用于多變量觀測(如風數據)時，需要分別對u/v分量進行異常值檢查，當其中一個向量被認定為離群值，則該樣本二維向量被處理為離群值；另一方面，IRMCD作為一種多變量離群點檢測方法，可以直接應用于多變量數據集Y，即可以同時檢測u/v分量的離群點，在用于風廓線雷達小時觀測增量數據后，從波形指標、統計指標和離群值數量上都表明IRMCD更有效。

(2)它們的穩健均值和標準差是以不同的方式計算的，它們的識別規則也是如此。在IRMCD中，通過比較穩健距離的平方與具有形狀參數分布的參考值進行比較，這些參數隨著應用IRMCD的不同數據集而變化，獲得的非離群值的分布近似“橢圓形”。在雙權重檢查中，設定距離雙權重標準差的預定倍數作為識別離群值的閾值，獲得的非離群值分布近似“方形”，這其中必然存在一定量的誤判，也表明IRMCD比BSD方法有優勢。

(3)IRMCD具有防止假陽性的機制。在IRMCD中，測試的第四步(公式12)是專門設計來防止在任何好的數據集中出現錯誤判定離群值情況[9]，因為誤報是傳統MCD規則的明顯缺點。在沒有步驟4的情況下，IRMCD相當于正常的有限樣本重加權MCD，直接執行第五步會導致錯誤地識別正確的數據集，因此，傳統MCD和雙權重標準差都存在著同樣的缺陷。即使對于一個完美的數據集，離群值也或多或少被錯誤地檢測到。這一點在ZHANG,et al[13]中已經有所驗證，但是在本次樣本執行同樣的過程發現，利用兩種方法都能識別出的非離群值進行試驗發現兩種方法都不能再識別出額外的離群值，因此，并不能完全通過這種方式來說明IRMCD方法比BSD方法更有效果，對于不同的數據集可能會呈現不同效果。

3 結論

本文選取了2018年2月2—11日福建9部風廓線雷達的小時水平風觀測數據與相應的模式數據之差，即觀測增量，利用IRMCD和BSD兩種方法分別進行質量控制，并對質量控制結果以不同的形式進行比較分析。主要總結如下：

(1)制定了IRMCD和BSD質控方法獲得最優解的判定指標峰偏值KS，同時通過KS指標的大小判斷兩種方法的優劣性，IRMCD的KS指標明顯小于BSD方法的KS指標，說明IRMCD比BSD方法更接近正態分布。

(2)IRMCD方法可以同時應用在多維變量的離群值檢測，而BSD方法只能應用在一維變量的離群值檢測中，BSD應用在二維變量離群值檢測的時候必須分別進行離群值檢測，對于具有相關性的兩個變量是不利的。從波形指標、統計指標和離群值數量上都說明IRMCD比BSD更有優越。

(3)IRMCD和BSD的穩健均值和標準差是以不同的方式計算的，它們的識別規則也是如此。在IRMCD中，通過比較穩健距離的平方與具有形狀參數分布的參考值進行比較，這些參數隨著應用IRMCD的不同數據集而變化，獲得的非離群值的分布近似“橢圓形”。在雙權重檢查中，設定距離雙權重標準差的預定倍數作為識別離群值的閾值，獲得的非離群值分布近似“方形”，這其中必然存在一定量的誤判，同時IRMCD具有防止假陽性的機制，這也減少了離群值的誤判，也表明IRMCD比BSD方法有優勢。

從多個方面都表明了IRMCD的在風廓線數據質量控制的優勢，特別是對于二維向量離群值檢測具有普遍意義，也可以應用在激光測風雷達、探空雷達、天氣雷達等設備的風場離群值檢測。也將為下一步在同化業務應用中提供了依據，同時今后也將該方法質控后同化應用于福建區域數值預報模式中，是否能改進數值預報效果，也是下一步的工作目標。