Takashi Unuma's blog.

Home About Archive

降水強度のKDP推定アルゴリズムについて

はじめに

気象レーダーが単偏波の場合は Z-R 関係と呼ばれる関係式が用いられ、反射強度から降水強度が推定される。ところが、Signal Processing Methods で述べたように、水平偏波のみでは主に降雨減衰などの影響を受けやすく、その降水強度推定時の精度が低下する場合がある。一方、二重偏波の場合、水平偏波に加えて垂直偏波による観測があり、それらの位相差を用いると前述の影響を受けにくく高精度に降水強度を推定できうるという利点がある。加えて、二重偏波要素の特徴で具体例を示したように、二重偏波要素を用いた反射強度等の品質向上にも繋がる。
　降水強度の推定については二重偏波要素を用いた様々な方法が提案されている。近年では、偏波間位相差変化率 ( $K_{\rm DP}$ ) を用いた降水強度の推定手法が提案されており、主流となってきている。偏波間位相差変化率 ( $K_{\rm DP}$ )の算出においては、ライブラリ化されたものやオープンソースとなっているものがある。Reimel and Kumjian (2021) は、Known-Truth Framework を用い、公開されている幾つかの $K_{\rm DP}$ 推定アルゴリズムの性能やそれぞれの特徴について調べられている。本稿では、この論文をレビューしつつ各 $K_{\rm DP}$ 推定アルゴリズムの特徴を把握し理解を深めることを目的とする。

用語

gate: レンジ方向 (レーダーのアンテナからビームが放射される方向) のデータグリッド単位
$\Phi_{\rm DP}$ : 偏波間位相差
$K_{\rm DP}$ : 偏波間位相差変化率
$Z_{H}$ : 水平偏波反射強度
$\delta$ $δ$ : 後方散乱偏波間位相差
- Backscatter Differential Phase
- ノイズ全般を指す場合もある
FIR: 有限インパルス応答
- Finite Impulse Response
KF: カルマンフィルター
- Kalman Filter

$K_{\rm DP}$ 推定アルゴリズム

NWS OPERATIONAL

Ryzhkov et al. (2005)

National Weather Service の現業で利用されているアルゴリズム。まず、2 つ (9-gate 移動平均・25-gate 移動平均) の平滑化された $\Phi_{\rm DP}$ プロファイルを用意。それぞれの平滑化された $\Phi_{\rm DP}$ から $K_{\rm DP}$ を最小二乗法で求める。各 gate で 2 つの $K_{\rm DP}$ をそれぞれ保持。 $Z_{H}$ が各 gate で計算された上で、 $Z_{H} \ge$ 40 dBZ ならば 9-gate $K_{\rm DP}$ が、 $Z_{H} \lt$ 40 dBZ ならば 25-gate $K_{\rm DP}$ がそれぞれ用いられる、というもの。

HUBBERT AND BRINGI FIR FILTER

Hubbert and Bringi (1995)

ノイズの平滑化と $\delta$ 除去をそれぞれ行うために、有限インパルス応答フィルターを $\Phi_{\rm DP}$ に繰り返し適用する手法。まず、 $\delta$ が無いと仮定し、ノイズを含む $\Phi_{\rm DP}$ へ一度だけ 21-gate フィルターを適用。21 の係数をレンジ方向に連続した $\Phi_{\rm DP}$ に掛け合わせ、適用 (表は省略)。プロダクトは積算の後、FIR フィルター係数の和で規格化される。値は、21-point フィルターの中心 gate に適用される。フィルターは 1-gate づつ前方へ移動し、全ての場に対して繰り返し実施される。以上の処理済み $\Phi_{\rm DP}$ を用い $K_{\rm DP}$ を求める、というもの。

CALC_KDP_BRINGI (CSU_RADAR_TOOLS)

Lang et al. (2007)
- CSU RadarTools: calc_kdp_bringi function

観測された $\Phi_{\rm DP}$ を平滑化し、それを用い $K_{\rm DP}$ を推定する手法。まず、 $\Phi_{\rm DP}$ の標準偏差 s( $\Phi_{\rm DP}$ ) をユーザーが定義する windows 幅で計算。s( $\Phi_{\rm DP}$ ) がユーザー定義閾値よりも大きければ $\Phi_{\rm DP}$ はマスクされる。次に、Hubbert and Bringi (1995) に似た FIR フィルターをかける。オリジナルとの違いは、フィルター幅と繰り返し回数をユーザーが定義出来るようになっていること。 $\Phi_{\rm DP}$ の勾配の半分を用い、 $Z_{H}$ に逆比例する window 幅で $K_{\rm DP}$ を計算・推定する、というもの。

PHASE_PROC_LP (Py-ART)

Giangrande et al. (2013)
- Py-Art: phase_proc_lp function

単調な $\Phi_{\rm DP}$ に対する線型計画問題を定義し、偏波 self-consistency 拘束条件を与えるアルゴリズム。0 $^{\circ}{\rm C}$ 層以下で $\delta$ の影響を受けていない、なめらかな $\Phi_{\rm DP}$ を出力することを想定。Sobel フィルターを用い平滑化した場から $K_{\rm DP}$ が計算される。

KDP_SCHNEEBELI (Py-ART)

Schneebeli et al. (2014)
- Py-Art: kdp_schneebeli function

\Phi_{\rm DP}

の処理に Kalman filter (KF) を用い、その結果から

K_{\rm DP}

を算出する手法。カルマンフィルタでは、処理前の

\Phi_{\rm DP}

データが観測値として用いられ、

K_{\rm DP}

\Phi_{\rm DP}

\delta

間の所定の関係が、観測値から場の変数への遷移に用いられる。カルマンフィルタはそれぞれの gate に対して 2 回行われる。1 回目はレーダーから遠ざかる方向(前方処理)、2 回目は観測範囲終端からレーダーに向う方向(後方処理)。このプロセスはそれぞれの gate に対して複数回繰り返される。各々の実行時には、固有の共分散行列の値が含まれ、個々のカルマンフィルタのアンサンブルを推定する。それぞれの gate における最終的な

K_{\rm DP}

値は、前方・後方処理のアンサンブル平均、もしくは 2 つのカルマンフィルタ推定値の平均が適用される。このとき、アンサンブル平均が gate 間で

0.15^{\circ}\, {\rm km}^{-1}

を上回り増加する(下回り減少する)場合は、前方(後方)処理によるカルマンフィルタ推定値が用いられる。勾配がある場所で

-0.15^{\circ}\, {\rm km}^{-1}

と

+0.15^{\circ}\, {\rm km}^{-1}

の間の値をとる場合は、前方・後方処理双方によるカルマンフィルタ推定値の平均が用いられる。アンサンブルサイズは

K_{\rm DP}

の標準偏差に依存する。

K_{\rm DP}

の標準偏差が大きいほどアンサンブルサイズの大きな値が最終推定値に用いられる。

KDP_VULPIANI (Py-ART)

Vulpiani et al. (2012, 2015)
- Py-Art: kdp_vulpiani function

下記 4 step を指定回数分繰り返し、 $K_{\rm DP}$ を推定するアルゴリズム。1) $\Phi_{\rm DP}$ の移動平均した値から $K_{\rm DP}$ を有限差分によって求める。2) $K_{\rm DP}$ の有効値を閾値と比較。波長帯によって上・下限値が異なる。有効値の範囲外となる場合は、全て $0^{\circ}\, {\rm km}^{-1}$ とする。3) $\Phi_{\rm DP}$ を処理済み $K_{\rm DP}$ から復元。4) 最終 $K_{\rm DP}$ 推定値を $\Phi_{\rm DP}$ の有限差分からもう一度計算。

KDP_MAESAKA (Py-ART)

Maesaka et al. (2012)
- Py-Art: kdp_maesaka function

平滑化した $\Phi_{\rm DP}$ を元の $\Phi_{\rm DP}$ に一致させるように Cost function を用いると同時に、low-pass フィルターを適用するアルゴリズム。フィルターの重みは $C_{\rm lpf}$ として与え、 $C_{\rm lpf}$ の値が大きいほど積極的に滑らかにするよう作用する。 $K_{\rm DP}$ の値は cost function が最小値となる場合に一意に求まる。0 $^{\circ}{\rm C}$ 層より下で雨についてのみ適用可能。

各アルゴリズムの比較

$\Phi_{\rm DP}$ $Φ_{D P}$ を平滑化する際、平滑化のやり方 (移動平均・フィルター) や度合い (適用 gate 数・window 幅) がそれぞれ異なる。
- 移動平均: NWS OPERATIONAL
- FIR フィルター: HUBBERT AND BRINGI FIR FILTER, CALC_KDP_BRINGI
- Sobel フィルター: PHASE_PROC_LP
- Low-pass フィルター: KDP_MAESAKA
$\delta$ $δ$ の除去を念頭に置いたアルゴリズム。
- HUBBERT AND BRINGI FIR FILTER, PHASE_PROC_LP, KDP_SCHNEEBELI.
観測された $\Phi_{\rm DP}$ $Φ_{D P}$ に処理済み $\Phi_{\rm DP}$ $Φ_{D P}$ の値を近づけるために、Kalman Filter や Cost Function を用いるアルゴリズムがある。
- KDP_SCHNEEBELI, KDP_MAESAKA.
$K_{\rm DP}$ $K_{D P}$ 推定後に $\Phi_{\rm DP}$ $Φ_{D P}$ を復元するといった self-consistency を念頭に置いたアルゴリズムがある。
- PHASE_PROC_LP, KDP_VULPIANI.

$K_{\rm DP}$ 推定結果の確認

基本的に $K_{\rm DP}$ は $\Phi_{\rm DP}$ の有限差分 (微分) から計算される。このため、直感的にはレンジ (ray) 方向のグラフで $\Phi_{\rm DP}$ の値の変化が大きくなるとともに $K_{\rm DP}$ の値も大きく変化していることが期待される。

図は、Py-ART のサンプルデータを用い、このページで紹介されているサンプルスクリプトを改変した下記で描画したものである。

# 
# plot_ray.py
#

import numpy as np
import matplotlib.pyplot as plt
import pyart

# Read sample data which got from the URL: https://engineering.arm.gov/~jhelmus/pyart_example_data/095636.mdv.
radar = pyart.io.read_mdv('095636.mdv')

# Specify the first sweep in the volume scan
radar = radar.extract_sweeps([0])

# Set a field for a single ray
display = pyart.graph.RadarDisplay(radar)
fig = plt.figure(figsize=[10, 5])
ax = fig.add_subplot(111)

# Set a ray number
ray_num = 200

# Plot original phidp
display.plot_ray('differential_phase', ray_num, format_str='b-',
                 axislabels_flag=False, title_flag=False, ax=ax)

# Set labels
ax.set_ylim(-180, 180)
ax.set_ylabel('Differential phase shift (degrees)')
ax.set_xlabel('Range (km)')

# Plot KDP on second axis
ax2 = ax.twinx()
ax2.set_ylabel('Specific differential phase shift (degrees/km)')
display.plot_ray('specific_differential_phase', ray_num, format_str='r-',
                 axislabels_flag=False, title_flag=False, ax=ax2)

# Decorate
ax2.yaxis.grid(color='gray', linestyle='dashed')
ax.legend(display.plots, ["phiDP", 'KDP'], loc='upper right')

# Display
plt.show()

青実線が $\Phi_{\rm DP}$ 、赤実線が $K_{\rm DP}$ を示す。サンプルデータの $\Phi_{\rm DP}$ は特にフィルター処理されていないため、グラフはギザギザになっている。レンジ方向 (横軸) 70–80 km のあたりを見ると、 $\Phi_{\rm DP}$ の勾配が大きい部分と $K_{\rm DP}$ の値が大きい部分とが一致している。一方、レンジ方向 30–50 km のあたりを見ると、特に 40–50 km 付近で $\Phi_{\rm DP}$ の変動があるもののほぼ平滑化されてしまっており、 $K_{\rm DP}$ の値はほぼ 0 になっている。こういった部分 ( $\Phi_{\rm DP}$ の平滑化や $K_{\rm DP}$ 負値の取り扱い) を適切に処理する必要があると考えられる。そのために、様々な工夫がなされつつ多様な $K_{\rm DP}$ 推定方法が提案されているものと推察する。

参考文献

Giangrande, S. E., R. McGraw, and L. Lei, 2013: An application of linear programming to polarimetric radar differential phase processing. J. Atmos. Oceanic Technol., 30, 1716–1729, doi:10.1175/JTECH-D-12-00147.1.
Hubbert, J., and V. N. Bringi, 1995: An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements. J. Atmos. Oceanic Technol., 12, 643–648, doi:10.1175/1520-0426(1995)012<0643:AIFTFT>2.0.CO;2.
Lang, T. J., D. A. Ahijevych, S. W. Nesbitt, R. E. Carbone, S. A. Rutledge, and R. Cifelli, 2007: Radar-observed characteristics of precipitating systems during NAME 2004. J. Climate, 20, 1713–1733, doi:10.1175/JCLI4082.1.
Maesaka, T., K. Iwanami, and M. Maki, 2012: Non-negative KDP estimation by monotone increasing ΦDP assumption below melting layer. Seventh European Conf. on Radar in Meteorology and Hydrology, Toulouse, France, ERAD, PDF.
Reimel, K. J., and M. R. Kumjian, 2021: Evaluation of KDP estimation algorithm performance in rain using a known-truth framework. Journal of Atmospheric and Oceanic Technology, 38, 587–605, doi:10.1175/JTECH-D-20-0060.1.
Ryzhkov, A. V., T. J. Schuur, D. W. Burgess, P. L. Heinselman, S. E. Giangrande, and D. S. Zrnić, 2005: The Joint Polarization Experiment: Polarimetric rainfall measurements and hydrometeor classification. Bull. Amer. Meteor. Soc., 86, 809–824, doi:10.1175/BAMS-86-6-809.
Schneebeli, M., J. Grazioli, and A. Berne, 2014: Improved estimation of the specific differential phase shift using a compilation of Kalman filter ensembles. IEEE Trans. Geosci. Remote Sens., 52, 5137–5149, doi:10.1109/TGRS.2013.2287017.
Vulpiani, G., M. Montopoli, L. D. Passeri, A. G. Gioia, P. Giordano, and F. S. Marzano, 2012: On the use of dual-polarized C-band radar for operational rainfall retrieval in mountainous areas. J. Appl. Meteor. Climatol., 51, 405–425, doi:10.1175/JAMC-D-10-05024.1.
Vulpiani, G., L. Baldini, and N. Roberto, 2015: Characterization of Mediterranean hail-bearing storms using an operational polarimetric X-band radar. Atmos. Meas. Tech., 8, 4681–4698, doi:10.5194/amt-8-4681-2015.

更新履歴

2021-07-28: 初稿
2021-08-05: 誤字脱字を修正
2021-08-07: $K_{\rm DP}$ 推定結果の確認について追記・文章を少し修正
2021-08-19: 体裁微修正・アルゴリズムのリンク先修正

Miscellaneous — Jul 28, 2021

View / Make Comments

Made with ❤ and at Japan.