基于天氣雷達反射率特征匹配的降水估測方法
【專利摘要】本發(fā)明公開了基于天氣雷達反射率特征匹配的降水估測方法，首先構建雷達反射率特征庫，然后確立實況雷達數(shù)據(jù)與雷達反射率特征庫的匹配規(guī)則，并在此基礎上動態(tài)確立Z?R關系中最優(yōu)的參數(shù)A和參數(shù)b。本發(fā)明在短臨降水預報業(yè)務應用中實現(xiàn)Z?R關系的動態(tài)最優(yōu)，進一步提高雷達估測降水的準確性。
【專利說明】
基于天氣雷達反射率特征匹配的降水估測方法
技術領域
[0001] 本發(fā)明屬于大氣科學領域，特別涉及了基于天氣雷達反射率特征匹配的降水估測 方法。
【背景技術】
[0002] 天氣雷達（以下簡稱雷達)是強對流天氣探測的重要手段，它能夠探測其掃描半徑 內的雨強和雨量分布情況，進一步可以估測總降水量和降水率，尤其在短時臨近預報等業(yè) 務中已得到廣泛應用。
[0003] 基于Z-R關系的計算是雷達估測降水量的關鍵方法，而估測降水量的準確率很大 程度上依賴于Z-R關系（Z = ARb)中參數(shù)A和參數(shù)b的選擇。很多學者通過對地面實測降水量 與雷達回波的分析研究，得出特定區(qū)域和(或)特定時期內相對最優(yōu)的參數(shù)A和參數(shù)b，由此 作為未來對同區(qū)域、同氣象條件下雷達估測降水Z-R關系計算時的重要參考。
[0004] 在短時臨近預報業(yè)務軟件中，其參數(shù)A和參數(shù)b或是針對當?shù)貧庀髼l件和降水類型 設定為一組相對最優(yōu)值，或是按時期、按降水類型(如大陸強對流降水型、熱帶降水型等)設 定多組相對最優(yōu)值。
[0005] 然而，由于雷達回波的時空變化性很大，從雷達回波圖像上看，幾乎每次降水過程 的回波圖像都是獨一無二的，其Z-R關系的最優(yōu)參數(shù)A和參數(shù)b也并非固定不變的。盡管部分 氣象業(yè)務軟件中Z-R關系的參數(shù)是可調的，但在某一降水事件發(fā)生前，無論人工或是軟件自 動地確定適當?shù)膮?shù)都是非常困難的。因此，如何在短臨降水預報業(yè)務應用中，實現(xiàn)Z-R關 系的動態(tài)最優(yōu)，進一步提高雷達估測降水的準確性需要特別關注。

【發(fā)明內容】

[0006] 為了解決上述【背景技術】提出的技術問題，本發(fā)明旨在提供基于天氣雷達反射率特 征匹配的降水估測方法，通過構建雷達反射率特征庫，以及將反射率特征與雷達反射率特 征庫進行匹配，動態(tài)確立Z-R關系中最優(yōu)的參數(shù)A和參數(shù)b，提高了降水估測的準確性。
[0007] 為了實現(xiàn)上述技術目的，本發(fā)明的技術方案為：
[0008] 基于天氣雷達反射率特征匹配的降水估測方法，包括以下步驟：
[0009] 步驟一、構建天氣雷達反射率特征庫：
[0010] (1)提取天氣雷達基數(shù)據(jù)中的基本反射率數(shù)據(jù)，并將基本反射率數(shù)據(jù)從原來的極 坐標形式轉換為平面直角坐標系形式;分別計算P千米和q千米高度處的等高平面位置顯示 圖?△??1_口和0厶?？1_9,0厶?？1_口和0厶??1_9的結構相同，均為~\袖勺二維網(wǎng)格，求取0厶?？1_口 和CAPPI_q各網(wǎng)格點在同一水平位置的最大回波強度值CAPPI_MAX，即CAPPI_MAX(x，y) = Max(CAPPI_p(x，y)，CAPPI_q(x，y))，其中x、y分別表示平面直角坐標系下的橫坐標和縱坐 標，Max為取最大值函數(shù)；
[0011] (2)根據(jù)得到的最大回波強度值CAPPI_MAX (x，y)，計算雷達回波的灰度圖像Rad (x,y)：
[0012]
[0013] (3)計算灰度圖像Rad(x，y)的一階梯度值，從而強化回波圖像的變化特征；
[0014] (4)新建一個MXM的雷達回波特征矩陣？68，]\1〈1'1/2，？63中每個單元的值計算方式 為：
[0015] Fea(x'，y '）= Hex(Fea'（X，，y '））
[0016]
[0017] 其中，STEP = |三|，1'￡[0，-1]，7'￡[0，-1]此叉為對括號中的值作16進制轉換 Μ 的函數(shù)；
[0018] (5)將得到的雷達回波特征矩陣Fea，從Fea(0,0)到Fea(M-l，M-l)逐個添加到一個 字符串變量FeaStr中；
[0019] (6)選取雷達有效探測范圍內的地面雨量觀測設備，以當前雷達探測時間為基準， 提取這些地面雨量觀測設備在基準后的X小時內的累積觀測降水量，形成數(shù)據(jù)集0RF，其中 〇RF[c]表示第c個地面雨量觀測設備的X小時降水量，xe[l，24]，將當前雷達有效探測范圍 內地面雨量觀測設備的總數(shù)記為MSum，即數(shù)據(jù)集0RF的大小為MSum，則ce [0，MSum-l];
[0020] (7)根據(jù)數(shù)據(jù)集0RF中每個地面雨量觀測設備的經(jīng)煒度信息，以及當前雷達的地理 坐標信息，提取各個地面雨量觀測設備所對應地理位置的CAPPI_MAX值；定義CM為CAPPI_ MAX的一個子集，CM中的每個值CM[ c ]為與0RF[ c ]地理位置最相鄰的網(wǎng)格點的回波強度值； 設置兩個遞增變量PA和Pb，其中PAe [ 100，400]，以10為步長遞增，Pb ￡[1.0,2.0]，以0.1為 步長遞增，分別計算0RF與CM中各記錄項的誤差平方和：
[0021]
[0022] (8)記錄下所得D (PA，Pb)取值最小時的一組PA和Pb，分別記為;
[0023] (9)將步驟(5)中得到的FeaStr和步驟(8)得到的添加到數(shù)據(jù)庫；
[0024] (10)依據(jù)上述步驟（1)-(9)，將有資料以來，雷達探測數(shù)據(jù)和地面雨量觀測數(shù)據(jù)逐 一進行計算并入庫，從而構建起天氣雷達反射率特征庫；
[0025] 步驟二:將實況天氣雷達反射率數(shù)據(jù)與建立的天氣雷達反射率特征庫進行匹配：
[0026] (11)將需要匹配的天氣雷達反射率數(shù)據(jù)按照上述步驟（1)-(5)，計算出能夠反映 該雷達回波特征的字符串FeaStrCur;
[0027] (12)從天氣雷達反射率特征庫中遍歷同一部天氣雷達的FeaStr，形成數(shù)據(jù)集FS， FS中的每條記錄FS[c]與FeaStrCur作相關性計算，計算方法為：
[0028]
[0029] 其中，F(xiàn)eaStrCur和FS[c]分別表示FeaStrCur與FS[c]中所有記錄的算術平均值， FS [ c ] (k)表示FS [ c ]的第k項記錄值；
[0030] (13)找到R(c)取值最大時，F(xiàn)S[c]所對應的數(shù)據(jù)庫記錄，該記錄中對應的 Para_b即為當前天氣雷達資料下Z-R關系的最優(yōu)參數(shù)A和參數(shù)b，并據(jù)此估測降水。
[0031] 進一步地，在步驟(1)中，取p = l .5，q = 3.0。
[0032]進一步地，在步驟(3)中，灰度圖像Rad(x，y)最邊緣一圈的各像素的一階梯度為0， 其余像素的一階梯度采用Sobel檢測算子來計算。
[0033] 進一步地，橫向Sobel檢測算弓
縱向Sobel檢測算子
。
[0034] 進一步地，在步驟(5)中，字符串變量FeaStr中每2個字節(jié)代表1個Fea(x'，y'）。
[0035] 進一步地，在步驟(13)中，若找到的最大R( c)〈0.5時，則Z-R關系中的最有參數(shù)A和 參數(shù)b的取值改用預先設定值。
[0036] 采用上述技術方案帶來的有益效果：
[0037] (1)本發(fā)明對雷達回波進行了一階梯度信息提取，強化雷達回波強弱區(qū)間的邊緣 特征，為不同時刻雷達反射率的比較、匹配提供了分析的數(shù)據(jù)基礎。
[0038] (2)本發(fā)明將雷達回波特征信息以字符串形式加以存儲，便于數(shù)據(jù)庫的實現(xiàn)以及 特征匹配時的快速計算。
[0039] (3)本發(fā)明中基于字符串和16進制數(shù)據(jù)格式的特征值的相關性計算，具有較高的 計算性能，有利于本方法業(yè)務化實施時，提高雷達估測降水預報的時效性。
【附圖說明】
[0040] 圖1是本發(fā)明中構建雷達反射率特征庫的流程圖，對應步驟1-9,圖中實線表示程 序流，虛線表示數(shù)據(jù)流，當程序流與數(shù)據(jù)流重疊時，虛線省略；
[0041] 圖2是本發(fā)明中實況雷達數(shù)據(jù)與雷達反射率特征庫的匹配的流程圖，對應步驟11-13，圖中FS. Count表示數(shù)據(jù)集FS中的總記錄數(shù)。
【具體實施方式】
[0042] 以下將結合附圖，對本發(fā)明的技術方案進行詳細說明。
[0043]本實施例選用2010年至2015年期間，南京地區(qū)多普勒天氣雷達9層體積掃的基本 反射率數(shù)據(jù)。該雷達數(shù)據(jù)格式為CINRAD-SA，體掃模式為VCP21，時間分辨率為6分鐘，徑向分 辨率為lkm，所處地理位置為東經(jīng)118.698°、北煒32.191°。剔除沒有明顯降水回波的雷達數(shù) 據(jù)，并對剩余待計算的雷達數(shù)據(jù)作必要的質量控制。然后，對這些雷達數(shù)據(jù)逐一作如下處 理：
[0044] 步驟1、不妨設當前待處理的雷達基數(shù)據(jù)文件名為RF.bin，提取該文件中的基本反 射率數(shù)據(jù)，將原有極坐標形式的數(shù)據(jù)結構轉換為平面直接坐標形式。分別計算1.5km和 3.01〇11高度的04??1，得到04??1_15和04??1_30。再求取這兩個高度，各網(wǎng)格點在同一水平位 置上的最大回波強度值，即CAPPI_MAX( X，y)=MaX(CAPPI_p(X，y)，CAPPI_q(X，y))，其中 X、y 分別表示平面直角坐標系下的橫坐標和縱坐標，Max表示取最大值的函數(shù)。根據(jù)該部雷達的 性能技術指標，定義N=230,因此xe[0,229]，ye[0,229]。
[0045] 步驟2、由CAPP I_MAX，輸出雷達回波的灰度圖像，記為Rad，Rad中的值滿足：
[0046]
[0047]步驟3、采用Sobel檢測算子計算Rad的一階梯度值，計算方法為：
[0048]
[0049] 其中，*表示卷積計算，展開后的計算方法形如：
[0050]
[0051]
[0052] 步驟4、新建一個MXM的雷達回波特征矩陣Fea，M=46，SETP = 5，F(xiàn)ea中每個單元的 值計算方式為：
[0053] Fea(x'，y '）= Hex(Fea'（X，，y '））
[0054]
[0055] 其中，SETP = 5，x' e[0,45]，y' e[0,45]，Hex表示一個函數(shù)，功能是對其括號中的 值作16進制轉換。
[0056]步驟5、將Fea 中的各個元素，從Fea [ 0，0 ]、Fea [ 0，1 ]、Fea [ 0，2 ]依次到Fea [ 45，45 ] 逐個添加到一個字符串變量FeaStr中。例如，假設?6&[0,0]、？6&[0，1]、？6 &[0,2]至?6&[0， 5]值分別為：0、61、18、42、59、29，則？6&3釷中的記錄形如 ：003012243810。
[0057] 步驟6、選取以該雷達地理坐標(東經(jīng)118.698°北煒32.191°)為中心，230km為半徑 范圍內的地面雨量觀測設備(包括常規(guī)地面觀測站、加密自動氣象站、單要素雨量計），以該 雷達探測時間為基準，提取這些設備在此時間后的1小時內累積觀測降水量，形成數(shù)據(jù)集 ORFARFk]則表示第c個地面雨量觀測設備的1小時雨量觀測記錄。不妨假設當前0RF數(shù)據(jù) 集的記錄總數(shù)為312個，即雷達有效覆蓋區(qū)域內存在312個地面雨量觀測設備。
[0058]步驟7、根據(jù)0RF中每個設備的經(jīng)煒度信息，以及當前雷達的地理坐標信息，提取各 個設備所對應地理位置的CAPPI_MAX值。為便于表述，不妨定義CM為CAPPI_MAX的一個子集， CM中的每個值CM[ c ]為與0RF[ c ]最相鄰地理位置的回波強度值。顯然，CM數(shù)據(jù)集中的記錄數(shù) 與0RF相同。設置兩個遞增變量PA和Pb，其中PAe [ 100，400]，以10為步長遞增，Pb e [ 1. 0， 2.0]，以0.1為步長遞增，分別計算0RF與CM中各記錄項的誤差平方和，計算方法為：
[0059]
[0060] 步驟8、對每一組PA和Pb，都可計算出一個D(PA，Pb)，記錄下D(PA，Pb)最小時的一 組PA和Pb，記為。
[0061 ] 步驟9、將計算得到的？6&31：1'和?&抑_4和?&抑_13，添加到數(shù)據(jù)庫的？6&1：11代111；1^〇數(shù) 據(jù)表中，該數(shù)據(jù)表的結構如表1所;
[0062] 表 1
[0063]
[0064] 上述步驟1-9的過程如圖1所示。
[0065] 步驟10、依據(jù)上述步驟1-步驟9,將本實施例所選用的雷達和地面雨量觀測的數(shù)據(jù) 逐一進行計算并入庫。完成后，即得到該部雷達的反射率特征庫。
[0066] 建立好雷達反射率特征庫后，就可以在實際業(yè)務應用中，將實時、最新的雷達基數(shù) 據(jù)文件通過反射率特征提取、特征匹配，找到最佳的Z-R關系的參數(shù)A和參數(shù)b，從而進行雷 達降水估測。仍以這部南京的雷達為例，具體流程如下：
[0067] 步驟11、將最新接收到的雷達基數(shù)據(jù)文件按上述步驟1-步驟5,計算出回波特征字 符串，記為FeaStrCur。
[0068] 步驟12、從雷達反射率特征庫中遍歷具有相同RadarlD(即這部南京的雷達）的 FeaStr，形成數(shù)據(jù)集FS，F(xiàn)S [ c ]表示FS中的c條記錄。將FS中的逐條記錄FS [ c ]與FeaStrCur作 相關性計算，計算方法為：
[0069]
[0070] 其中，F(xiàn)eaStrCur和F'S.[c]分別表示FeaStrCur與FS[c]中所有記錄的算術平均值。 由于FeaStr字符串的記錄方式是每2個字節(jié)為一個特征值，用FS[c](k)表示第k項記錄值， 且該值是16進制的數(shù)值，因此，在進行相關性R(c)計算前，先將數(shù)據(jù)FS[c](k)和FeaStrCur (k)轉換為10進制，再進行計算。
[0071 ] 步驟13、找到R( c)值最大時FS [ c ]所對應的記錄，該記錄中對應的 即為當前雷達資料估測降水時Z-R關系的最優(yōu)參數(shù)A和參數(shù)b，由此計算出的降水量具有更 高的準確性。需要補充說明的是，考慮到本方法投入實施前期，雷達反射率特征庫中的特征 信息記錄可能并不完善，又或者當前天氣現(xiàn)象和雷達回波信息比較特殊，從雷達反射率特 征庫中匹配計算出的相關系數(shù)R(c)都非常低，致使特征庫中的參數(shù)A和參數(shù)b不適用于當前 Z-R關系。因此，限定當最大的R(c)〈0.5時，Z-R關系的參數(shù)A和參數(shù)b改為使用預先設定好的 參數(shù)值，如 2 = 3001^^2 = 2301^25 或 2 = 2001^·4 等。
[0072] 上述步驟11-13的過程如圖2所示。
[0073]以上實施例僅為說明本發(fā)明的技術思想，不能以此限定本發(fā)明的保護范圍，凡是 按照本發(fā)明提出的技術思想，在技術方案基礎上所做的任何改動，均落入本發(fā)明保護范圍 之內。
【主權項】
1.基于天氣雷達反射率特征匹配的降水估測方法，其特征在于，包括W下步驟： 步驟一、構建天氣雷達反射率特征庫： (1) 提取天氣雷達基數(shù)據(jù)中的基本反射率數(shù)據(jù)，并將基本反射率數(shù)據(jù)從原來的極坐標 形式轉換為平面直角坐標系形式；分別計算P千米和q千米高度處的等高平面位置顯示圖 CAPPI_p和CAPPI_q，CAPPI_p和CAPPI_q的結構相同，均為NXN的二維網(wǎng)格，求取CAPPI_p和 CAPPI_q各網(wǎng)格點在同一水平位置的最大回波強度值CAPPI_MAX，即CAPPI_MAX(X，y) =Max (CAPPI_p(X，y)，CAPPI_q(X，y))，其中X、y分別表示平面直角坐標系下的橫坐標和縱坐標， Max為取最大值函數(shù)； (2) 根據(jù)得到的最大回波強度值CAPPI_MAX(x，y)，計算雷達回波的灰度圖像Rad(x，y):(3) 計算灰度圖像Rad(x，y)的一階梯度值，從而強化回波圖像的變化特征； (4) 新建一個MXM的雷達回波特征矩陣化a，M<N/2Jea中每個單元的值計算方式為：其中，x'e陽;，Μ - 1]巧'￡[〇，1-1]，胎義為對括號中的值作16進制轉換的 函數(shù)； (5) 將得到的雷達回波特征矩陣化曰，從Fea(0,0巧化ea(M-l，M-l)逐個添加到一個字符 串變量化aStr中； (6) 選取雷達有效探測范圍內的地面雨量觀測設備，W當前雷達探測時間為基準，提取 運些地面雨量觀測設備在基準后的X小時內的累積觀測降水量，形成數(shù)據(jù)集ORF，其中ORF k]表示第C個地面雨量觀測設備的X小時降水量，xe [1,24]，將當前雷達有效探測范圍內 地面雨量觀測設備的總數(shù)記為MS皿，即數(shù)據(jù)集ORF的大小為MS皿，則C e [0 ,MS皿-1 ]; (7) 根據(jù)數(shù)據(jù)集ORF中每個地面雨量觀測設備的經(jīng)締度信息，W及當前雷達的地理坐標 信息，提取各個地面雨量觀測設備所對應地理位置的CAPPI_MAX值;定義CM為CAPPI_MAX的 一個子集，CM中的每個值CMk]為與ORFk]地理位置最相鄰的網(wǎng)格點的回波強度值;設置兩 個遞增變量PA和Pb，其中?4￡[100,400]，^10為步長遞增，？6￡[1.0,2.0]，^0.1為步長遞 增，分別計算0RF與CM中各記錄項的誤差平方和：(8) 記錄下所得D (PA，Pb)取值最小時的一組PA和Pb，分別記為?日拘_4和?日拘_6; (9) 將步驟(5)中得到的化aStr和步驟(8)得到的？日拘_4和?日拘_6添加到數(shù)據(jù)庫； (10) 依據(jù)上述步驟(1)-(9),將有資料w來，雷達探測數(shù)據(jù)和地面雨量觀測數(shù)據(jù)逐一進 行計算并入庫，從而構建起天氣雷達反射率特征庫； 步驟二:將實況天氣雷達反射率數(shù)據(jù)與建立的天氣雷達反射率特征庫進行匹配： (11) 將需要匹配的天氣雷達反射率數(shù)據(jù)按照上述步驟（1)-(5)，計算出能夠反映該雷 達回波特征的字符串化aStr化。 (12) 從天氣雷達反射率特征庫中遍歷同一部天氣雷達的FeaS化，形成數(shù)據(jù)集FS，!^S中 的每條記錄FSk ]與化aStr化r作相關性計算，計算方法為：其中，晚a託蛇此和RS [C]分別表示FeaStrCur與FS [ C ]中所有記錄的算術平均值，F(xiàn)S k ] 化)表示FSk]的第k項記錄值； (13) 找到R(c)取值最大時，F(xiàn)SU]所對應的數(shù)據(jù)庫記錄，該記錄中對應的化^_八和 Para_b即為當前天氣雷達資料下Z-R關系的最優(yōu)參數(shù)A和參數(shù)b，并據(jù)此估測降水。2. 根據(jù)權利要求1所述基于天氣雷達反射率特征匹配的降水估測方法，其特征在于:在 步驟(1)中，取P二 1.5,q = 3.0。3. 根據(jù)權利要求1所述基于天氣雷達反射率特征匹配的降水估測方法，其特征在于:在 步驟(3)中，灰度圖像Rad(x，y)最邊緣一圈的各像素的一階梯度為0,其余像素的一階梯度 采用Sobel檢測算子來計算。4. 根據(jù)權利要求3所述基于天氣雷達反射率特征匹配的降水估測方法，其特征在于:橫 向Sobel檢測算子5. 根據(jù)權利要求1所述基于天氣雷達反射率特征匹配的降水估測方法，其特征在于:在 步驟(5)中，字符串變量化aStr中每2個字節(jié)代表1個化a (X '，y '）。6. 根據(jù)權利要求1所述基于天氣雷達反射率特征匹配的降水估測方法，其特征在于:在 步驟(13)中，若找到的最大R(cK〇.5時，則Z-R關系中的最有參數(shù)A和參數(shù)b的取值改用預先 設定值。
【文檔編號】G01S13/95GK105974418SQ201610540931
【公開日】2016年9月28日
【申請日】2016年7月8日
【發(fā)明人】王興, 苗春生, 王堅紅, 錢代麗, 王麗娟
【申請人】南京信息工程大學