diff --git a/.DS_Store b/.DS_Store index 4afb6ce..a9b1919 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/2021/03/20/index.html b/2021/03/20/index.html index 9b7218d..7b38f4f 100644 --- a/2021/03/20/index.html +++ b/2021/03/20/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/03/21/index.html b/2021/03/21/index.html index a0cf711..8a875ba 100644 --- a/2021/03/21/index.html +++ b/2021/03/21/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/03/22/index.html b/2021/03/22/index.html index 7284d34..4cf3b12 100644 --- a/2021/03/22/index.html +++ b/2021/03/22/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/03/23/index.html b/2021/03/23/index.html index 9ec8fc7..9a8e71e 100644 --- a/2021/03/23/index.html +++ b/2021/03/23/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/04/28/index.html b/2021/04/28/index.html index 980ca2b..f5d143e 100644 --- a/2021/04/28/index.html +++ b/2021/04/28/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/07/28/index.html b/2021/07/28/index.html index 373f22b..e3962da 100644 --- a/2021/07/28/index.html +++ b/2021/07/28/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/08/16/index.html b/2021/08/16/index.html index f5002c5..6ae5c18 100644 --- a/2021/08/16/index.html +++ b/2021/08/16/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/08/19/index.html b/2021/08/19/index.html index 95eb874..2d5d5c4 100644 --- a/2021/08/19/index.html +++ b/2021/08/19/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/08/28/index.html b/2021/08/28/index.html index dc45a3b..2597c07 100644 --- a/2021/08/28/index.html +++ b/2021/08/28/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/09/04/index.html b/2021/09/04/index.html index bde1738..b1305a5 100644 --- a/2021/09/04/index.html +++ b/2021/09/04/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/09/20/index.html b/2021/09/20/index.html index fc5f249..730fc17 100644 --- a/2021/09/20/index.html +++ b/2021/09/20/index.html @@ -128,7 +128,7 @@ - + diff --git a/2021/10/01/index.html b/2021/10/01/index.html index 993aeee..77d38bb 100644 --- a/2021/10/01/index.html +++ b/2021/10/01/index.html @@ -128,7 +128,7 @@ - + diff --git a/2022/01/26/index.html b/2022/01/26/index.html index c3201e9..1273914 100644 --- a/2022/01/26/index.html +++ b/2022/01/26/index.html @@ -128,7 +128,7 @@ - + diff --git a/2022/06/18/index.html b/2022/06/18/index.html index 6b913fa..e440cd6 100644 --- a/2022/06/18/index.html +++ b/2022/06/18/index.html @@ -128,7 +128,7 @@ - + diff --git a/2022/08/05/index.html b/2022/08/05/index.html index 8f9efa7..5c1a406 100644 --- a/2022/08/05/index.html +++ b/2022/08/05/index.html @@ -128,7 +128,7 @@ - + diff --git a/2022/10/12/index.html b/2022/10/12/index.html index d7ea814..d74350e 100644 --- a/2022/10/12/index.html +++ b/2022/10/12/index.html @@ -128,7 +128,7 @@ - + diff --git a/2022/11/19/index.html b/2022/11/19/index.html index d3c8170..ca7f8e2 100644 --- a/2022/11/19/index.html +++ b/2022/11/19/index.html @@ -128,7 +128,7 @@ - + diff --git a/2024/01/01/index.html b/2024/01/01/index.html index 595f9a2..95f8ebb 100644 --- a/2024/01/01/index.html +++ b/2024/01/01/index.html @@ -128,7 +128,7 @@ - + diff --git a/2024/01/02/index.html b/2024/01/02/index.html index 4c07001..96e8b14 100644 --- a/2024/01/02/index.html +++ b/2024/01/02/index.html @@ -128,7 +128,7 @@ - + diff --git a/2024/01/19/index.html b/2024/01/19/index.html index 9d1bc4a..dba2b00 100644 --- a/2024/01/19/index.html +++ b/2024/01/19/index.html @@ -128,7 +128,7 @@ - + @@ -327,10 +327,10 @@

出力結果から粒径分布を計算

# 粒径 (半径) [unit: cm] 定数データ読み込み # index とカテゴリの対応は下記。 -# 0: aerosol -# 1: cloud + rain -# 2: ice + snow -# 3: graupel + hail +# 0: cloud + rain +# 1: ice/columns +# 2: ice/plates +# 3: ice/dendrites # 4: snow # 5: graupel # 6: hail @@ -460,6 +460,7 @@

更新履歴

  • 2024-01-20: 誤字を修正
  • 2024-01-23: 表記を修正
  • 2024-02-01: Registory の URL を修正
    • +
  • 2024-05-26: 定数ファイル読み込み時の index とカテゴリの対応表を修正

    • — Jan 19, 2024

    diff --git a/2024/01/20/index.html b/2024/01/20/index.html index cfb279d..77aeb36 100644 --- a/2024/01/20/index.html +++ b/2024/01/20/index.html @@ -132,7 +132,7 @@ - + diff --git a/2024/01/27/index.html b/2024/01/27/index.html index ff46116..094f531 100644 --- a/2024/01/27/index.html +++ b/2024/01/27/index.html @@ -128,7 +128,7 @@ - + diff --git a/about/index.html b/about/index.html index 809583f..fed926e 100644 --- a/about/index.html +++ b/about/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/03/index.html b/archives/2021/03/index.html index 26480d2..54168ce 100644 --- a/archives/2021/03/index.html +++ b/archives/2021/03/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/04/index.html b/archives/2021/04/index.html index 3a9dcc4..190f2d5 100644 --- a/archives/2021/04/index.html +++ b/archives/2021/04/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/07/index.html b/archives/2021/07/index.html index f4a0834..c63acec 100644 --- a/archives/2021/07/index.html +++ b/archives/2021/07/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/08/index.html b/archives/2021/08/index.html index 7fda93e..09aef23 100644 --- a/archives/2021/08/index.html +++ b/archives/2021/08/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/09/index.html b/archives/2021/09/index.html index 6547f6f..6d1016a 100644 --- a/archives/2021/09/index.html +++ b/archives/2021/09/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/10/index.html b/archives/2021/10/index.html index 42a7706..9afed32 100644 --- a/archives/2021/10/index.html +++ b/archives/2021/10/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/index.html b/archives/2021/index.html index 3af9ec9..e3ad1cb 100644 --- a/archives/2021/index.html +++ b/archives/2021/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2021/page/2/index.html b/archives/2021/page/2/index.html index 2ad49e0..a03e8ff 100644 --- a/archives/2021/page/2/index.html +++ b/archives/2021/page/2/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/01/index.html b/archives/2022/01/index.html index bc34fc6..0192fd6 100644 --- a/archives/2022/01/index.html +++ b/archives/2022/01/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/06/index.html b/archives/2022/06/index.html index eabf2e3..60fc7a1 100644 --- a/archives/2022/06/index.html +++ b/archives/2022/06/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/08/index.html b/archives/2022/08/index.html index b4a3040..ebdaaff 100644 --- a/archives/2022/08/index.html +++ b/archives/2022/08/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/10/index.html b/archives/2022/10/index.html index 5953495..67344f2 100644 --- a/archives/2022/10/index.html +++ b/archives/2022/10/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/11/index.html b/archives/2022/11/index.html index e42c176..afb75ba 100644 --- a/archives/2022/11/index.html +++ b/archives/2022/11/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2022/index.html b/archives/2022/index.html index c0c41fb..250548a 100644 --- a/archives/2022/index.html +++ b/archives/2022/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2024/01/index.html b/archives/2024/01/index.html index f4afd05..0aa8434 100644 --- a/archives/2024/01/index.html +++ b/archives/2024/01/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/2024/index.html b/archives/2024/index.html index 894c3b5..f9235a9 100644 --- a/archives/2024/index.html +++ b/archives/2024/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/index.html b/archives/index.html index 8ba4a42..258ae66 100644 --- a/archives/index.html +++ b/archives/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/page/2/index.html b/archives/page/2/index.html index 608aa2b..ccb1f4e 100644 --- a/archives/page/2/index.html +++ b/archives/page/2/index.html @@ -128,7 +128,7 @@ - + diff --git a/archives/page/3/index.html b/archives/page/3/index.html index a438097..71dcd84 100644 --- a/archives/page/3/index.html +++ b/archives/page/3/index.html @@ -128,7 +128,7 @@ - + diff --git a/atom.xml b/atom.xml index 139d3a9..0598bfa 100644 --- a/atom.xml +++ b/atom.xml @@ -6,7 +6,7 @@ - 2024-03-09T11:08:09.734Z + 2024-04-15T14:13:04.492Z https://takashiunuma.github.io/blog/ @@ -21,7 +21,7 @@ https://takashiunuma.github.io/blog/2024/01/27/ 2024-01-27T06:15:47.000Z - 2024-03-09T11:08:09.734Z + 2024-04-15T14:13:04.492Z はじめに

    雨滴粒径分布の話題は、レーダー気象学と関連するものとして Self-consistency Principle、粒径分布とレーダー変数を繋ぐ関係式、或いはその逆で粒径分布パラメータのレーダー変数による推定について取り上げてきた。これらの記事では、仮定する粒径分布関数の形は基本的に既知であるものとして扱ってきた。ところが、実際に観測される雨滴粒径分布というのは、我々が思っている以上に多様である。これは、Spectral Bin Microphysics Model で示したように、数値モデル内で計算される粒径分布からもその多様性が認知されうる。そこで、本稿では複数のピークを持つ雨滴粒径分布の形成過程に着目し、主に McFarquhar (2010) のレビューを行うことで理解を深める。

    複数のピークを持つ雨滴粒径分布

    前述した通り、雨滴粒径分布は複数のピークを持つ場合があることが知られている (Valdez and Young 1985; Brown 1986; Feingold et al. 1988; List and McFarquhar 1990; Chen and Lamb 1994)。中でも、2 つのピークをもつ雨滴粒径分布について室内実験或いは地上観測により調べられてきた (Steiner and Waldvogel 1987; List et al. 1987; Asselin de Beauville et al. 1988; Zawadzki and de Agostinho Antonio 1988)。これらの複数ピークを持つ雨滴粒径分布の形成には、雨滴同士の衝突・併合・分裂が主な要因として重要視されている。これらの物理過程を数値的に扱い十分な時間を経た結果、「定常的な分布 (stationary distribution)」が得られる。この定常的な分布はあくまで数値的に表現されるものであって、自然界で観測的或いは物理的に表現されるものとは異なると考えられている (McFarquhar 2010)。

    雨滴同士の衝突・併合が継続して生じれば、より大きな雨滴が生じることは容易に想像できるだろう。しかしながら、雨滴にも大きさの限界があり、ある一定の大きさになると分裂してしまう。この分裂にはいくつかのパターンがあることが知られている。具体的には、filament (neck) 型・sheet 型・disk 型・Bag 型である (McTaggart-Cowan and List 1975; Low and List 1982a)。前述の数値的に得られる定常的な分布は、この分裂過程と併合過程とのバランスによって表現される (Srivastava 1971; Young 1975; List and Gillespie 1976)。

    雨滴の分裂過程を取り入れた数値的研究の多くは、Low and List (1982b) による breakup kernel が用いられた。この kernel を用いると、定常あるいは非定常な降雨に対して 3 つのピークを持つ定常的な分布が得られる。ところが、より物理的な分裂過程を加味した kernel を用いると、2 つのピークを持つ定常的な分布が得られる (McFarquhar 2004a)。さらに、Schlottke et al. (2010) が 32 個の雨滴ペアの衝突過程を数値的に調べ、その結果を踏まえて構築した kernel を用いると、ピークの位置がやや異なるものの McFarquhar (2004a) と同様に 2 つのピークを持った定常的な分布が得られる (Straub et al. 2010)。いずれの場合も、観測事実として特有の分布が得られている訳では無い。

    雨滴の衝突頻度は、弱雨の場合に定常分布に関連するピークを作り出すほど頻繁ではないと考えられる (e.g., Zawadzki et al. 1994)。一方、雨滴の衝突頻度は非常に激しい雨量の場合には数値的に得られたような定常的な分布を作り出せる可能性がある (McFarquhar and List 1991a, 1991b; McFarquhar 2004b)。このようなピークが生じた観測結果は、いくつかの研究 (Willis 1984; McFarquhar et al. 1996; Garcia-Garcia and Gonzalez 2000) で指摘されているものの、定常的な分布でのピークは一貫して観測される訳では無い。このことから、数値モデルと観測結果との不一致の要因を調べる必要がある。

    なぜ数値的に得られる定常的な分布は自然界で一貫して観測されないのか？

    McFarquhar (2004b) は、衝突する雨滴の特定のペアによって生成されるフラグメント分布の大きな広がりが、空間内の特定の場所と時間でサンプリングされた粒径分布に大きな変動を引き起こすという仮説を立てた。このことは、自然界で定常分布が観察されない理由や、衝突相互作用率を高める雨滴のクラスターの役割を考慮しても、大雨で同じ特定の場所にピークが見られない理由を説明することになる (McFarquhar 2004b)。また、定常分布に近づくのに必要な時間は、初期の粒径分布によって大きく変化する (e.g., Prat and Barros 2007)。加えて、雨滴の蒸発 (Brown 1993; Hu and Srivastava 1995) や周辺大気環境条件としての上昇気流による size sorting (Kollias et al. 2001) など、粒径分布への影響が局所的な条件に依存する他の要因もある。

    ここで、定式化された数値モデルの準統一的な仮定を正当化するために必要なサンプル数を考える。50 mm h−1^{-1} を超える降水強度の条件下で、いくつかのイベントにわたって積分された地上の粒形分布観測データを数値実験の結果と比較するためには、約 6 時間のデータが必要であることが推定されている (McFarquhar 2004b)。さまざまなフィールドキャンペーンからのデータを組み合わせれば、そのようなデータセットは存在しうる。今後、個々のフィールドキャンペーンからのデータ解析だけでなく、様々なプログラムから収集されたデータの統合的な解析にも注力が必要である (McFarquhar 2010)。

    最近の研究動向

    ここからは、McFarquhar (2010) で議論された内容の内、近年の研究で明らかになってきていることについて述べる。

    D’Adderio et al. (2015) は、数値的に表現される定常的な分布の一つの特徴である衝突分裂のシグナルを捉えるアルゴリズムを開発し、いくつかの地域で行われた観測データを用いてその特徴が調べられた。さらに、D’Adderio et al. (2018) は、弱い雨に対してもアルゴリズムを適用し、頻度は低いものの衝突分裂の特徴を有する定常的な分布が自然界でも観測される場合があることを明らかにした。加えて、弱雨や層状性の降雨の場合には定常的な分布になりにくく、強雨や対流性の降雨ほど定常的な分布になりやすいことを示した。

    Kumjian and Prat (2014) は、鉛直一次元の Spectral bin microphysics model を用いることで粒径分布の鉛直変化に対する影響を調べ、個々の雲微物理過程は二重偏波レーダー変数の空間パラメーター上で遷移することを明らかにした。このように、粒径分布が鉛直方向に変化していることに対しても考慮が必要と考えられる。Kobayashi and Adachi (2001) は、400 MHz 帯のウィンドプロファイラーを用いた粒径分布の推定結果から地上で強雨がもたらされる際には粒径分布の鉛直分布に違いがあり、地上付近では衝突分裂のシグナルが捉えられている可能性があることを示した。近年の二重偏波レーダーによる観測結果では、地上付近で衝突分裂が支配的なシグナルを示す一方、その上空では衝突併合のシグナルが支配的となることが確認されてきている (Jung and jou 2023; Unuma 2024)。

    雨滴の蒸発は、周辺の熱力学場における気温・湿度の条件に大きく左右されると考えられる。乾燥している場合には雨滴に対する蒸発の効果は高く、逆に飽和に近い条件下では蒸発しにくいであろう。日本の暖候期、特に梅雨期には大気下層の相対湿度が 100% に近い状態がしばしば観測され、そのような条件下で大雨がもたらされることが知られている (Tsuguti and Kato 2014; Unuma and Takemi 2016)。そのような大雨をもたらす降水系の内部構造を、Unuma (2024) は世界で初めて現業二重偏波レーダーと地上のディスドロメーターにより捉えることに成功した。二重偏波レーダーの解析から、高度 1.5 km より上層では雨滴の衝突併合、それより下層では衝突分裂のシグナルが捉えられていた (Unuma 2024)。さらに、地上のディスドロメーターデータでも、降水強度が 100 mm h−1^{-1} を超えた時間帯に 2 つのピークを持つ粒形分布が観測され、雨滴の衝突分裂のシグナルが捉えられていることを示した (Unuma 2024)。Unuma (2024) で得られた地上での粒形分布のピーク位置は、数値的に得られる定常的な分布 (Low and List 1982b; McFarquhar 2004a; Straub et al. 2010) のいずれとも異なっていた。この結果を受け、近傍の高層気象観測データを調べたところ、地上付近から対流圏中層付近までほぼ飽和していることが分かった。このことから、非常に湿った環境下では雨滴の蒸発が起こりにくいだけでなく、雨粒になる前の雲粒と雨粒との衝突併合 (accretion) が生じていることでピーク位置が異なっていると考察している (Unuma 2024)。

    以上の研究動向を鑑みると、降水強度が強い場合に限っては周辺大気の環境条件のいくつかを数値実験と同じ条件として見做すことが出来、McFarquhar (2010) で指摘されているように数値実験と観測データとで定常的な分布の一貫性を調べられるのかもしれない。

    おわりに

    複数ピークを持つ雨滴粒径分布について、McFarquhar (2010) を基に理解を深めた。特に、複数のピークを持つ定常的な分布が観測事実として得られるかどうかは決着がついておらず、議論は現在も続いているようである。今後の研究動向に注視するとともに、自身でも研究を進めていきたい。

    参考文献

    更新履歴

    ]]>     @@ -44,7 +44,7 @@ https://takashiunuma.github.io/blog/2024/01/20/ 2024-01-20T11:00:45.000Z - 2024-01-20T13:18:15.693Z + 2024-04-15T14:13:04.493Z はじめに

    前回、Spectral Bin Microphysics Model にて、気象数値モデルを用いた粒形分布の計算例を示した。以前、粒径分布とレーダー変数で示したように、レーダー変数もまた粒形分布から計算可能である。そこで、今回は Spectral Bin Microphysics Model で出力したデータを用い、レーダー変数を計算・描画することを試みる。

    Cloud Resolving Model Radar Simulator

    雲解像モデルを用いたレーダー変数の計算は、いくつかの方法が提案されている。最近では、パッケージ化されたものが開発されてきている。中でも Cloud Resolving Model Radar Simulator (CR-SIM; Oue et al. 2020) は、GNU PUBLIC LICENSE で利用可能なツールの一つである。2020年6月20日時点の最新版は 3.33 となっている。Stony Brook University の Web ページからダウンロードが可能。ただし、ファイルサイズが約 13 GB と少々大きいので、ダウンロード時の回線速度やパケット容量に注意が必要。今回は、レーダー変数の計算に必要なコアパッケージ crsim-3.33.tar.gz のみをダウンロードした。

    プログラムの修正

    tar.gz ファイルを展開後、そのままコンパイルするだけではエラーを吐くので、下記を修正する。

    diff --git a/share/crsim/test/Makefile.in b/share/crsim/test/Makefile.inindex 3b862a3..1fb76f5 100644--- a/share/crsim/test/Makefile.in+++ b/share/crsim/test/Makefile.in@@ -195,21 +195,19 @@ target_alias = @target_alias@ top_build_prefix = @top_build_prefix@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@-TESTS = run_test1.sh run_test2.sh run_test3.sh run_test4.sh+TESTS = run_test1.sh run_test2.sh run_test3.sh EXTRA_DIST = run_test1.sh.in \              run_test2.sh.in \-             run_test3.sh.in \-             run_test4.sh.in+             run_test3.sh.in auxdir = $(datadir)/crsim/test ddir = .$(pwd) edit = sed -e 's|@nccmp[@]|$(NCCMP)|g' nobase_dist_aux_DATA = run_test1.sh \                         run_test2.sh \-                        run_test3.sh \-                        run_test4.sh+                        run_test3.sh-CLEANFILES = run_test1.sh run_test2.sh run_test3.sh run_test4.sh+CLEANFILES = run_test1.sh run_test2.sh run_test3.sh all: all-am .SUFFIXES:@@ -522,21 +520,16 @@ run_test2.sh: run_test2.sh.in Makefile run_test3.sh: run_test3.sh.in Makefile        $(edit) '$(srcdir)/$@.in' > $@        chmod 755 run_test3.sh-run_test4.sh: run_test4.sh.in Makefile-       $(edit) '$(srcdir)/$@.in' > $@-       chmod 755 run_test4.sh install-data-hook:        chmod a+rx $(auxdir)/run_test1.sh        chmod a+rx $(auxdir)/run_test2.sh        chmod a+rx $(auxdir)/run_test3.sh-       chmod a+rx $(auxdir)/run_test4.sh installcheck-local:        cd $(auxdir) && sh run_test1.sh        cd $(auxdir) && sh run_test2.sh        cd $(auxdir) && sh run_test3.sh-       cd $(auxdir) && sh run_test4.sh # Tell versions [3.59,3.63) of GNU make to not export all variables. # Otherwise a system limit (for SysV at least) may be exceeded.diff --git a/src/crsim_subrs.f90 b/src/crsim_subrs.f90index 4539fe4..9c9120b 100644--- a/src/crsim_subrs.f90+++ b/src/crsim_subrs.f90@@ -4023,6 +4023,7 @@       if (conf%MP_PHYSICS==9)  den_string='800' !changed by oue because mpl lut does not includes 840 2017/12/08       if (conf%MP_PHYSICS==10) den_string='500'       if (conf%MP_PHYSICS==8)  den_string='900' !Added by oue 2016/09/16: Note that Thopmson scheme assumes ice density of 890+      if (conf%MP_PHYSICS==20)  den_string='400' !Added for HUCM SBM       if (conf%MP_PHYSICS==30)  den_string='400' !Added by oue 2017/07/17:       if (conf%MP_PHYSICS==40)  den_string='400' !Added by oue 2017/07/21:       if (conf%MP_PHYSICS==50)  den_string='900' !Added by DW for P3

    Makefile.in の run_test4.sh はテスト用のサンプルデータが同梱されておらず実行不可だったので、除外した。また、WRF-SBM の結果を用いる場合、条件分岐に抜けがあったので暫定値を追加した。

    動作環境

    コンパイル

    make check 時に NCCMP が必要となるので、事前にインストールしておく。また、tar.gz ファイルを展開すると、数十 GB 程度に肥大するので、ディスク容量に十分な空きスペースがあることも確認しておく。また、configure の前に FCFLAGS と CFLAGS の環境変数を設定しておく。

    設定出来たら、NetCDF, NetCDF-Fortran, NCCMP, インストール先のそれぞれのディレクトリを指定した上で configure を実行。

    Makefile が作成されたら、make でコンパイル。

    続いて、test case を実行し、結果が nccmp でバイナリ一致することを確認。

    問題なければ、make install を実行する。

    指定したディレクトリ以下に bin/ etc/ share/ がそれぞれあれば OK。

    パラメーターファイル

    PARAMETERS に設定を書き込む仕様となっている。サンプルファイルは、インストール先の etc 以下にある。以下、今回利用するパラメーター例。

    # PARAMETERS CR-SIM v3.2x# speccify the number of OpenMP threads used in the code (integer, OMPThreadNum)1# specify the starting and ending indices of domain and desired saved time step: ix_start,ix_end, iy_start,iy_end, iz_start,iz_end,it. If negative, all points in specific direction are used( 1-> nx, 1-> ny, 1->nz, respectivelly.1,411,411,40    # iz_start,iz_end1       # it# specify the microphysics Morrison 10, JFan 20, Milbranndt-Yau 9 with snow spherical,optionally 901 for Milbranndt-Yau with snow not spherical, P3 50, and Thompson 8 (modified by DW)20# specify the indices of radar position ixc,iyc. if ixc=-999, radar locations are every grid box at iyc, iyc=-999, radar locations are every grid box at ixc.20,20# specify the height of radar in meters1000.d0#Specify the file of ice PSD parameters for P3 microphysics scheme (added by DW)#"./../aux/p3_ice_psd_para_lookup_table.dat" #(added by DW)#assign the scattering type to each hydrometeor (5 lines for MP_PHYSICS=10,20; 6 lines for MP_PHYSICS=9 and 50). Order has to be 1-cloud, 2-rain, 3-ice, 4-snow, 5-graupel and 6 -hail. (For P3, order has to be 1-cloud, 2-rain, 3-small ice, 4-unrimed ice, 5-graupel and 6-partial rimed ice.) (modified by DW)"cloud""rainb""ice_ar0.20""snow_ar0.60""gh_ryzh"#assign the minimum thresholds of the input mixing ratio for the radar simulation. The input values of mixing ratios <= specified threshold will be set to 0. Order has to be 1-cloud, 2-rain, 3-ice, 4-snow, 5-graupel and 6 -hail. (For P3, order has to be 1-cloud, 2-rain, 3-small ice, 4-unrimed ice, 5-graupel and 6-partial rimed ice.) (modified by DW)1.d-81.d-81.d-111.d-111.d-8# horienID (integer),sigma(double precis.) for each hydrometeor category  - choice of orientation distribution (1=chaotic,2=random orientation in the horizontal plane, 3=2D Gaussian with zero mean and sigma standard deviation). In the case those values are negative, the default value are used. If horienID /=3, sigma values are disregarded.Order has to be 1-cloud, 2-rain, 3-ice, 4-snow, 5-graupel (and 6 -hail). (For P3, order has to be 1-cloud, 2-rain, 3-small ice, 4-unrimed ice, 5-graupel and 6-partial rimed ice.) (modified by DW)3,10.d03,10.d03,10.d03,40.d03,40.d0# Specify radar frequency (3.0d0, 5.5d0, 9.5d0, 35.0d0, 94.0d0)5.5d0# Specify elevation, if this ranges from -90 to 90, elevation has a fixed value; if -999, elevation of each pixel is relative to radar origin given by indices (ixc,iyc)0.7d0##turn off the polaraimetric variables : yes ==1 , any other number no0# Specify radar beamwidth ( i.e one-way angular resolution Theta1) in degrees0.1d0# Specify the radar range resolution dr in meters (= (c Tau) / 2 )250.d0# Specify value of coefficient ZMIN in relation dBZ_min(dBZ)=ZMIN(dBZ)+20 log10(range in km)-50.d0# Introduce cloud lidar (ceilometer,fr=905 nm) ceiloID=0 no; =1 yes0# Introduce micropulse lidar (MPL) mplID =1 for 353 nm, mplID =2 for 532 nm, or do not include MPL simulations,mplID<=00# Whether to include an average aerosol profile in MPL simulation (aeroID=1), or not (aeroID=0)2# Value of desired  optical  to normalize the aerosol ext profile, if negative, do not normalize0.1d0# fixed value for aerosol lidar ratio aero_lidar_ratio30.d0# Introduce Doppler spectra (spectraID=1: yes; /=1:  no)0# Post Processing: Produce ARSCL variables : arsclID=1 yes; any other number no. arsclID=1 works with ceiloID=1 and mplID>00# Post Processing: Produce microwave radiometer liquid water path taking account of field of view ARSCL variables : mwrID=1 yes; any other number no.0# Post Processing: Microwave radiometer field of view in degree5.9# Post Processing: Altitude of Microwave radiometer location in m0.0#

    計算領域は、入力に使用する WRF-SBM の mass point でのグリッド数に合わせる。レーダーの設定は、気象庁の C バンドレーダーを想定し、5.5 GHz に変更。Micropulse lidar, Doppler spectra, Microwave radiometer の計算は無しに変更。

    実行

    下記を実行する。引数には、上記で用意したパラメーターファイル、WRF-SBM の出力ファイル、CR-SIM 計算結果の出力ファイル名をそれぞれ指定する。

    ./crsim PARAMETERS wrfout_d01_0001-01-01_00:50:00 wrfout_d01_005000.nc

    実行すると、計算結果の出力ファイル名で、全カテゴリをまとめたもの、雲・雨・雪・氷・雹の各カテゴリに分けたものがそれぞれ出力される。詳細は User Guide を参照されたし。

    データの描画

    今回は、雨 (rain) のデータを可視化してみる。データの読み込みには xarray を用いた。以下、サンプルコード。

    import xarray as xrimport pyart# データを xarray で読み込みds = xr.open_dataset("wrfout_d01_005000_rain.nc")# 水平偏波反射強度 Zhh (dBZ) の水平断面 (k = 0)fig, ax = plt.subplots()ds.Zhh[0,:,:].plot(ax=ax, cmap=pyart.graph.cm_colorblind.HomeyerRainbow, vmin=-25, vmax=65)ax.set_xlim(15,24)ax.set_ylim(15,24)ax.set_title('$Z_{H}$ (k = 0)')ax.grid()# 水平偏波反射強度 ZH (dBZ) の鉛直断面 (j = 22)fig, ax = plt.subplots()ds.Zhh[:,22,:].plot(ax=ax, cmap=pyart.graph.cm_colorblind.HomeyerRainbow, vmin=-25, vmax=65)ax.set_xlim(16,22)ax.set_ylim(0,14)ax.set_title('$Z_{H}$ (j = 22)')ax.grid()# 反射因子差 ZDR (dB) の鉛直断面 (j = 22)fig, ax = plt.subplots()ds.Zdr[:,22,:].plot(ax=ax, cmap=pyart.graph.cm.RefDiff, vmin=-3, vmax=6)ax.set_xlim(16,22)ax.set_ylim(0,14)ax.set_title('$Z_\mathrm{DR}$ (j = 22)')ax.grid()fig, axs = plt.subplots(ncols=2)ds.Zhh[:11,22,19].plot(y='nz', ax=axs[0])axs[0].grid()ds.Zdr[:11,22,19].plot(y='nz', ax=axs[1])axs[1].set_ylabel('')axs[1].grid()fig.suptitle('Vertical profiles of polarimetric retreivals')

    まず、Spectral Bin Microphysics Model で確認した雨水混合比と同じ領域で、計算した水平偏波反射強度の水平分布を確認した。反射強度で見てもだいたい同じ位置で値のピークがある。

    次に、y = 22 に沿った断面を示す。ここでは、水平偏波反射強度と反射因子差を示す。z = 10 がだいたい高度 5 km に相当する。Zhh > 30 dBZ の領域が x = 19 付近で 5 km を超えて分布している。同じ断面を反射因子差で見ると、反射強度の高い領域とは少しずれた部分に ZDR > 3 dB となる分布が計算されている。

    Spectral Bin Microphysics Model の粒径分布の高度変化を確認した grid で水平偏波反射強度と反射因子差の鉛直分布を確認する。z < 4 では、反射強度は地上に向かって減少、反射因子差は地上に向かって微増していることがわかる。ここで、粒径分布とレーダー変数で示した Zhh と ZDR の式の形をよく見ると、Zhh は粒径の 6 乗、ZDR は粒子のアスペクト比にそれぞれ比例する。2 mm 未満の粒径の頻度が地上に向かって減少していたことは、Zhh が地上に向かって減少していることと整合的である。加えて、2 mm 以上の粒径の頻度が地上に向かってやや増加していたことは、ZDR が地上に向かって微増していることに対応していると考えることができる。これを言い換えると、雨粒は大きくなるほど扁平になるので、アスペクト比が大きくなるということを意味する。

    以上のように、粒形分布に加えて二重偏波レーダー変数が算出できると、それぞれの特徴を相互に比較することが可能となる。加えて、観測データと直接比較も容易になる。

    おわりに

    今回は、CR-SIM を用いたレーダー変数の計算、それらの描画について試行した。有益なツールの開発・管理を行っている開発者の方々に謝意を示す。

    参考文献

    更新履歴

    ]]>     @@ -67,9 +67,9 @@ https://takashiunuma.github.io/blog/2024/01/19/ 2024-01-18T21:33:33.000Z - 2024-02-01T01:33:30.051Z + 2024-05-26T05:58:30.897Z - はじめに

    これまでに、粒径分布とレーダー変数や粒径分布パラメータのレーダー変数による推定といった記事で、主に観測データに関連した話題を扱ってきた。観測データはスナップショットではあるものの、真値に近い状態を示す。一方、時間・空間分解能に制約があるために観測結果がどうしてそうなったのか、という要因の特定には不向きな場合がある。その点を補うものとして、数値モデルによる素過程のシミュレーションや解析が挙げられる。中でも雲微物理過程は、雨や雲を扱う重要な物理過程である。計算コストの観点から、粒径分布は特定の関数で近似されたものが使われることが多い。一方、Spectral bin microphysics (SBM) model は、粒径毎の数濃度等を陽に解くものである。最近では、Hebrew University で開発された Cloud Model (Khain et al. 2004, 2008; HUCM) の spectral (bin) microphysics 部分を他の大気モデルへ移植したものが、利用可能になっている (e.g., Khain et al. 2011; Igel and van den Heever 2017)。本稿では、粒径分布の特徴を把握するためのツールの一つとして、WRF-SBM の利用と簡単な出力結果の描画を試みる。

    WRF model

    Weather Research and Frecasting Model (WRF) は、主に米国大気研究センターが開発している領域数値気象予報モデルである。WRF のバージョンは、2023年12月22日時点で 4.5.2 が最新版である (参考)。HUCM SBM は mp_physics = 32 or 30 として実装されている。これらは、それぞれ FULL-SBM、FAST-SBM と呼ばれている。前者は計算時間を要するものの物理過程の厳密さを重視したバージョン、後者は前者を特定の現象に対してやや簡略化し計算コストを減少させたバージョンである。バージョン 4.4.1 以降、FAST-SBM のみ利用可能であり、FULL-SBM は利用不可となっている (参考)。ここでは、計算速度よりも物理過程の正確さを重視することを目的とし、FULL-SBM が利用可能なバージョン 4.2.2 を用いることにした (参考)。

    動作環境

    必要なライブラリは、dnf install した。

    コンパイル

    configure.wrf の LIB_EXTERNAL に -L/usr/lib64 を追記し、コンパイル。

    main/ 以下に ideal.exe と wrf.exe が生成されていれば、コンパイル成功。

    出力変数の追加

    今回は、粒径分布の特徴を把握することを意図し、各 bin の変数を出力できるように修正を施す。出力変数を追加するには、Registry を編集するのが一番簡単。方法は下記。

    定数ファイルの準備

    WRF で SBM を実行するには、追加の定数ファイルが必要となる。下記に従い、定数ファイルをダウンロード・展開・配置する。

    実行

    準備が出来たら、namelist.input を編集し、ideal.exe、wrf.exe の順に実行する。
    今回は、ideal case = 2 の supercell 理想化実験を用いる。namelist.input で mp_physics = 32 に変更し、雲微物理過程のみ用いる。計算領域は、mass point で 41x41x40 grid (水平方向 82 km、鉛直方向に 20 km)。モデル上端から 5 km までダンピング層を配置。ideal.exe を実行すると wrfinput_d01 が出力され、その後 wrf.exe を実行すると wrfout_d01_* が出力される。

    データの描画

    まず、初期プロファイルの input_sounding データを Skew-T LogP plot by MetPy の方法で描画する。ただし、read_fwf の引数に infer_nrows=222 を追加する。理由は、行数のデフォルト値が 100 に制限されているため。

    横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。Skew-T LogP plot by MetPy で示した例よりも赤色の陰影で示される面積が大きく、また、中・上層は乾燥している点も異なる。

    次に、WRF の出力結果 wrfout_d01_* を wrf-python により描画する。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import (to_np, getvar)# WRF 出力ファイル指定ncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# QRAIN 変数を取得qr = getvar(ncfile, "QRAIN")# 描画開始 (xy 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[0,15:25,15:25].west_east)+15,             to_np(qr[0,15:25,15:25].south_north)+15,             to_np(qr[0,15:25,15:25]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title("QRAIN (k = "+str(kk)+")")plt.xlabel('west_east grid number')plt.ylabel('south_north grid number')plt.grid()plt.show()# 描画開始 (xz 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[:15,22,16:23].west_east)+16,             to_np(qr[:15,22,16:23].bottom_top),             to_np(qr[:15,22,16:23]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title(qr.description + " (" + qr.units + ")")plt.xlabel('west_east grid number')plt.ylabel('bottom_top grid number')plt.grid()plt.show()

    計算開始 50 分後の雨水混合比の水平分布。グリッド位置 (19,21) で雨水混合比の極大値が計算されている。

    雨水混合比の y = 22 における x-z 断面。地上付近に計算されていた雨水混合比は、鉛直方向には一様ではなく高度とともに値が変化している。

    以後、QRAIN の極大値が計算されていたグリッド近傍の (19,22) に着目する。

    出力結果から粒径分布を計算

    方針は下記の通り。

    ソースコードは下記。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import getvar# 粒径 (半径) [unit: cm] 定数データ読み込み# index とカテゴリの対応は下記。#  0: aerosol#  1: cloud + rain#  2: ice + snow#  3: graupel + hail#  4: snow#  5: graupel#  6: hailbulkradii = [[0] * 33 for i in range(7)]with open("bulkradii.asc_s_0_03_0_9") as f1:    k = 0    j = 0    for line in f1:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkradii[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 密度 [unit: g cm^-3] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。bulkdens = [[0] * 33 for i in range(7)]with open("bulkdens.asc_s_0_03_0_9") as f2:    k = 0    j = 0    for line in f2:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkdens[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 質量 [unit: g] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。masses = [[0] * 33 for i in range(7)]with open("masses.asc") as f3:    k = 0    j = 0    for line in f3:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                masses[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# WRF 出力結果を読み込みncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# 出力結果の cloud and rain bin データを xarray dataset に変換。vlist = ['ff1i01', 'ff1i02', 'ff1i03', 'ff1i04', 'ff1i05',         'ff1i06', 'ff1i07', 'ff1i08', 'ff1i09', 'ff1i10',         'ff1i11', 'ff1i12', 'ff1i13', 'ff1i14', 'ff1i15',         'ff1i16', 'ff1i17', 'ff1i18', 'ff1i19', 'ff1i20',         'ff1i21', 'ff1i22', 'ff1i23', 'ff1i24', 'ff1i25',         'ff1i26', 'ff1i27', 'ff1i28', 'ff1i29', 'ff1i30',         'ff1i31', 'ff1i32', 'ff1i33']xdsff1i = getvar(ncfile, "ff1i01")xdsff1i = xdsff1i.to_dataset()for ivar in vlist:    xdsff1i[ivar] = getvar(ncfile, ivar)# 粒径を半径から直径に、単位を cm から mm に、それぞれ変換。diameter = np.zeros((33), float)for i in range(0,33):    diameter[i] = bulkradii[0][i] * 2. * 10. # unit: diameter in mm# 特定グリッドの指定 (水平分布で極大値のあったグリッド)jj = 22ii = 19for kk in [3, 2, 1, 0]:    # 粒径分布データ用の numpy array を初期化。    ff1i = np.zeros((33), float)    # 33 bin に分かれた変数を一つの numpy array に結合    i = 0    for ivar in vlist:        ff1i[i] = xdsff1i[ivar][kk,jj,ii].to_numpy()        i = i + 1    # 粒形分布 N(D) [m^-3 mm^-1] 計算    for i in range(0,33):        # unit: (kg kg^-1) (g^-1) (g cm^-3) (mm^-1)        #    -> (kg kg^-1) (kg^-1) (kg m^-3) (mm^-1)        #    -> m^-3 mm^-1 (多分これであっているはず)        ff1i[i] = (ff1i[i] / (masses[0][i] * 0.001)) * (bulkdens[0][i] * 1000.) / (diameter[i] * 1000.)    # 粒径分布データ描画。ただし、rain bin のみ。    plt.plot(diameter[16:33], ff1i[16:33], label="k = "+str(kk))    plt.scatter(diameter[16:33], ff1i[16:33])plt.yscale('log')plt.xlim(0,6)plt.ylim(0.01, 10000)plt.title('N(D) for rain at ('+str(ii)+','+str(jj)+')')plt.xlabel('Diameter (mm)')plt.ylabel('Number concentration (mm$^{-1}$ m$^{-3}$)')plt.grid()plt.legend()plt.show()

    粒形分布を確認すると、地上に近い (i.e., k = 0) ほど粒径 2 mm 未満の頻度が低く、上空 (i.e., k = 3) ほど頻度が高くなっている。一方、粒径 2 mm よりも大きい粒径は地上ほど頻度が高い。また、頻度のピークは、地上に近いほど複数 (粒径 0.26 mm, 0.51 mm, 1.62 mm) 現れる傾向にある。

    鉛直分布を確認した時、雨水混合比の値は k = 8 を極大とし、地上に向かって減少していた。これは、粒形分布では第一近似として粒径 2 mm 未満の頻度が低くなっていたことに対応していると考えられる。雨水混合比としては減少する一方、粒径 2 mm 以上の頻度が僅かながら高くなっていた。このことは、雨滴同士の衝突併合によって雨粒が地上に向かって成長している可能性を示唆する。このように、粒形分布を確認することで雨の特徴をより多角的に捉えることが可能となる。

    おわりに

    今回は、WRF-SBM を用いた粒形分布の解析例を示した。実は、十数年前に WRF モデルの version 3 系列を利用したことがあった。その当時よりも利用可能なオプションが格段に増えており、使い勝手もとてもよくなっていると感じた。特に、描画を Python で行える点は、近年の流行を追っていて舌を巻く。当時は大雨を再現するだけでなく出力結果の解析にも時間と労力を費やしていたが、今後はそんな苦労も少なくなるだろうか。

    参考文献

    更新履歴

    ]]>     + はじめに

    これまでに、粒径分布とレーダー変数や粒径分布パラメータのレーダー変数による推定といった記事で、主に観測データに関連した話題を扱ってきた。観測データはスナップショットではあるものの、真値に近い状態を示す。一方、時間・空間分解能に制約があるために観測結果がどうしてそうなったのか、という要因の特定には不向きな場合がある。その点を補うものとして、数値モデルによる素過程のシミュレーションや解析が挙げられる。中でも雲微物理過程は、雨や雲を扱う重要な物理過程である。計算コストの観点から、粒径分布は特定の関数で近似されたものが使われることが多い。一方、Spectral bin microphysics (SBM) model は、粒径毎の数濃度等を陽に解くものである。最近では、Hebrew University で開発された Cloud Model (Khain et al. 2004, 2008; HUCM) の spectral (bin) microphysics 部分を他の大気モデルへ移植したものが、利用可能になっている (e.g., Khain et al. 2011; Igel and van den Heever 2017)。本稿では、粒径分布の特徴を把握するためのツールの一つとして、WRF-SBM の利用と簡単な出力結果の描画を試みる。

    WRF model

    Weather Research and Frecasting Model (WRF) は、主に米国大気研究センターが開発している領域数値気象予報モデルである。WRF のバージョンは、2023年12月22日時点で 4.5.2 が最新版である (参考)。HUCM SBM は mp_physics = 32 or 30 として実装されている。これらは、それぞれ FULL-SBM、FAST-SBM と呼ばれている。前者は計算時間を要するものの物理過程の厳密さを重視したバージョン、後者は前者を特定の現象に対してやや簡略化し計算コストを減少させたバージョンである。バージョン 4.4.1 以降、FAST-SBM のみ利用可能であり、FULL-SBM は利用不可となっている (参考)。ここでは、計算速度よりも物理過程の正確さを重視することを目的とし、FULL-SBM が利用可能なバージョン 4.2.2 を用いることにした (参考)。

    動作環境

    必要なライブラリは、dnf install した。

    コンパイル

    configure.wrf の LIB_EXTERNAL に -L/usr/lib64 を追記し、コンパイル。

    main/ 以下に ideal.exe と wrf.exe が生成されていれば、コンパイル成功。

    出力変数の追加

    今回は、粒径分布の特徴を把握することを意図し、各 bin の変数を出力できるように修正を施す。出力変数を追加するには、Registry を編集するのが一番簡単。方法は下記。

    定数ファイルの準備

    WRF で SBM を実行するには、追加の定数ファイルが必要となる。下記に従い、定数ファイルをダウンロード・展開・配置する。

    実行

    準備が出来たら、namelist.input を編集し、ideal.exe、wrf.exe の順に実行する。
    今回は、ideal case = 2 の supercell 理想化実験を用いる。namelist.input で mp_physics = 32 に変更し、雲微物理過程のみ用いる。計算領域は、mass point で 41x41x40 grid (水平方向 82 km、鉛直方向に 20 km)。モデル上端から 5 km までダンピング層を配置。ideal.exe を実行すると wrfinput_d01 が出力され、その後 wrf.exe を実行すると wrfout_d01_* が出力される。

    データの描画

    まず、初期プロファイルの input_sounding データを Skew-T LogP plot by MetPy の方法で描画する。ただし、read_fwf の引数に infer_nrows=222 を追加する。理由は、行数のデフォルト値が 100 に制限されているため。

    横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。Skew-T LogP plot by MetPy で示した例よりも赤色の陰影で示される面積が大きく、また、中・上層は乾燥している点も異なる。

    次に、WRF の出力結果 wrfout_d01_* を wrf-python により描画する。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import (to_np, getvar)# WRF 出力ファイル指定ncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# QRAIN 変数を取得qr = getvar(ncfile, "QRAIN")# 描画開始 (xy 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[0,15:25,15:25].west_east)+15,             to_np(qr[0,15:25,15:25].south_north)+15,             to_np(qr[0,15:25,15:25]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title("QRAIN (k = "+str(kk)+")")plt.xlabel('west_east grid number')plt.ylabel('south_north grid number')plt.grid()plt.show()# 描画開始 (xz 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[:15,22,16:23].west_east)+16,             to_np(qr[:15,22,16:23].bottom_top),             to_np(qr[:15,22,16:23]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title(qr.description + " (" + qr.units + ")")plt.xlabel('west_east grid number')plt.ylabel('bottom_top grid number')plt.grid()plt.show()

    計算開始 50 分後の雨水混合比の水平分布。グリッド位置 (19,21) で雨水混合比の極大値が計算されている。

    雨水混合比の y = 22 における x-z 断面。地上付近に計算されていた雨水混合比は、鉛直方向には一様ではなく高度とともに値が変化している。

    以後、QRAIN の極大値が計算されていたグリッド近傍の (19,22) に着目する。

    出力結果から粒径分布を計算

    方針は下記の通り。

    ソースコードは下記。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import getvar# 粒径 (半径) [unit: cm] 定数データ読み込み# index とカテゴリの対応は下記。#  0: cloud + rain#  1: ice/columns#  2: ice/plates#  3: ice/dendrites#  4: snow#  5: graupel#  6: hailbulkradii = [[0] * 33 for i in range(7)]with open("bulkradii.asc_s_0_03_0_9") as f1:    k = 0    j = 0    for line in f1:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkradii[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 密度 [unit: g cm^-3] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。bulkdens = [[0] * 33 for i in range(7)]with open("bulkdens.asc_s_0_03_0_9") as f2:    k = 0    j = 0    for line in f2:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkdens[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 質量 [unit: g] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。masses = [[0] * 33 for i in range(7)]with open("masses.asc") as f3:    k = 0    j = 0    for line in f3:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                masses[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# WRF 出力結果を読み込みncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# 出力結果の cloud and rain bin データを xarray dataset に変換。vlist = ['ff1i01', 'ff1i02', 'ff1i03', 'ff1i04', 'ff1i05',         'ff1i06', 'ff1i07', 'ff1i08', 'ff1i09', 'ff1i10',         'ff1i11', 'ff1i12', 'ff1i13', 'ff1i14', 'ff1i15',         'ff1i16', 'ff1i17', 'ff1i18', 'ff1i19', 'ff1i20',         'ff1i21', 'ff1i22', 'ff1i23', 'ff1i24', 'ff1i25',         'ff1i26', 'ff1i27', 'ff1i28', 'ff1i29', 'ff1i30',         'ff1i31', 'ff1i32', 'ff1i33']xdsff1i = getvar(ncfile, "ff1i01")xdsff1i = xdsff1i.to_dataset()for ivar in vlist:    xdsff1i[ivar] = getvar(ncfile, ivar)# 粒径を半径から直径に、単位を cm から mm に、それぞれ変換。diameter = np.zeros((33), float)for i in range(0,33):    diameter[i] = bulkradii[0][i] * 2. * 10. # unit: diameter in mm# 特定グリッドの指定 (水平分布で極大値のあったグリッド)jj = 22ii = 19for kk in [3, 2, 1, 0]:    # 粒径分布データ用の numpy array を初期化。    ff1i = np.zeros((33), float)    # 33 bin に分かれた変数を一つの numpy array に結合    i = 0    for ivar in vlist:        ff1i[i] = xdsff1i[ivar][kk,jj,ii].to_numpy()        i = i + 1    # 粒形分布 N(D) [m^-3 mm^-1] 計算    for i in range(0,33):        # unit: (kg kg^-1) (g^-1) (g cm^-3) (mm^-1)        #    -> (kg kg^-1) (kg^-1) (kg m^-3) (mm^-1)        #    -> m^-3 mm^-1 (多分これであっているはず)        ff1i[i] = (ff1i[i] / (masses[0][i] * 0.001)) * (bulkdens[0][i] * 1000.) / (diameter[i] * 1000.)    # 粒径分布データ描画。ただし、rain bin のみ。    plt.plot(diameter[16:33], ff1i[16:33], label="k = "+str(kk))    plt.scatter(diameter[16:33], ff1i[16:33])plt.yscale('log')plt.xlim(0,6)plt.ylim(0.01, 10000)plt.title('N(D) for rain at ('+str(ii)+','+str(jj)+')')plt.xlabel('Diameter (mm)')plt.ylabel('Number concentration (mm$^{-1}$ m$^{-3}$)')plt.grid()plt.legend()plt.show()

    粒形分布を確認すると、地上に近い (i.e., k = 0) ほど粒径 2 mm 未満の頻度が低く、上空 (i.e., k = 3) ほど頻度が高くなっている。一方、粒径 2 mm よりも大きい粒径は地上ほど頻度が高い。また、頻度のピークは、地上に近いほど複数 (粒径 0.26 mm, 0.51 mm, 1.62 mm) 現れる傾向にある。

    鉛直分布を確認した時、雨水混合比の値は k = 8 を極大とし、地上に向かって減少していた。これは、粒形分布では第一近似として粒径 2 mm 未満の頻度が低くなっていたことに対応していると考えられる。雨水混合比としては減少する一方、粒径 2 mm 以上の頻度が僅かながら高くなっていた。このことは、雨滴同士の衝突併合によって雨粒が地上に向かって成長している可能性を示唆する。このように、粒形分布を確認することで雨の特徴をより多角的に捉えることが可能となる。

    おわりに

    今回は、WRF-SBM を用いた粒形分布の解析例を示した。実は、十数年前に WRF モデルの version 3 系列を利用したことがあった。その当時よりも利用可能なオプションが格段に増えており、使い勝手もとてもよくなっていると感じた。特に、描画を Python で行える点は、近年の流行を追っていて舌を巻く。当時は大雨を再現するだけでなく出力結果の解析にも時間と労力を費やしていたが、今後はそんな苦労も少なくなるだろうか。

    参考文献

    更新履歴

    ]]>     @@ -90,7 +90,7 @@ https://takashiunuma.github.io/blog/2024/01/02/ 2024-01-01T15:00:08.000Z - 2024-01-01T03:12:10.778Z + 2024-04-15T14:13:04.496Z はじめに

    学術論文を読む際、タイトルやアブストラクトを一読し、どんな内容か把握することが多い。学術誌によっては、キーワードが書いてあったりするので、それを参考にするのも良いだろう。また、文章内に登場する単語毎の頻度分布を作成し、その頻出する語によって文章の特徴を把握することも可能である。最近では、WordCloud というものが Python ベースで利用可能となってきている。それを応用し、大気科学系の学術誌が出版されたら、URL を入力として論文内の単語を解析、その単語の頻度分布を図化したものとともに X に投稿するといった取り組みをしている方も居るようだ。ただし、XML や plain HTML を入力の想定としている。このため、既に持っている PDF 形式の文章や論文では解析が出来ないという課題があった。そこで、本稿では PDF 形式の文章でも行える WordCloud によるテキストマイニングについて紹介する。

    準備

    ソースコード

    基本的には、こちらを元にした。本家と違うところは、除外する単語のリストを読み込んでいるところと、PDF 形式の文章を pdftotext で TEXT 形式にしているところである。

    import osimport refrom glob import globimport pdftotextfrom wordcloud import WordCloud, STOPWORDSimport matplotlib.pyplot as plt# 除外する単語のリストの設定exclude_words = list(STOPWORDS)stopword_files = glob(os.path.join('data/', 'exclude_words/', "*.txt"))for fname in stopword_files:    with open(fname, "r") as f:        exclude_words += f.read().split("\n")exclude_words = set(exclude_words)# PDF データの取得 (例として Unuma et al. 2023 をダウンロード)!wget 'https://www.jstage.jst.go.jp/article/sola/19/0/19_2023-020/_pdf/-char/en' -O test.pdf# PDF を読み込み、text 整形。with open("test.pdf", "rb") as f:    pdf = pdftotext.PDF(f)text = []for page in pdf:    text.append(re.sub('[0-9][0-9]', ' ', page).replace("\n", " "))text = ''.join(text)text = re.sub('[0-9]', ' ', text).replace('      ', " ").replace('     ', " ").replace('    ', " ").replace('   ', " ")# WordCloud による解析・可視化wordcloud = WordCloud(stopwords=exclude_words, colormap='cividis').generate(text)plt.imshow(wordcloud, interpolation='bilinear')plt.axis("off")plt.savefig('skewt.png', bbox_inches='tight')

    出力結果

    図のような結果となる。使っている単語の数が多いものほど単語が大きく表示され、少ないものほど小さくなる。

    本当は raindrop size distribution で解析してほしかったのだが、raindrop size と size distribution で分割されてしまっている。ほか、著者名は除外リストに入れたほうが良いかもしれない。

    おわりに

    PDF 形式の文章でも、Wordcloud を用いてテキストマイニングをする方法について紹介した。

    参考文献/URL

    更新履歴

    ]]>     @@ -113,7 +113,7 @@ https://takashiunuma.github.io/blog/2024/01/01/ 2024-01-01T04:15:24.000Z - 2024-01-01T05:37:46.000Z + 2024-01-01T05:37:46.502Z はじめに

    大気の鉛直構造を把握するため、日々、高層気象観測が行われている。これらは観測としてのみならず、数値予報モデルの初期値・解析値としても用いられる重要な観測データである。本稿では、その観測データを地点毎に把握するための Skew-T LogP plot について述べる。ここでは、MetPy で提供されているものを利用し、特に理想化した数値実験の初期値に用いる一点鉛直プロファイルの可視化について紹介する。

    データ

    とある時刻・地点のデータを利用する。
    1行目は、地上の気圧 (hPa)・温位 (K)・水蒸気混合比 (g/kg)。
    2行目以降は、高度 (m)・温位 (K)・水蒸気混合比 (g/kg)・東西風 (m/s)・南北風 (m/s)。
    これを、input_sounding という名前で保存。

     1000.00  298.81   18.08  277.04  299.83   16.13   -5.29    1.54  504.73  299.94   16.02   -5.96    1.84  736.79  300.06   16.00   -6.08    2.16  973.81  301.25   15.48   -5.57    3.11 1465.20  303.95   14.38   -2.57    5.48 1982.15  306.88   12.90    1.13    5.17 3104.82  312.11    9.51    5.20    4.45 4369.76  318.33    6.41    5.66    4.46 5823.29  325.50    4.00   11.84   -0.45 7542.22  334.36    1.78   15.32    3.46 9655.91  343.85    0.65   15.49   -5.9212000.00  351.53    0.11   15.49   -5.9221000.00  500.53    0.01   15.49   -5.92

    ソースコード

    上記のデータを、下記で可視化する。

    import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport matplotlibimport matplotlib.dates as mdatesimport metpy.calc as mpcalcfrom metpy.calc import temperature_from_potential_temperature, height_to_pressure_std, dewpoint_from_specific_humidityfrom metpy.plots import SkewTfrom metpy.units import units, pandas_dataframe_to_unit_arraysimport pint_xarray# ヘッダーを取り除いた上でカラムを指定、pandas dataframe として格納。df = pd.read_fwf('input_sounding', header=0, usecols=[0, 1, 2, 3, 4], names=['height', 'theta', 'qv', 'uwnd', 'vwnd'])# pandas dataframe を xarray dataset に変換ds = df.to_xarray()# MetPy.units を用い、単位を指定。ds['height'] = ds['height'] * units.meterds['theta'] = ds['theta'] * units.kelvinds['qv'] = ds['qv'] * units('g/kg')ds['uwnd'] = ds['uwnd'] * units.meter_per_secondds['vwnd'] = ds['vwnd'] * units.meter_per_second# SkewT logP plot の描画に必要な変数に変換。ds['pressure'] = height_to_pressure_std(ds['height']).pint.to('hPa')ds['temperature'] = temperature_from_potential_temperature(ds['pressure'], ds['theta']).pint.to('degC')ds['dewpoint'] = dewpoint_from_specific_humidity(ds['pressure'], ds['temperature'], ds['qv']).pint.to('degC')# ある高度から持ち上げたパーセルの気温・露点温度のプロファイルを計算。parcel_prof = mpcalc.parcel_profile(ds['pressure'], ds['temperature'][0], ds['dewpoint'][0]).pint.to('degC')# 持ち上げ凝結高度の気圧・気温・露点温度を計算。lcl_pressure, lcl_temperature = mpcalc.lcl(ds['pressure'][0], ds['temperature'][0], ds['dewpoint'][0])# 描画fig = plt.figure(figsize=(8, 6))skew = SkewT(fig)skew.ax.set_ylim(1050,100)skew.ax.set_xlim(-30,35)skew.ax.set_title("A sample sounding")skew.plot(ds['pressure'], ds['temperature'], color="k", linestyle='-')skew.plot(ds['pressure'], ds['dewpoint'], linestyle='-')skew.shade_cape(ds['pressure'], ds['temperature'], parcel_prof)skew.plot(ds['pressure'], parcel_prof, linestyle=':', linewidth=2)skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')skew.plot_moist_adiabats(linewidth=0.5)fig.savefig('skewt.png', bbox_inches='tight')# 各種指数・パラメーターの計算p = ds['pressure']T = ds['temperature']Td = ds['dewpoint']height = ds['height']u = ds['uwnd']v = ds['vwnd']ctotals = mpcalc.cross_totals(p, T, Td)kindex = mpcalc.k_index(p, T, Td)showalter = mpcalc.showalter_index(p, T, Td)total_totals = mpcalc.total_totals_index(p, T, Td)vert_totals = mpcalc.vertical_totals(p, T)lift_index = mpcalc.lifted_index(p, T, parcel_prof)cape, cin = mpcalc.cape_cin(p, T, Td, parcel_prof)lclp, lclt = mpcalc.lcl(p[0], T[0], Td[0])lfcp, _ = mpcalc.lfc(p, T, Td)el_pressure, _ = mpcalc.el(p, T, Td, parcel_prof)ml_t, ml_td = mpcalc.mixed_layer(p, T, Td, depth=50 * units.hPa)ml_p, _, _ = mpcalc.mixed_parcel(p, T, Td, depth=50 * units.hPa)mlcape, mlcin = mpcalc.mixed_layer_cape_cin(p, T, parcel_prof, depth=50 * units.hPa)mu_p, mu_t, mu_td, _ = mpcalc.most_unstable_parcel(p, T, Td, depth=50 * units.hPa)mucape, mucin = mpcalc.most_unstable_cape_cin(p, T, Td, depth=50 * units.hPa)(u_storm, v_storm), *_ = mpcalc.bunkers_storm_motion(p, u, v, height)critical_angle = mpcalc.critical_angle(p, u, v, height, u_storm, v_storm)# 各種指数・パラメーターの値を確認print(f'        CAPE: {cape:.2f}')print(f'         CIN: {cin:.2f}')print(f'LCL Pressure: {lclp:.2f}')print(f'LFC Pressure: {lfcp:.2f}')print(f' EL Pressure: {el_pressure:.2f}')print()print('Mixed Layer - Lowest 50-hPa')print(f'     ML Temp: {ml_t:.2f}')print(f'     ML Dewp: {ml_td:.2f}')print(f'     ML CAPE: {mlcape:.2f}')print(f'      ML CIN: {mlcin:.2f}')print()print('Most Unstable - Lowest 50-hPa')print(f'     MU Temp: {mu_t:.2f}')print(f'     MU Dewp: {mu_td:.2f}')print(f' MU Pressure: {mu_p:.2f}')print(f'     MU CAPE: {mucape:.2f}')print(f'      MU CIN: {mucin:.2f}')print()print(f'   Lifted Index: {lift_index:.2f}')print(f'        K-Index: {kindex:.2f}')print(f'Showalter Index: {showalter:.2f}')print(f'   Cross Totals: {ctotals:.2f}')print(f'   Total Totals: {total_totals:.2f}')print(f'Vertical Totals: {vert_totals:.2f}')print()print('Bunkers Storm Motion Vector')print(f'  u_storm: {u_storm:.2f}')print(f'  v_storm: {v_storm:.2f}')print(f'Critical Angle: {critical_angle:.2f}')

    出力結果

    図が Skew-T logP plot。横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。
    図の左下から右上に向かう灰色の線は等温線を、図の右下から左上に伸びる青色点線は等湿潤断熱線をそれぞれ示す。

    図の 925 hPa から少なくとも 500 hPa までは、気温と露点温度の線が重なっている。これは、ほぼ相対湿度 100% であることを示す。加えて、800 hPa から 600 hPa までの層では、気温減率の傾きが湿潤断熱線の傾きよりも大きくなっており、湿潤絶対不安定 (Moist Absolutely Unstable Layer, Bryan and Fritsch 2000) の可能性がある。

    下記は、一点鉛直プロファイルから得られる指数やパラメーターの出力である。一地点の鉛直プロファイルから、これだけ多様な診断パラメーターを算出することが可能となる。

            CAPE: 729.54 joule / kilogram         CIN: -3.21 joule / kilogramLCL Pressure: 926.10 hectopascalLFC Pressure: 885.54 hectopascal EL Pressure: 260.51 hectopascalMixed Layer - Lowest 50-hPa     ML Temp: 22.89 degree_Celsius     ML Dewp: 20.64 degree_Celsius     ML CAPE: 1997.43 joule / kilogram      ML CIN: -0.15 joule / kilogramMost Unstable - Lowest 50-hPa     MU Temp: 20.55 degree_Celsius     MU Dewp: 20.17 degree_Celsius MU Pressure: 927.77 hectopascal     MU CAPE: 754.89 joule / kilogram      MU CIN: -1.08 joule / kilogram   Lifted Index: -3.48 delta_degree_Celsius        K-Index: 41.32 degree_CelsiusShowalter Index: -3.25 delta_degree_Celsius   Cross Totals: 24.40 delta_degree_Celsius   Total Totals: 48.64 delta_degree_CelsiusVertical Totals: 24.24 delta_degree_CelsiusBunkers Storm Motion Vector  u_storm: 0.25 meter_per_second  v_storm: -4.63 meter_per_secondCritical Angle: 179.85 degree

    指数・診断パラメーターは、発案された地域や特性に依存するため、計算結果をそのまま異なる地域に当てはめることは得策ではない。上記で調べたように、可視化することで指数が意図した計算となっているか、愚直に確認することが重要である (時に時間を要するが)。ここでは、計算方法の紹介に留める。

    おわりに

    大気の状態を把握するための基礎的な観測の一つである、高層気象観測データの可視化と関連する指数・パラメータの計算について紹介した。

    参考文献

    更新履歴

    ]]>     @@ -136,7 +136,7 @@ https://takashiunuma.github.io/blog/2022/11/19/ 2022-11-19T07:01:44.000Z - 2024-01-20T10:18:58.957Z + 2024-04-15T14:13:04.504Z はじめに粒径分布とレーダー変数の記事では、粒径分布からどのような仮定の下でレーダー変数を計算しているかについてまとめた。今回はその逆で、レーダー変数から如何にして粒径分布パラメータを推定するか、について述べる。

    粒径分布

    粒径分布は、対象とする粒径の範囲によって当てはめうる関数が異なる。ここでは、地上に到達する範囲の粒径、すなわち約 0.5 mm 以上の雨滴を対象に考える。多くの場合、ガンマ分布で粒径分布を近似しうることが知られている。ガンマ分布は次式で与えられる (Ulbrich 1983)。

    N(D)=N0Dμexp⁡(−ΛD)\quad N(D) = N_{\rm 0} D^{\mu} \exp(-\Lambda D)

    ここで、N0N_{\rm 0}, μ\mu, Λ\Lambdaは、それぞれ Number of concentration, Shape parameter, Slope parameter でありガンマ分布を特徴づけるパラメータである。しかしながら、これらのパラメータは物理的な意味合いを持たないために、気象学的にも直感的には分かりにくいという課題がある。例えば、N0N_{0} の単位は mm−1^{-1} m−(1+μ)^{-(1 + \mu)} で次元が μ\mu に依存するため一定ではない、等である。

    そこで、このガンマ分布を規格化したものとして、次式が提案されている。

    N(D)=Nwf(μ)(DD0)μexp⁡[−(3.67+μ)DD0]\quad N(D) = N_{\rm w} f(\mu) \left( \cfrac{D}{D_{0}} \right)^{\mu} \exp\left[ -(3.67 + \mu) \cfrac{D}{D_{0}} \right]

    ここで、

    Nw=(3.67)4πρw(103WD04)\quad N_{\rm w} = \cfrac{(3.67)^{4}}{\pi \rho_{\rm w}} \left( \cfrac{10^{3} W}{D_{0}^{4}} \right)

    は Intercept parameter と呼ばれ、規格化した数濃度 (mm−1^{-1} m−3^{-3}) と見做せる。この場合、N0N_{0} と異なり単位の次元が一定となる。D0D_{0} は中央粒径 (mm)、ρw\rho_{\rm w} は雨滴の密度 (kg m−3^{-3})、WW は Rainwater content (g m−3^{-3}) である。

    f(μ)=6(3.67)4(3.67+μ)μ+4Γ(μ+4)\quad f(\mu) = \cfrac{6}{(3.67)^{4}} \cfrac{(3.67 + \mu)^{\mu + 4}}{\Gamma(\mu + 4)}

    は、μ\mu (mm−1^{-1}) に依存する関数である。DD は粒径 (mm) を示す。また、

    Dm=(4+μ3.67+μ)D0\quad D_{\rm m} = \left( \cfrac{4 + \mu}{3.67 + \mu} \right) D_{0}

    の関係を用いると、質量平均した粒径 DmD_{\rm m} を用いて規格化したガンマ分布を表現することも可能である。この時、(3.67+μ)(3.67 + \mu) の部分は、ΛD0≅3.67+μ\Lambda D_{0} \cong 3.67 + \mu である。これは、D0D_{0} と DmaxD_{\rm max} との関係が Dmax/D0>2.5D_{\rm max}/D_{0} > 2.5 となる場合に ΛDmax\Lambda D_{\rm max} と ΛD0\Lambda D_{0} との関係がほぼ一定値となることを利用した近似である (Ulbrich 1983, Fig. 1; Sekhon and Srivastava 1970)。

    このように、規格化したガンマ分布の場合、多少近似が介入するものの、粒径分布を直感的分かりやすいパラメータに置き換えられるという利点がある。

    粒径分布パラメータの推定

    粒径分布とレーダー変数で述べたように、気象レーダーではレイリー近似による仮定が前提となる。気象レーダーで得られる信号から直接粒径分布を得ることは出来ないものの、推定する粒径分布の関数や束縛条件を与えることで粒径分布パラメータを推定する手法が考案されてきた。代表的な手法として、β\beta method (Gorgucci et al. 2000, 2001, 2002) が知られている。以下、β\beta method について概観する。

    粒径分布として、前述の規格化したガンマ分布を仮定する。Raindrop shape を変数 β\beta として扱い (名前の由来)、ZhZ_{\rm h}, KDPK_{\rm DP}, ZDRZ_{\rm DR} から次式で β\beta を求める (Gorgucci et al. 2000, eq. A1)。

    β=aZhbKDPc×10−dZDR\quad \beta = a Z_{\rm h}^{b} K_{\rm DP}^{c} \times 10^{-d Z_{\rm DR}}

    ここで、a, b, c, d は定数、ZhZ_{\rm h} (mm6^{6} m−3^{-3}) は水平偏波反射強度の linear unit (真数表記)、ZDRZ_{\rm DR} (dB) は反射因子差の logarithm unit (対数表記) である (以下、同様に扱う)。求めた β\beta と ZhZ_{\rm h}, ZDRZ_{\rm DR} から D0D_{0} NwN_{\rm w} μ\mu をそれぞれ推定する (Gorgucci et al. 2002, eq (24), (34), and (39))。

    D0=a1Zhb1(100.1ZDR)c1log⁡10Nw=a2Zhb2(100.1ZDR)c2μ=a3Zhb3(100.1ZDR)c3\quad D_{0} = a_{1} Z_{\rm h}^{b_{1}} (10^{0.1 Z_{\rm DR}})^{c_{1}} \\\quad \log_{10} N_{\rm w} = a_{2} Z_{\rm h}^{b_{2}} (10^{0.1 Z_{\rm DR}})^{c_{2}} \\\quad \mu = a_{3} Z_{\rm h}^{b_{3}} (10^{0.1 Z_{\rm DR}})^{c_{3}}

    ここで、a1_{1}, b1_{1}, c1_{1}, a2_{2}, b2_{2}, c2_{2} は定数、a3_{3}, b3_{3}, a3_{3} はβ\beta から求まる定数である。具体的な値は、散乱計算から求められる。この点が β\beta method の一つの利点と言える。各定数の S バンドレーダーでの具体的な値は、Gorgucci et al. (2002) を参照されたし。

    Gorgucci et al. (2002) では、粒径分布パラメータの推定の式を ZhZ_{\rm h} から KDPK_{\rm DP} 変えた場合に、得られる粒径分布パラメータの推定精度への影響も調べられている。さらに、Gorgucci and Baldini (2009) では、S, C, X バンドレーダーを想定したシミュレーションを行い、各送信周波数毎に異なる β\beta が与えられた場合の降水強度の推定精度が調べられている。以上のように、散乱計算で係数を求めさえすれば何らかの粒径分布パラメータの値や降水強度の値が算出可能となる。ただし、どの論文でも必ずと言っていいほど地上観測との比較が行われており、推定精度の確認は必須と言える。

    他には、Exponential method (Seliga and Bringi, 1976)、Normalized Gamma method (Testud et al. 2001)、Constrained Gamma method (Zhang et al. 2001)、偏波間相関係数を用いた粒径分布パラメータの推定手法 (Thurai et al. 2008)、Double-moment model (Anagnostou et al. 2009)、Self-consistent with optical parameterization attenuation correction and microphysics estimation (Raupach and Beme 2017)、Inverse model (Alcoba et al. 2022) 等の手法が提案されている。これらの手法については紹介のみに留め、詳細は省略する。

    おわりに

    レーダー変数から粒径分布パラメータを推定する、様々な手法と考え方について述べた。歴史的な背景としては、S バンドレーダーを用いた推定手法が考案され、その後 C バンドレーダーや X バンドレーダーに応用されているようである。恐らく、送信周波数帯毎に降雨減衰の度合いが異なるため、それらの影響が小さいものから検討が進められたと考えられる。加えて、推定結果はそれだけに留まらず、観測結果との比較検証を通して推定精度を担保し、結果を吟味することの重要性が暗に示唆された。

    参考文献

    更新履歴

    ]]>     @@ -204,7 +204,7 @@ https://takashiunuma.github.io/blog/2022/06/18/ 2022-06-18T01:15:10.000Z - 2024-01-20T10:07:47.797Z + 2024-04-15T14:13:04.504Z はじめに二重偏波レーダー変数の算出とその補正では、降水強度算出に関連するレーダー変数の算出方法とレーダー変数推定時にバイアスとなりうる要素に対する補正方法について具体例を示した。この中で降雨減衰補正を利用したものの、その内容を十分に把握できていなかった。そこで、本稿では降雨減衰補正 (主に ZPHI method) の内容を把握することを目的とする。

    ZPHI method の変遷

    降水によるマイクロ波放射の減衰は、反射強度 ZZ と 反射因子差 ZDRZ_{\rm DR} の観測時 (特に C や X などの短い波長帯で) に顕著なバイアスとなりうる。降雨による減衰 (以後、降雨減衰) は、反射強度 ZZ を用いた Hitschfeld and Bordan (1954) の解として単偏波レーダーの時代から推定されてきた。しかしながら、減衰が大きい場合には Hitschfeld and Bordan (1954) で得られる解は、非常に不安定となる。

    そこで、解を安定させるための制約条件として、独立した観測量により推定される path-integrated attenuation (PIA) が注目されるようになった。PIA が既知ならば、減衰量の視線方向のプロファイルが得られる (Meneghini and Nakamura 1990)。PIA は、山岳や背の高いタワー、海面といった比較的安定した物体を対象に用いることで、航空機や衛星観測において推定可能である。

    二重偏波要素が観測可能な場合は、偏波間の位相差を元にした代替 PIA 推定が推奨される (反射強度と異なり、偏波間の位相差は降雨による減衰の影響を受けにくいため)。AhA_{h} と ADRA_{\rm DR} が KDPK_{\rm DP} に対し線形に比例すると仮定した場合、α\alpha と β\beta は次のように定義出来る。

    α=AhKDP(1)β=ADPKDP(2)\quad \alpha = \cfrac{A_{h}}{K_{\rm DP}} \quad (1) \\\quad \beta = \cfrac{A_{\rm DP}}{K_{\rm DP}} \quad (2) \\

    これらを range 方向に沿って一定と仮定する場合、ZZ と ZDRZ_{\rm DR} の減衰量は偏波間位相差の総距離に比例すると考え、それぞれ下記のように定義出来る。

    ΔZh(r)=2∫0rAh(s)ds=2α∫0rKDP(s)ds=αΦDP(r)(3)ΔZDR(r)=2∫0rADP(s)ds=2β∫0rKDP(s)ds=βΦDP(r)(4)\quad \Delta Z_{h}(r) = 2 \displaystyle\int_{0}^{r} A_{h}(s){\rm d}s = 2 \alpha \displaystyle\int_{0}^{r} K_{\rm DP}(s){\rm d}s = \alpha\Phi_{\rm DP}(r) \quad\quad (3) \\\quad \Delta Z_{\rm DR}(r) = 2 \displaystyle\int_{0}^{r} A_{\rm DP}(s){\rm d}s = 2 \beta \displaystyle\int_{0}^{r} K_{\rm DP}(s){\rm d}s = \beta\Phi_{\rm DP}(r) \quad (4) \\

    この手法は Balakrishnan and Zrnic (1989) と Bringi et al. (1990) により提案された。Testud et al. (2000) では、式の導出過程がやや異なるものの同様に減衰量が求められる ZPHI method が提案されている。後に Le Bouar et al. (2001) で詳細に調べられ、拡張された。

    降水での α\alpha と β\beta は、ZDRZ_{\rm DR} と気温の関数である。特に C バンドレーダーでは、α\alpha と β\beta の値が大きく変動することに課題があった (Carey et al. 2000 など)。このため、最適な α\alpha と β\beta を求めることは、降雨減衰量の推定において非常に重要な役割を果たす。

    この課題に対し、大きく分けて 2 つの異なるアプローチが提案されている。1 つ目は、self-consistent method with constraints (Bringi et al. 2001) である。2 つ目は、“hot-spot” algorithm (Ryzhkov et al. (2007) で提案、Gu et al. (2011) で拡張) である。以下、これら 2 つの手法についてそれぞれ述べる。

    ZPHI method の2つの系統

    A self-consistent method with constraints

    ΦDP\Phi_{\rm DP} の range 方向のプロファイルを考えた時に、減衰する降水セルの遠方での ZDRZ_{\rm DR} 補正値に対し制約条件を課すことで、α\alpha と β\beta の値を最適化する手法 (Bringi et al. 2001) である。

    C バンドレーダーの場合、式 (3) と 式 (4) における α\alpha と β\beta の代表的な値は、それぞれ 0.08 dB deg−1^{\rm -1} と 0.02 dB deg−1^{\rm -1} である (see Ryzhkov and Zrnic 2019, Table 6.4)。しかしながら、その変動幅は大きく、range 毎に最適な α\alpha と β\beta の値は異なるであろう。

    そこで、観測される ΦDP\Phi_{\rm DP} と事前計算する ΦDP(mod)\Phi_{\rm DP}^{\rm (mod)} との差 Δ\Delta が、レンジ方向において最小となる α(opt)\alpha^{\rm (opt)} を求める (Bringi et al. 2001)。事前計算する ΦDP(mod)\Phi_{\rm DP}^{\rm (mod)} と差 Δ\Delta は、それぞれ次式で与えられる。

    ΦDP(mod)(r,α)=2∫0rAh(s,α)αds(5)Δ=∑i=1N∣ΦDP(mod)(ri,α)−ΦDP(ri)∣(6)\quad \Phi_{\rm DP}^{\rm (mod)}(r,\alpha) = 2 \displaystyle\int_{0}^{r} \displaystyle\cfrac{A_{h}(s, \alpha)}{\alpha}{\rm d}s \quad (5) \\\quad \Delta = \displaystyle\sum_{i=1}^{N} | \Phi_{\rm DP}^{\rm (mod)}(r_{\rm i},\alpha) - \Phi_{\rm DP}(r_{\rm i}) | \quad (6) \\

    ここで i\rm i は、range 方向の gate 位置を示す index である。Le Bouar et al. (2001) は、α(opt)\alpha^{\rm (opt)} を求める際に、range 方向の隣り合う偏波間位相差の差が 6 ∘^{\circ} 以上にすることを提案している。これにより、結果として得られる降水強度のノイズを減らせることが報告されている。

    “Hot-spot” algorithm

    α\alpha β\beta を non hot spot と hot spot とで分けることで、最適値を求める手法 (Gu et al. 2011) である。

    まず、non hot spot で α0\alpha_{\rm 0} β0\beta_{\rm 0}、hot spot で Δα\Delta\alpha Δβ\Delta\beta とし、次式を定義する。

    α=α0+Δα(7)β=β0+Δβ(8)\quad \alpha = \alpha_{\rm 0} + \Delta\alpha \quad (7) \\\quad \beta = \beta_{\rm 0} + \Delta\beta \quad (8) \\

    次に、α0\alpha_{\rm 0} β0\beta_{\rm 0} を次式で求める (Gu et al. 2011)。

    ΔZh(r)=α0ΦDP(r)(9)ΔZDR(r)=β0ΦDP(r)(10)\quad \Delta Z_{h}(r) = \alpha_{\rm 0} \Phi_{\rm DP}(r) \quad (9) \\\quad \Delta Z_{\rm DR}(r) = \beta_{\rm 0} \Phi_{\rm DP}(r) \quad (10) \\

    そして、ZZ と ZDRZ_{\rm DR} の減衰量は次式で表される。

    ΔZh(r)=α0ΦDP(r)+ΔαΔΦDP(HS)(11)ΔZDR(r)={β0ΦDP(r)if r<r1(12)β0ΦDP(r)+Δβ[ΦDP(r)−ΦDP(r1)]if r1<r<r2(13)β0ΦDP(r)+ΔβΔΦDP(HS)if r>r2(14)\quad \Delta Z_{h}(r) = \alpha_{\rm 0} \Phi_{\rm DP}(r) + \Delta\alpha \Delta\Phi_{\rm DP}(HS) \quad (11) \\\quad \Delta Z_{\rm DR}(r) = \begin{cases} \beta_{\rm 0} \Phi_{\rm DP}(r) &\text{if } r < r_{\rm 1} \quad \quad \quad (12) \\ \beta_{\rm 0} \Phi_{\rm DP}(r) + \Delta\beta [\Phi_{\rm DP}(r) - \Phi_{\rm DP}(r_{\rm 1})] &\text{if } r_{\rm 1} < r < r_{\rm 2} \quad (13) \\ \beta_{\rm 0} \Phi_{\rm DP}(r) + \Delta\beta \Delta\Phi_{\rm DP}(HS) &\text{if } r > r_{\rm 2} \quad \quad \quad (14) \\\end{cases}

    ここで、r1r_{\rm 1} 及び r2r_{\rm 2} は、range 方向の r1r_{\rm 1} から r2r_{\rm 2} の範囲に強雨が観測されたと仮定した場合の gate 位置である。この時、次式を満たす最適な Δα\Delta\alpha Δβ\Delta\beta をそれぞれ求める。

    ∫OHSAh(s,Δα)ds=α02ΔΦDP(OHS)(15)ZDR(a)(rm)+β0ΔΦDP(r0,rm)+ΔβΔΦDP(HS)=ZDR′(rm)(16)\quad \displaystyle\int_{\rm OHS} \displaystyle A_{h}(s, \Delta\alpha){\rm d}s = \displaystyle\cfrac{\alpha_{\rm 0}}{2} \Delta\Phi_{\rm DP}(OHS) \quad (15) \\\quad Z_{\rm DR}^{(a)}(r_{\rm m}) + \beta_{\rm 0} \Delta\Phi_{\rm DP}(r_{\rm 0}, r_{\rm m}) + \Delta\beta \Delta\Phi_{\rm DP}(HS) = Z_{\rm DR}^{'}(r_{\rm m}) \quad (16) \\

    ここで、ΔΦDP(OHS)=ΔΦDP(r0,rm)+ΔΦDP(HS)\Delta\Phi_{\rm DP}(OHS) = \Delta\Phi_{\rm DP}(r_{\rm 0}, r_{\rm m}) + \Delta\Phi_{\rm DP}(HS) で与えられる。また、OHS が non hot spot、HS が hot spot の意。r0r_{\rm 0} はレーダーサイトの gate 位置で、rmr_{m} は range 方向の r2r_{\rm 2} よりも遠方の適当な gate 位置である。

    この時、次の 1. から 4. までの条件を全て満たすものを hot spot として扱う (Gu et al. 2011)。言い換えると、この条件を満たさない gate は全て non hot spot として扱う、ということになる。

    1. ZZ > Z(th)Z^{\rm (th)} かつ ρhv\rho_{\rm hv} > 0.7
    2. hot spot 内での ZDRZ_{\rm DR} の最大値 > ZDR(th)Z_{\rm DR}^{\rm (th)} (通常 3 dB を超える)
    3. hot spot の range 方向の長さが十分に長い (通常 2 km 以上)
    4. hot spot 内の ΦDP\Phi_{\rm DP} が ΔΦDP(th)\Delta\Phi_{\rm DP}^{\rm (th)} よりも大きい (典型的には ΔΦDP(th)\Delta\Phi_{\rm DP}^{\rm (th)} = 10∘^{\circ})。

    以上により、hot spot 内外で最適な係数のペアが求められる。なお、Gu et al. (2011) の “Hot-spot” algorithm は、Py-ART の calculate_attenuation として実装されている。実例はこちら。

    おわりに

    今回は、降雨減衰補正 (主に ZPHI method) について変遷と系統を把握した。

    参考文献

    更新履歴

    ]]>     @@ -227,7 +227,7 @@ https://takashiunuma.github.io/blog/2022/01/26/ 2022-01-26T14:11:06.000Z - 2024-01-20T10:05:02.447Z + 2024-04-15T14:13:04.503Z はじめに

    気象レーダーによって観測された要素は、そのまま利用するのではなく校正を要する。多くの場合、レーダー設置時や点検の際に校正が行われる。しかしながら、校正が不十分な場合や経年劣化によって、再校正が必要になる場合がある。校正がサイト側で行われるのならば、データの利用者は特に気にしなくて良いだろう。そうではないならば、データの利用者が適切に校正を行うとともにバイアス補正をする必要がある。
    といっても真値が分からないので、多くの場合散乱シミュレーションで得た結果との差をバイアスとするようだ (Goddard et al. 1994; 増田 2016)。今回は、主に Goddard et al. (1994) の手法により、データ利用者で行う校正とそれからバイアス補正値の算出について述べる。

    方針

    反射強度のバイアス補正値算出

    まず、散乱シミュレーションで反射強度・反射因子差・偏波間位相差変化率を算出する。Self consistency principleと同様に、Pytmatrix を利用した。帯域はC、温度 10°C^{\degree}\rm C、ガンマ分布の修正パラメータ 5.0、雨滴の縦横比関数 Thurai et al. (2007) とした。
    その後、実データの偏波間位相差変化率を用いて反射強度を算出。式は下記。

    ZH(sim)=10C2log⁡10(KDP(obs)C1)\quad Z_{\rm H(sim)} = \cfrac{10}{C_{\rm 2}} \log_{\rm 10} \left( \cfrac{\rm K_{DP(obs)}}{C_{\rm 1}} \right)

    ここで、添字の obs はレーダー観測値、sim は散乱シミュレーションによる値。NNは観測データ数。C1,C2C_{\rm 1}, C_{\rm 2} は散乱シミュレーションにより求めた定数。C1=0.0003,C2=0.8C_{\rm 1} = 0.0003, C_{\rm 2} = 0.8。

    バイアスが存在しないと仮定した偏波間位相差変化率を用い、散乱シミュレーションから得た平均的な反射強度と偏波間位相差変化率との自己整合性から反射強度のバイアス補正値を算出する。式は下記。

    BIASZH=−1N∑i=1N(ZH(obs),i−ZH(sim),i)\quad BIAS_{\rm Z_{\rm H}} = - \cfrac{1}{N} \displaystyle\sum_{\rm i=1}^{N} \left( Z_{\rm H(obs), i} - Z_{\rm H(sim), i} \right)

    反射因子差のバイアス補正値算出

    反射因子差のバイアス補正値については、後に述べるバイアス補正済の反射強度を用い、算出。式は下記。

    ZDR(sim)=D1Z∗D2\quad Z_{\rm DR(sim)} = D_{\rm 1} Z_{\ast}^{D_{\rm 2}}

    ここで、Z∗=10ZH(obs)−BIASZHZ_{\ast} = 10^{Z_{\rm H(obs)} - BIAS_{Z_{\rm H}}}。添字の obs はレーダー観測値、sim は散乱シミュレーションによる値。NNは観測データ数。D1,D2D_{\rm 1}, D_{\rm 2} は散乱シミュレーションにより求めた定数。D1=0.11,D2=0.25D_{\rm 1} = 0.11, D_{\rm 2} = 0.25。

    バイアス補正済みの反射強度を用い、反射因子差と反射強度との自己整合性から反射因子差のバイアスを算出する。式は下記。

    BIASZDR=−1N∑i=1N(ZDR(obs),i−ZDR(sim),i)\quad BIAS_{\rm Z_{\rm DR}} = - \cfrac{1}{N} \displaystyle\sum_{\rm i=1}^{N} \left( Z_{\rm DR(obs), i} - Z_{\rm DR(sim), i} \right)

    おわりに

    利用者側で行う、反射強度・反射因子差のバイアス補正について述べた。

    参考文献

    ]]>     @@ -250,7 +250,7 @@ https://takashiunuma.github.io/blog/2021/10/01/ 2021-10-01T12:42:38.000Z - 2024-01-20T10:18:09.537Z + 2024-04-15T14:13:04.503Z はじめにSelf consistency principle の記事の中で、雨滴粒径分布のシミュレーション・その結果を用いたレーダー変数推定を扱った。文献や教科書を追っていると、そもそも気象レーダーがどんなものを対象とし、どのような仮定の下で対象を観測しているのかについて、適切に理解できていなかった。そこで、気象レーダーが粒径分布からレーダー変数を計算するために必要な事項を Bringi and Chandrasekar (2001) と深尾・浜津 (2005) を基にまとめることにした。

    散乱領域

    レイリー領域・ミー領域・光学領域がある。観測可能な標的のサイズは、レーダーの周波数帯によって異なる。

    リモートセンシングの対象標的

    ハードターゲット: 孤立型標的、離散分布型標的
    ソフトターゲット: 連続分布型標的

    気象レーダーが対象とする雨・雲粒のような微小な離散分布型標的は、散乱体積内の微小散乱断面積の総和として得られる。雨粒などの微小な降水粒子は、レイリー領域に位置づけられ、レイリー近似が成り立つ (例えば、5 GHz 帯の場合、標的のサイズが 3 cm 程度まで)。

    粒径分布

    降水強度をレーダー断面積の関数として捉える場合、粒径分布が密接に関係する。一般的な粒径分布を表現するものとして、修正ガンマ分布がある。

    N(D)=N0Dμexp⁡(−ΛD)\quad N(D) = N_{\rm 0} D^{\mu} \exp(-\Lambda D) \\
    N0N_{0} (分布曲線の切片, intercept parameter)・μ\mu (ガンマ分布の修正パラメータ)・Λ\Lambda (分布曲線の傾斜, slope parameter)。

    後方散乱行列 (Backscattering matrix)

    単位面積当たりの散乱断面積と後方散乱行列の共分散行列の各要素とを関連付けることができる。後方散乱受信信号電圧を散乱体積全体で総和をとり、それとその複素共役の積から2次モーメントを得る。この2次モーメントの一部が可逆なので、共分散行列自体は 3x3 になる。これらの要素から偏波パラメータを得る。偏波パラメータには、後方散乱によるものと伝搬によるものとがある。

    まず、散乱波を水平偏波と垂直偏波の合成によるものと仮定すると、アンテナ点での後方散乱電界 Eb\bm{E}^{\rm b} は入射電界 Ei\bm{E}^{\rm i} と直線偏波の後方散乱行列 s\bm{s} により、次式のように表される (後方散乱はアンテナでの受信、入射はアンテナからの送信と言い換えることができる)。

    [EhEv]b=s[EhEv]ie−jkrr\qquad\begin{bmatrix} E_{\rm h} \\ E_{\rm v}\end{bmatrix}^{\rm b}=\bm{s}\begin{bmatrix} E_{\rm h} \\ E_{\rm v}\end{bmatrix}^{\rm i}\cfrac{e^{-jkr}}{r}

    ここで、

    s=[shhshvsvhsvv]\qquad\bm{s} =\begin{bmatrix} s_{\rm hh} & s_{\rm hv} \\ s_{\rm vh} & s_{\rm vv}\end{bmatrix}

    とし、kk は波数、rr はレーダーと標的の間の距離。電界ベクトルや行列要素の添字のうち、h\rm h は水平偏波を、v\rm v は垂直偏波を示す。第1添字は散乱電界の偏波を、第2添字は入射電界の偏波を示す。

    散乱波を右旋円偏波と左旋円偏波の合成と見做せるときは、s\bm{s} を円偏波の後方散乱行列 sc\bm{s}_{\rm c} に置き換えることができる。

    sc=[srrsrlslrsll]\qquad\bm{s}_{\rm c}=\begin{bmatrix} s_{\rm rr} & s_{\rm rl} \\ s_{\rm lr} & s_{\rm ll} \\\end{bmatrix}
    s\bm{s} と sc\bm{s}_{\rm c} には下記のような関係がある。
    srr=12(svv−shh−j2svh)sll=12(svv−shh+j2svh)srl=slr=12(svv+shh)\qquad\begin{aligned}s_{\rm rr} &= \cfrac{1}{2} (s_{\rm vv} - s_{\rm hh} - j2s_{\rm vh} ) \\s_{\rm ll} &= \cfrac{1}{2} (s_{\rm vv} - s_{\rm hh} + j2s_{\rm vh} ) \\s_{\rm rl} &= s_{\rm lr} = \cfrac{1}{2} (s_{\rm vv} + s_{\rm hh} )\end{aligned}

    次に、単位入射電界に対する個々の降水粒子による散乱を考えるとき、レーダーから距離 rnr_n にある nn 番目の粒子による後方散乱受信信号電圧 vijv_{ij} は、下記のように表される。

    vij(rn)=sij(n)F(rn)e−j2πkrn\qquadv_{ij}(\bm{r}_{n}) = s_{ij}(n)F(\bm{r}_{n})e^{-j 2 \pi k r_{n}}

    散乱体積全体の受信信号 VijV_{ij} は上記の総和となり、下記で与えられる。

    Vij(r)=∑nsij(n)F(rn)e−j2πkrn\qquadV_{ij}(r) = \displaystyle\sum_{n} s_{ij}(n) F(\bm{r}_{n}) e^{-j 2 \pi k r_{n}}
    VijV_{ij} とその複素共役の積で与えられる2次モーメントは、⟨VijVkl∗⟩=⟨∑n∑msij(n)skl∗(m)e−j2πk(rn−rm)F(rn)F∗(rm)⟩=∑n⟨sij(n)skl∗(n)⟩∣F(r)∣2=∫⟨n(r)sijskl∗⟩∣F(r)∣2dV\qquad\begin{aligned}\lang V_{ij} V_{kl}^{\ast} \rang &= \left\lang \displaystyle\sum_{n} \displaystyle\sum_{m} {s_{ij}(n) s^{\ast}_{kl}(m) e^{-j 2 \pi k (r_{n} - r_{m})} F(\bm{r}_{n}) F^{\ast}(\bm{r}_{m})} \right\rang \\ &= \displaystyle\sum_{n} \lang {s_{ij}(n) s^{\ast}_{kl}(n) \rang |F(\bm{r})|^{\rm 2}} \\ &= \displaystyle\int \lang n(r) s_{ij} s^{\ast}_{kl} \rang |F(\bm{r})|^{\rm 2} {\rm d}V\end{aligned}

    となる。ここで ⟨∗⟩\lang \ast \rang は期待値、n(r)n(r) は位置 rr での単位体積あたりの粒子の大きさの分布。上式で示される2次モーメントは、一般に 4x4 の共分散行列で表される。ここで、VijV_{ij} と VjiV_{ji} とが可逆となる性質を利用すると共分散行列は 3x3 となる。したがって、各要素は

    [⟨n∣shh∣2⟩⟨nshvshh∗⟩⟨nsvvshh∗⟩⟨nshhshv∗⟩⟨n∣shv∣2⟩⟨nsvvshv∗⟩⟨nshhsvv∗⟩⟨nshvsvv∗⟩⟨n∣svv∣2⟩]\qquad\begin{bmatrix} \lang n |s_{\rm hh}|^{\rm 2} \rang & \lang n s_{\rm hv} s^{\ast}_{\rm hh} \rang & \lang n s_{\rm vv} s^{\ast}_{\rm hh} \rang \\ \lang n s_{\rm hh} s^{\ast}_{\rm hv} \rang & \lang n |s_{\rm hv}|^{\rm 2} \rang & \lang n s_{\rm vv} s^{\ast}_{\rm hv} \rang \\ \lang n s_{\rm hh} s^{\ast}_{\rm vv} \rang & \lang n s_{\rm hv} s^{\ast}_{\rm vv} \rang & \lang n |s_{\rm vv}|^{\rm 2} \rang\end{bmatrix}

    となる。

    粒子の形状と偏波パラメータ

    上記までで粒径分布と共分散行列がそれぞれ分かった。以前、粒径と扁平度の関係を Self consistency principle の雨滴粒径分布の節で述べ、その結果からレーダー変数を計算した。得られた共分散行列の各要素を粒子形状と関連付けることで、それらの関係を明らかにする。
    扁平な雨滴の長軸及び短軸の径をそれぞれ aa および bb としたとき、軸比を b/ab/a、扁平率に依存する因子 λz\lambda_{z} を

    λz=1+f2f2(1−1ftan⁡−1f)\qquad\lambda_{z} = \cfrac{1 + f^{\rm 2}}{f^{\rm 2}} \left(1 - \cfrac{1}{f} \tan^{\rm -1} f \right)

    とする。ただし、

    f2=(ab)2−1\qquadf^{\rm 2} = \left( \cfrac{a}{b} \right)^{\rm 2} - 1

    とすると、shhs_{\rm hh} および svvs_{\rm vv} は Rayleigh-Gunn の理論により、それぞれ下記のように表される。

    shh=k24πV(εr−1)1+12(1−λz)(εr−1)svv=k24πV(εr−1)1+λz(εr−1)\qquad\begin{aligned}s_{\rm hh} &= \cfrac{k^{\rm 2}}{4 \pi} \cfrac{V (\varepsilon_{\rm r} - 1)}{1 + \frac{1}{2}(1 - \lambda_{z})(\varepsilon_{\rm r} - 1)} \\s_{\rm vv} &= \cfrac{k^{\rm 2}}{4 \pi} \cfrac{V (\varepsilon_{\rm r} - 1)}{1 + \lambda_{z}(\varepsilon_{\rm r} - 1)}\end{aligned}

    ここで、kk はレーダー電波の波数で k=2π/λk = 2 \pi/\lambda、VV は粒子体積で V=(π/6)D3V = (\pi/6) D^{3}、ε\varepsilon は比誘電率 (=ε/ε0)(= \varepsilon/\varepsilon_{0})。

    粒径分布が N(D)N(D) になる場合を考える。shhs_{\rm hh} と svvs_{\rm vv} は共に b/ab/aとDDの関数であることを加味すると、レーダー反射因子は下記のように書き換えられる。

    Zhh=4λ4π4∣Kw∣2⟨n∣shh∣2⟩=∫D6N(D){1+λz(ε−1)}2dDZvv=4λ4π4∣Kw∣2⟨n∣svv∣2⟩=∫D6N(D){1+12 (1−λz)(ε−1)}2dD\qquad\begin{aligned}Z_{\rm hh} &= \cfrac{4 \lambda^{\rm 4}}{\pi^{\rm 4} |K_{\rm w}|^{\rm 2}} \lang n |s_{\rm hh}|^{\rm 2} \rang \\ &= \displaystyle\int D^{6} N(D) \{1 + \lambda_{z} (\varepsilon - 1) \}^{2} {\rm d}D \\Z_{\rm vv} &= \cfrac{4 \lambda^{\rm 4}}{\pi^{\rm 4} |K_{\rm w}|^{\rm 2}} \lang n |s_{\rm vv}|^{\rm 2} \rang \\ &= \displaystyle\int D^{6} N(D) \left\{1 + \cfrac{1}{2} \, (1 - \lambda_{z})(\varepsilon - 1) \right\}^{2} {\rm d}D \\\end{aligned}

    反射因子差は、真数表現を ZdrZ_{\rm dr} とすると、レーダー反射因子差を用い下記のように表される。

    Zdr=4π⟨n∣shh∣2⟩4π⟨n∣svv∣2⟩=∫D6N(D){1+λz(ε−1)}2dD∫D6N(D){1+12 (1−λz)(ε−1)}2dD≃∫D6N(D)dD∫(ba)7/3D6N(D)dD\qquad\begin{aligned}Z_{\rm dr} &= \cfrac{4 \pi \lang n |s_{\rm hh}|^{\rm 2} \rang}{4 \pi \lang n |s_{\rm vv}|^{\rm 2} \rang} \\ &= \cfrac{\displaystyle\int D^{6} N(D) \{1 + \lambda_{z} (\varepsilon - 1) \}^{2} {\rm d}D}{\displaystyle\int D^{6} N(D) \left\{1 + \cfrac{1}{2} \, (1 - \lambda_{z})(\varepsilon - 1) \right\}^{2} {\rm d}D} \\ &\simeq \cfrac{\displaystyle\int D^{6} N(D) {\rm d}D}{\displaystyle\int \left(\cfrac{b}{a} \right)^{7/3} D^{6} N(D) {\rm d}D}\end{aligned}

    ここで、

    ∣1+λz(εr−1)∣2∣1+12 (1−λz)(εr−1)∣2≃1(ba)7/3\qquad\cfrac{|1 + \lambda_{z}(\varepsilon_{r} - 1)|^{2}}{\left| 1 + \cfrac{1}{2} \, (1 - \lambda_{z})(\varepsilon_r - 1) \right|^{2}} \simeq \cfrac{1}{\left(\cfrac{b}{a} \right)^{7/3}}

    の近似を用いた (Bringi and Chandrasekar 2001, eq. 7.19)。

    偏波間相関係数は、下記のように表される (Bringi and Chandrasekar 2001, eq. 7.41, 7.42, 7.43, 7.44, 7.55)。

    ρhv(0)=⟨nsvvshh∗⟩⟨n∣shh∣2⟩1/2⟨n∣svv∣2⟩1/2=ZdrZh ∫N(D)∫∣Shh∣2r7/6p(r)drdD≃rz7/6(rz7/6)1/2≃1−12 var(rz)(rˉz)2≃1−12 Λ2[f(μ)]−2 var(D){1−Λ[D/f(μ)]}2≃1−12 Λ2[f(μ)]−2(D2){1−Λ[D/f(μ)]}2(μ+3)\qquad\begin{aligned}\rho_{\rm hv} ({\rm 0}) &= \cfrac{\lang n s_{\rm vv} s^{\ast}_{\rm hh} \rang}{\lang n |s_{\rm hh}|^{2} \rang^{1/2} \lang n |s_{\rm vv}|^{2} \rang^{1/2}} \\ &= \cfrac{\sqrt{Z_{dr}}}{Z_{h}} \, \displaystyle\int N(D) \displaystyle\int |S_{\rm hh}|^{2} r^{7/6} p(r) dr dD \\ &\simeq \cfrac{r^{7/6}_{z}}{(r^{7/6}_{z})^{1/2}} \\ &\simeq 1 - \cfrac{1}{2} \, \cfrac{ {\rm var}(r_{z})}{(\bar{r}_z)^{2}} \\ &\simeq 1 - \cfrac{1}{2} \, \cfrac{\Lambda^{2}[f(\mu)]^{-2} \, {\rm var}(D)}{\{1 - \Lambda [D/f(\mu)]\}^{2}} \\ &\simeq 1 - \cfrac{1}{2} \, \cfrac{\Lambda^{2}[f(\mu)]^{-2}(D^{2})}{\{1 - \Lambda [D/f(\mu)]\}^{2} (\mu + 3)} \\\end{aligned}

    式変形には、

    rz=1−ΛDzvar(rz)=γ2var(Dz)D=Γ(μ+4)Γ(μ+3) Γ(μ+6)Γ(μ+7)Dz=f(μ)Dzvar(D)=D2(μ+3)\qquad\begin{aligned}r_{z} &= 1 - \Lambda D_{z} \\{\rm var}(r_{z}) &= \gamma^{2} {\rm var}(D_{z}) \\D &= \cfrac{\Gamma(\mu + 4)}{\Gamma(\mu + 3)} \, \cfrac{\Gamma(\mu + 6)}{\Gamma(\mu + 7)} D_{z} \\ &= f(\mu) D_{z} \\{\rm var}(D) &= \cfrac{D^{2}}{(\mu + 3)}\end{aligned}

    の関係を用いた。なお、完全ガンマ関数は、

    ∫0∞xa−1 e−bxdx=1ba Γ(a)\qquad\begin{aligned}\displaystyle\int^{\infty}_{0} x^{a-1} \, e^{-bx} {\rm d}x = \cfrac{1}{b^{a}} \, \Gamma(a)\end{aligned}

    である。ただし、a>1,b>1a > 1, b > 1 とする。

    偏波間位相差変化率は、下記のように表される。これまでのレーダー変数が後方散乱によるものであったのに対し、偏波間位相差変化率は伝搬によるものである。

    KDP=2πk∫N(D) Re [h⋅f(b/a,D)−v⋅f(b/a,D)] dD≃πk12∫D3 N(D) Re [F(b/a,εr)] dD=πkCk12∫D3(1−b/a) N(D)dD=(πλ)Ckρw∫π6ρwD3(1−b/a) N(D)dD\qquad\begin{aligned}K_{\rm DP} &= \cfrac{2\pi}{k} \displaystyle\int N(D) \, {\rm Re} \, [\bm{h} \cdot \bm{f}(b/a, D) - \bm{v} \cdot \bm{f}(b/a, D)] \, {\rm d}D \\ &\simeq \cfrac{\pi k}{12} \displaystyle\int D^{3} \, N(D) \, {\rm Re} \, [F(b/a, \varepsilon_{r})] \, {\rm d}D \\ &= \cfrac{\pi k C_{k}}{12} \displaystyle\int D^{3} (1 - b/a) \, N(D) {\rm d}D \\ &= \left(\cfrac{\pi}{\lambda} \right) \cfrac{C_{k}}{\rho_{\rm w}} \displaystyle\int \cfrac{\pi}{6} \rho_{\rm w} D^{3} (1 - b/a) \, N(D) {\rm d}D\end{aligned}

    ここで、前方へ散乱される電波のベクトル振幅を f(b/a,D)\bm{f}(b/a, D)、h,v\bm{h}, \bm{v} はそれぞれ水平偏波・垂直偏波の単位ベクトル、ρw\rho_{\rm w} は降水粒子の密度 [g m−1^{-1}]、CkC_{k} は (周波数依存だけどもほぼ) 定数。式変形には、Rayleigh-Gunn の理論により、

    Re [h⋅f(b/a,D)−v⋅f(b/a,D)]=Re [F(b/a,εr)]=Re [εr−11+12 (1−λz)(εr−1)−εr−11+λz(εr−1)]≃Ck(1−b/a)\qquad\begin{aligned}{\rm Re} \, [\bm{h} \cdot \bm{f}(b/a, D) - \bm{v} \cdot \bm{f}(b/a, D)] &= {\rm Re} \, [F(b/a, \varepsilon_{r})] \\ &= {\rm Re} \, \left[\cfrac{\varepsilon_{r} - 1}{1 + \cfrac{1}{2} \, (1 - \lambda_{z})(\varepsilon_r - 1)} - \cfrac{\varepsilon_{r} - 1}{1 + \lambda_{z}(\varepsilon_{r} - 1)} \right] \\ &\simeq C_{k}(1 - b/a)\end{aligned}

    の関係を用いた。

    偏波間位相差変化率の式の形をよく見ると、∫π6 ρwD3N(D)dD\displaystyle\int \cfrac{\pi}{6} \, \rho_{\rm w} D^{3} N(D) {\rm d}D の部分は単位体積中の含水量となっている。二重偏波レーダーとメソ対流系で考察したように、降水強度の推定に主として偏波間位相差変化率が利用されるようになった一因と解釈できうる。

    おわりに

    Bringi and Chandrasekar (2001)、深尾・浜津 (2005) を基に、粒径分布から二重偏波気象レーダーの各種レーダー変数を算出する過程について扱った。

    参考文献

    更新履歴

    ]]>     @@ -273,7 +273,7 @@ https://takashiunuma.github.io/blog/2021/09/20/ 2021-09-20T13:13:20.000Z - 2024-01-20T10:14:01.207Z + 2024-04-15T14:13:04.502Z はじめにSignal Processing Methodsや二重偏波要素の特徴では、主にレーダー変数の意味やその分布の具体例を示した。ところが、得られるレーダー変数による現象の解釈については十分に出来ていなかった。そこで、今回はメソ対流系を二重偏波レーダーで観た場合、どんなことが分かるのか調べてみた。

    メソ対流系

    メソ対流系の典型例として米国でしばしば観測される squall line をドップラーレーダーで観測され、その内部構造が鉛直構造とともに調べられた (Houze et al. 1989)。
    主に squall line の進行方向全面に上昇流を伴う対流域と、その後面に層状域とで構成される。気流構造としては、進行方向全面に冷気外出流による gust front とその冷気外出流を滑昇する流れ、後面には対流域で上昇した流れが移動方向後面に向かって緩やかに上昇しつつ層状域を形成する流れと補償流として後面から全面に向かって流れ込む流れとがある。
    これらの構造は、二重偏波レーダー変数ではどのように見えるのだろうか。

    二重偏波レーダー変数とメソ対流系

    Ryzhkov and Zrnic (2019) Sec. 8.1 を一部意訳し、特徴を列挙する。

    1つ目の反射強度が高いというのは、単偏波の場合でも同様の特徴である。反射因子差の高い部分というのが、二重偏波化によって新たに得られる情報となる。
    2つ目の対流の上昇流域で反射因子差が大きくなるというのは、鉛直構造を観ると ZDRZ_{DR} column として特徴づけられる (e.g., Kumjian et al, 2014; Snyder et al. 2015)。雨粒の粒径が異なるというのは、いわゆる上昇流域内で雨粒の衝突併合が生じ、大きさや特徴の異なる雨粒が多く存在しうる状態になることを意味する。二重偏波要素の特徴で述べたように ZDRZ_{DR} は雨粒の形状を示す特徴量であるから、衝突併合が起こると雨粒で言うと大きく扁平なものが多くなり結果として ZDRZ_{DR} が大きな値を示す。上昇流域であるからそれが鉛直方向に分布し、ZDRZ_{DR} column として観測される、という解釈である。
    3つ目・4つ目は、Signal Processing Methodsや二重偏波要素の特徴で述べたことと概ね同様であるため省略する。

    ZDRZ_{DR} column

    先に述べた ZDRZ_{DR} column は、対流雲内の上昇流と対応が良いこと (Snyder et al. 2015) ZDRZ_{DR} column の深さが上昇流の強さと比例関係にあること (Kumjian et al, 2014) がそれぞれ知られている。おしなべて、発達前の対流雲を ZDRZ_{DR} column で診断出来るのではないか、という考え方である。Snyder et al. (2015) では、ZDRZ_{DR} column の深さと高度 4 km の上昇流 (最大値) との時刻毎の相関関係から、ZDRZ_{DR} column が深くなってから 2 分後に高度 4 km の上昇流ピーク値が生じることを示した。Kumjian et al, (2014) は、ZDRZ_{DR} column に関するレビューをしつつ、ZDRZ_{DR} column の life cycle や ZDRZ_{DR} column の深さ or 体積の上昇が地上での降水 or 降雹の 10〜15 分前に生じうることをそれぞれ示した。対流雲が発達してから降水を伴うまでに時間差があるので、ZDRZ_{DR} column の発現からだいたい数分〜数10分程度ということだろうか。雹を伴う系だと ZDRZ_{DR} column のシグナルも大きく有効と考えられる。降水または降雹とあるので、単純に降水だけだと ZDRZ_{DR} column のシグナルが小さいかもしれない。

    ZDRZ_{DR} column で上昇流が診断出来うるのならば、もしかして下降流の診断も出来るのでは？と考えていたら KDPK_{DP} column というものがあった。

    KDPK_{DP} column

    KDPK_{DP} column は、基本的には対流雲内の降水域 (下降流域) と対応することが知られている (Loney et al. 2002)。KDPK_{DP} column の起源は ZDRZ_{DR} column と似た場所となることが多いものの、時間と共に別々の場所で観測される (ZDRZ_{DR} column は上昇流域・KDPK_{DP} column は下降流域に対応する Ryzhkov and Zrnic, 2019: Sec. 8.2)。Loney et al. (2002) では、ZDRZ_{DR} column の上端 (融解層上部付近) で KDPK_{DP} column が観測されており、内部は一部凍った雨粒が対応しているそうだ。KDPK_{DP} column が観測されるタイミングなどは特に触れられていなかった。ただ、察するに雨粒として落下する際に KDPK_{DP} の値が大きいものとして観測され、結果として KDPK_{DP} column になっていると考えられる。降水強度の推定に KDPK_{DP} の利用が主流になったのは、落下する雨粒を捉えているという点が考慮されているからかもしれない。ZDRZ_{DR} だと先に述べたように上昇流に反応してしまうので、地上雨量として換算すると時刻ズレが大きくなり、地上雨量との対応が良くないのだろう。

    おわりに

    今回は二重偏波レーダー変数とメソ対流系について、ZDRZ_{DR}・KDPK_{DP} columns に着目しつつ文章だけで簡単に示した。

    参考文献

    ]]>     @@ -296,7 +296,7 @@ https://takashiunuma.github.io/blog/2021/09/04/ 2021-09-04T13:19:44.000Z - 2024-01-20T10:20:08.277Z + 2024-04-15T14:13:04.494Z はじめに

    Self-consistency principle とは、Goddard et al. (1994) で提唱された、水平偏波反射強度 (Z, 以降、反射強度と略す) の校正方法の一つ。雨に対し、KDP/Z = ZDRK_{\rm DP}/Z \, = \, Z_{\rm DR} の関係から様々な粒径分布をシミュレーションしたものを真値として反射強度の校正をする、というもの。ただし、Goddard et al. (1994) の手法は降雨減衰が考慮されていないため、降雨減衰の影響が無視出来ない X・C バンドレーダーでは利用しにくいという欠点があった。Adachi et al. (2015) ではその欠点を補う、降雨減衰を考慮した校正方法が提案されている。

    二重偏波レーダー変数の算出とその補正では、校正済の反射強度が得られていることを前提としていた。降水強度推定に KDPK_{\rm DP} を用いるのであれば、反射強度の校正は必要ないじゃないか、と考えることも出来る。しかしながら、現業レーダーの場合、世界気象機関 (WMO) でレーダー変数のバイアスをある決められた値以下にしましょうという基準 (努力目標に近い) があり、その変数に反射強度も入っているため校正は必須となる。加えて、受信電力の弱い領域では反射強度を用いて降水強度を推定した方が、降水強度の推定精度が総合的に向上する場合が多い。なので、反射強度の校正は、気象レーダーそのものひいては降水強度の量的精度向上にとって必要かつ重要なプロセスなのである。本稿では Self-consistency principle による反射強度の校正について述べる。

    実施環境

    今回は、こちらの ipython notebook に沿ってシミュレーションを行う。その後、こちらの実データを用い、実際に校正が必要かどうか確認する。

    シミュレーション

    実施環境で環境構築したものに加え、下記を導入する。

    conda install pytmatrix

    https://github.com/jleinonen/pytmatrix/wiki

    雨滴粒径分布

    まず、異なる形状をした雨滴データを用意する。雨滴の縦横比を近似的に求める方法が幾つか提案されていて、pytmatrix に実装されている。

    図は下記の通り。横軸が雨滴の直径、縦軸が雨滴の縦横比。どの方法でも雨滴の直径が大きくなるにつれて縦横比が概ね小さくなる分布となる。

    次に、上記で用意した異なる形状をした雨滴データから、雨滴粒径分布 (Drop Size Distribution) を算出する。雨滴粒径分布の推定は、Bringi and Chandrasekar (2001) の式 7.62 に則り規格化したガンマ分布を仮定する。その際、shape parameter (μ\mu)、雨滴の体積直径 (D0D_{\rm 0})、intercept parameter (NwN_{w}) をそれぞれ与える。上図は D0D_{\rm 0} の頻度分布、下図は NwN_{w} の頻度分布。それぞれ、雨滴粒径分布、1立方メートルあたりの密度と解釈できる。ピーク値でみると、D0D_{\rm 0} 1.7 mm 程度、NwN_{w} 2000 mm−1^{\rm -1} m−3^{\rm -3} 程度なので、層状性の降雨に近い分布である。もう少し値の大きい D0D_{\rm 0} も含まれているので、対流性のものも混じっているようだ。

    レーダー変数推定

    上記で得た雨滴粒径分布から、反射強度・反射因子差・偏波間位相差変化率をそれぞれ求める。確認のため、反射強度と反射因子差の頻度分布を下図に示す。反射強度の値が大きくなるとともに、反射因子差の値が大きくなるような分布をしている。反射強度の大きい雨粒ほど扁平な形状をしているようだ。

    f(KDP/Z) = ZDRf(K_{\rm DP}/Z) \, = \, Z_{\rm DR} 分布

    ようやく Self-consistency の関係を確認する準備が整った。横軸に ZDRZ_{\rm DR}、縦軸に KDP/ZK_{\rm DP}/Z をとり頻度分布を描いたものが下図。0.00007 degree km−1^{\rm -1} mm−6^{\rm -6} m3^{\rm 3} をピーク値とし、それ以降は ZDRZ_{\rm DR} の増加とともにべき乗で減少するような分布をしている。今回は雨滴粒径分布をシミュレーションで算出したけれども、可能であればディスドロメータやビデオ観測等を用いて実際の雨滴粒径分布を得た上で真値を算出すると、より高精度の真値が得られると考えられる。

    実例

    上記までで、シミュレーションによる真値が得られた。簡単に言うと、この真値とのズレが反射強度のバイアスと考えることが出来る (反射因子差・偏波間位相差変化率がそれぞれ正しく校正されていることが前提)。サンプルデータで反射強度の校正の余地があるかどうか実際に調べてみる。下図は、暖色系ほど頻度が高く寒色系ほど頻度が低くなっている。頻度の高い分布を辿ると、概ね先程確認したシミュレーションの分布に近いところにある。ただし、先程のシミュレーションでは気温が 10∘C^{\circ}C のみを対象としている一方で実例では全ての仰角分を描画している。このため、全てのデータがシミュレーションの分布に乗る訳ではない。正確に行うならば、用いるデータの仰角・高度情報から気温を推定し、条件に一致するデータのみを取り出す必要があるだろう。

    おわりに

    今回は、Self-consistency principle と銘打って反射強度の校正について扱った。

    参考文献

    更新履歴

    ]]>     @@ -319,7 +319,7 @@ https://takashiunuma.github.io/blog/2021/08/28/ 2021-08-28T13:47:33.000Z - 2024-01-20T10:16:52.027Z + 2024-04-15T14:13:04.503Z はじめに

    Signal Processing Methods では、気象レーダーの観測結果から得られるレーダー変数推定時にバイアスとなりうる要素やレーダー変数の算出・補正について概観した。その後、二重偏波要素の特徴にて ω\omegaradlib を用いて具体例を示した。降水強度推定については、偏波間位相差変化率 (KDPK_{\rm DP}) の推定アルゴリズムについて概観し (降水強度のKDP推定アルゴリズムについて) 、主として偏波間位相差 (ΦDP\Phi_{\rm DP}) の観測値を活かすようなフィルター・平滑化処理を行うことが重要であることを認識した。サイト毎に得られるレーダー変数の算出やその補正については二重偏波要素の特徴で言及したものの、それらの具体的な算出・補正方法については触れていなかった。そこで、本稿では降水強度算出に関連するレーダー変数の算出方法とレーダー変数推定時にバイアスとなりうる要素に対する補正方法について具体例を示す。それらを以て二重偏波レーダー変数の算出・補正について理解を深める。

    用語

    実施環境

    レーダー変数の算出・補正

    反射強度の降雨減衰補正

    Signal Processing Methods で挙げたレーダー変数推定時にバイアスとなりうる要素のうち、attenuation in atmospheric gases and precipitation (大気減衰・降雨減衰) を取り扱う。大気減衰は主に標準大気を仮定するため、多くの場合一定値が用いられ、補正は単純となる。一方、降雨減衰はその減衰量をどのように求めるかによって考え方や算出手法が異なるようだ。ここでは、Z-PHI 法 (Gu et al. 2011) と呼ばれる降雨減衰補正について述べる。Z-PHI 法では、ざっくり言うと反射強度 ZZ の減衰量が KDPK_{\rm DP} に比例することを仮定し、降雨減衰量を求める。

    Py-ART では、calculate_attenuation という関数を用いる。このページを参考にして、下記のようなスクリプトを用意。データはこちら。

    import matplotlib.pyplot as pltfrom matplotlib.colors import ListedColormap, BoundaryNormimport pyartRADAR_NAME = 'sgpcsaprsurcmacI7.c0.20110520.095101.nc'# read in the dataradar = pyart.io.read_cfradial(RADAR_NAME)# remove existing correctionsradar.fields.pop('specific_attenuation')radar.fields.pop('corrected_reflectivity_horizontal')# perform attenuation correctionspec_at, cor_z = pyart.correct.calculate_attenuation(    radar, 0, refl_field='reflectivity_horizontal',    ncp_field='norm_coherent_power', rhv_field='copol_coeff',    phidp_field='proc_dp_phase_shift')radar.add_field('specific_attenuation', spec_at)radar.add_field('corrected_reflectivity_horizontal', cor_z)# create colormapcmap = ListedColormap(["#0FF", "#07C", "#036", "#00F", "#097", "#0A4", "#0D1", "#4F1", "#7C0", "#FF0", "#F80", "#F8F", "#F4A", "#F00"])bounds = [23, 28, 33, 37, 40, 42, 45, 47, 49, 50, 51, 52, 53]norm = BoundaryNorm(bounds, cmap.N)# create the plotfig = plt.figure(figsize=(5, 10))ax1 = fig.add_subplot(311)ax1.set_aspect('equal', 'box')display = pyart.graph.RadarDisplay(radar)display.plot('reflectivity_horizontal', 0, ax=ax1,              vmin=min(bounds), vmax=max(bounds), cmap=cmap, norm=norm, ticks=bounds,             colorbar_label='', title='Raw Reflectivity', axislabels=('','North South distance from radar (km)'))ax2 = fig.add_subplot(312)ax2.set_aspect('equal', 'box')display.plot('specific_attenuation', 0, ax=ax2, vmin=0, vmax=1.0,             colorbar_label='', title='Specific Attenuation', axislabels=('','North South distance from radar (km)'))ax3 = fig.add_subplot(313)ax3.set_aspect('equal', 'box')display = pyart.graph.RadarDisplay(radar)display.plot('corrected_reflectivity_horizontal', 0, ax=ax3,             vmin=min(bounds), vmax=max(bounds), cmap=cmap, norm=norm, ticks=bounds,             colorbar_label='', title='Corrected Reflectivity')plt.suptitle('Attenuation correction using Py-ART', fontsize=16)plt.show()

    図は、上から降雨減衰補正前の反射強度 (dBZ)、降雨減衰量 (dB)、降雨減衰補正後の反射強度 (dBZ) をそれぞれ示す。
    サンプルデータは C バンドのようで、中央図で値が大きい部分で水平偏波反射強度の減衰が生じている。上図と下図を比べると、上図では下図で減衰している部分が補正された状態で計算されている。加えて、ρhv\rho_{hv}の足切り処理も施されており、一定値以下の反射強度は無効値となっている。このようにKDPK_{\rm DP}から降雨による減衰量を推定することが出来る。では、KDPK_{\rm DP}はどのように算出されているのだろうか。次節でその具体を見てみる。

    KDPK_{\rm DP} 算出

    降水強度のKDP推定アルゴリズムについてで述べたように、KDPK_{\rm DP} の算出は ΦDP\Phi_{\rm DP} とセットで行われることが多い。ただ、ΦDP\Phi_{\rm DP} から KDPK_{\rm DP} を算出する際に用いられる式は同じであり、KDP=12 dΦDPdrK_{\rm DP} = \cfrac{1}{2} \, \cfrac{d\Phi_{\rm DP}}{dr} で与えられる (Rauber and Nesbitt, 2018)。ここで、dr は ray 方向の gate 幅。計算する場合は有限差分として、ある gate と一つ前の gate との ΦDP\Phi_{\rm DP} 差分を gate 幅で割る。計算部分は KDP推定アルゴリズムの比較を参照することにして省略し、今回は計算済みのものを用いる。
    import matplotlib.pyplot as pltfrom matplotlib.colors import ListedColormap, BoundaryNormimport pyartRADAR_NAME = 'data/sgpcsaprsurcmacI7.c0.20110520.095101.nc'# read in the dataradar = pyart.io.read_cfradial(RADAR_NAME)# create the plotfig = plt.figure(figsize=(5, 10))ax1 = fig.add_subplot(311)ax1.set_aspect('equal', 'box')display = pyart.graph.RadarDisplay(radar)display.plot('proc_dp_phase_shift', 0, ax=ax1, cmap="viridis",             colorbar_label='', axislabels=('','North South distance from radar (km)'))ax2 = fig.add_subplot(312)ax2.set_aspect('equal', 'box')display.plot('recalculated_diff_phase', 0, ax=ax2, vmin=0, vmax=5, cmap="viridis",             colorbar_label='', title='Specific differential phase', axislabels=('','North South distance from radar (km)'))ax3 = fig.add_subplot(313)display = pyart.graph.RadarDisplay(radar)display.plot_ray('proc_dp_phase_shift', 225, format_str='b-',                 axislabels_flag=False, title='Elevation: 0.8 Azimuth: 225', ax=ax3)ax3.set_ylim(0, 300)ax3.set_ylabel('Differential phase shift (degrees)')ax3.set_xlabel('Range (km)')ax4 = ax3.twinx()display.plot_ray('recalculated_diff_phase', 225, format_str='g-',                 axislabels_flag=False, title_flag=False, ax=ax4)ax4.set_ylim(-1, 7)ax4.set_ylabel('Specific differential phase (degree/km)')ax3.legend(display.plots, ["PHIDP", "KDP"], loc='upper left')plt.suptitle('PHIDP and KDP', fontsize=16)plt.show()

    図は、上から ΦDP\Phi_{\rm DP} (degree)、KDPK_{\rm DP} (degree/km)、仰角 225 度の range 方向の ΦDP\Phi_{\rm DP} と KDPK_{\rm DP} の値をそれぞれ示す。上図と中央図とを見比べると、ΦDP\Phi_{\rm DP} の位相が 50〜350 degree に大きく変化している部分で、値の大きい KDPK_{\rm DP} が計算されている。方位角 225 度のみ取り出し、レーダーの中心から range 方向の値の変化を見てみる (下図) と、確かに ΦDP\Phi_{\rm DP} の勾配と KDPK_{\rm DP} の値の変化とが一致している。値の変化が一致しているのは分かるけれども、KDPK_{\rm DP} の算出が適切かどうか判断する術は無いだろうか。次節でその一例を示す。

    ΦDP\Phi_{\rm DP} 逆算

    KDPK_{\rm DP} は、定義式から分かるように積分すると数学的には ΦDP\Phi_{\rm DP} になる。ただし、KDPK_{\rm DP} 推定時に ΦDP\Phi_{\rm DP} にフィルターまたは平滑化処理をしているため、完全一致しない。観測された ΦDP\Phi_{\rm DP} と KDPK_{\rm DP} 推定後の値がどの程度一致しているかについては、ΦDP\Phi_{\rm DP} に値を揃えたほうが分かりやすいようだ。ということで、ΦDP\Phi_{\rm DP} の逆算は ΦDPret=2∫rminrmaxKDPdr\Phi_{\rm DP} ret = 2 \displaystyle\int_{rmin}^{rmax} K_{\rm DP} dr で得る。ここで、dr は ray 方向の gate 幅、rmin rmax はそれぞれ ray 方向の初めの gate 位置・終端 gate 位置。
    import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.colors import ListedColormap, BoundaryNormimport pyartfrom pyart.config import get_field_name, get_metadataRADAR_NAME = 'data/sgpcsaprsurcmacI7.c0.20110520.095101.nc'# read in the dataradar = pyart.io.read_cfradial(RADAR_NAME)# reconstract PHIDP from KDPdr = np.diff(radar.range['data'], n=1)drm = dr/1000.pdp_field = get_field_name('differential_phase')pdp_o = radar.fields['proc_dp_phase_shift']['data']nrays, nbins = np.shape(pdp_o)rcpdp = np.ma.zeros((nrays, nbins))kdp_field = radar.fields['recalculated_diff_phase']['data'].copy()for ray in range(nrays):    kdp_ray = kdp_field[ray, :-1]    rcpdp[ray, :-1] = np.ma.cumsum(kdp_ray)*2.*drmrcpdp_dict = get_metadata(pdp_field)rcpdp_dict['data'] = rcpdpradar.add_field('reconstrated_diff_phase', rcpdp_dict)# create the plotfig = plt.figure(figsize=(6, 5))ax = fig.add_subplot(111)display = pyart.graph.RadarDisplay(radar)display.plot_ray('proc_dp_phase_shift', 225, format_str='b-',                 axislabels_flag=False, ax=ax)display.plot_ray('reconstrated_diff_phase', 225, format_str='g-',                 axislabels_flag=False, title_flag=False, ax=ax)ax.set_ylim(0, 300)ax.set_ylabel('Differential phase shift (degrees)')ax.set_xlabel('Range (km)')ax.legend(display.plots, ["PHIDP", "Reconstracted PHIDP"], loc='upper left')plt.show()

    図は、仰角 225 度の range 方向の ΦDP\Phi_{\rm DP} の値を示す。青実線は元の ΦDP\Phi_{\rm DP}、緑実線は KDPK_{\rm DP} を range 方向に積分して ΦDP\Phi_{\rm DP} に直したもの。青実線と緑実線とがほぼ一致しており、計算された KDPK_{\rm DP} は元の ΦDP\Phi_{\rm DP} によく追従出来ている。

    ここまでで、減衰補正済みの反射強度と高精度な KDPK_{\rm DP} の準備が出来た。次は、これらを用いた降水強度推定について述べる。

    R(Zc)R(Z_c) R(KDP)R(K_{\rm DP}) 算出

    Z-R 関係を用いた降水強度算出は、単偏波の場合に広く用いられてきた手法である。日本では主に R(Z)=(1/200)Z1/1.6R(Z) = (1/200) Z^{\rm 1/1.6} (eq. 10.13: Ryzhkov and Zrnic, 2019) が用いられている。二重偏波の場合は、反射強度の降雨に対する減衰補正が可能。今回は、降雨減衰補正済の反射強度を用いた。計算には est_rain_rate_z 関数を用いた。

    KDPK_{\rm DP} を用いた降水強度算出は、最近になって優位性が確立されてきたようだ (Adachi et al. 2015)。日本では R(KDP)=29 KDP0.85R(K_{\rm DP}) = {\rm 29} \, K_{\rm DP}^{\rm 0.85} が用いられている模様 (table 10.3: Ryzhkov and Zrnic, 2019)。計算には est_rain_rate_kdp 関数を用いた。
    import matplotlib.pyplot as pltfrom matplotlib.colors import ListedColormap, BoundaryNormimport pyartRADAR_NAME = 'data/sgpcsaprsurcmacI7.c0.20110520.095101.nc'# read in the dataradar = pyart.io.read_cfradial(RADAR_NAME)# remove existing correctionsradar.fields.pop('specific_attenuation')radar.fields.pop('corrected_reflectivity_horizontal')# perform attenuation correctionspec_at, cor_z = pyart.correct.calculate_attenuation(    radar, 0, refl_field='reflectivity_horizontal',    ncp_field='norm_coherent_power', rhv_field='copol_coeff',    phidp_field='proc_dp_phase_shift')radar.add_field('specific_attenuation', spec_at)radar.add_field('corrected_reflectivity_horizontal', cor_z)# calculate R(Zc) R(KDP)rzc = pyart.retrieve.est_rain_rate_z(radar, alpha=0.005, beta=0.625, refl_field='corrected_reflectivity_horizontal')rkdp = pyart.retrieve.est_rain_rate_z(radar, alpha=29, beta=0.85, kdp_field='recalculated_diff_phase')radar.add_field('RZC', rzc)radar.add_field('RKDP', rkdp)# create colormapcmap = ListedColormap(["#0FF", "#07C", "#036", "#00F", "#097", "#0A4", "#0D1", "#4F1", "#7C0", "#FF0", "#F80", "#F8F", "#F4A", "#F00"])bounds = [23, 28, 33, 37, 40, 42, 45, 47, 49, 50, 51, 52, 53]norm = BoundaryNorm(bounds, cmap.N)bounds1 = [1, 2, 4, 8, 12, 16, 24, 32, 40, 48, 58, 64, 80]norm1 = BoundaryNorm(bounds1, cmap.N)# create the plotfig = plt.figure(figsize=(5, 10))ax1 = fig.add_subplot(311)ax1.set_aspect('equal', 'box')display = pyart.graph.RadarDisplay(radar)display.plot('corrected_reflectivity_horizontal', 0, ax=ax1,              vmin=min(bounds), vmax=max(bounds), cmap=cmap, norm=norm, ticks=bounds,             colorbar_label='', axislabels=('','North South distance from radar (km)'))ax2 = fig.add_subplot(312)ax2.set_aspect('equal', 'box')display.plot('RZC', 0, ax=ax2,             vmin=min(bounds1), vmax=max(bounds1), cmap=cmap, norm=norm1, ticks=bounds1,             colorbar_label='', title='R(Zc)', axislabels=('','North South distance from radar (km)'))ax3 = fig.add_subplot(313)ax3.set_aspect('equal', 'box')display = pyart.graph.RadarDisplay(radar)display.plot('RKDP', 0, ax=ax3,             vmin=min(bounds1), vmax=max(bounds1), cmap=cmap, norm=norm1, ticks=bounds1,             colorbar_label='', title='R(KDP)')plt.suptitle('Zc vs. R(Zc) vs. R(KDP)', fontsize=16)plt.show()

    図は、上から降雨減衰補正済の反射強度 (dBZ)、R(Zc)R(Z_c) (mm/hr) R(KDP)R(K_{\rm DP}) (mm/hr) をそれぞれ示す。R(Zc)R(Z_c) と R(KDP)R(K_{\rm DP}) の図を見比べると、ほぼ同じような分布となっている。仮に、R(Zc)R(Z_c) で反射強度に対して降雨減衰補正をしなかった場合は、もう少し値が下がると考えられる。上記で得られた降水強度の妥当性は、最終的には地上雨量計との比較が必要であろう。

    おわりに

    二重偏波レーダー変数の算出・補正について、Py-ART を用いつつ具体例を示した。二重偏波レーダーでは、レーダー変数が独立変数ではなく他の変数と関連しているため、それらの処理を遡って把握する必要がある。ひとたび処理の過程を把握し適切な補正済データが得られたならば、そのデータを用い新たな知見の創出が可能となるであろう。そのためにも二重偏波レーダーは事前の準備がとても大切なのである。

    参考文献

    更新履歴

    ]]>     @@ -342,7 +342,7 @@ https://takashiunuma.github.io/blog/2021/08/19/ 2021-08-18T21:17:37.000Z - 2024-01-20T12:36:41.242Z + 2024-04-15T14:13:04.492Z はじめに

    降水強度のKDP推定アルゴリズムについてでは、主に Reimel and Kumjian (2021) で用いられた偏波間位相差変化率 (KDPK_{\rm DP}) 推定アルゴリズムについて概観した。各アルゴリズムの内容が分かってきたので、どのような分布が得られ処理前後でどのように変化するかを実際に手を動かしながら確認してみる。今回は、Py-ART にお世話になった。

    比較方法

    下記スクリプトの 10, 12 行目の内、いづれかをコメントアウトし実行する。デフォルトは kdp_maesaka を用いて処理・描画する。

    import numpy as npimport matplotlib.pyplot as pltimport pyart# Read sample data and select sweep id = 0radar = pyart.io.read_mdv('095636.mdv')radar = radar.extract_sweeps([0])# phase_proc_lp#phidp, kdp = pyart.correct.phase_proc_lp(radar, 0.0, debug=True, LP_solver='cvxopt')# kdp_maesakakdp, phidpf, phidp = pyart.retrieve.kdp_maesaka(radar, gatefilter=None, debug=True)radar.add_field('corrected_differential_phase', phidp)radar.add_field('corrected_specific_diff_phase', kdp)# create a plot of the various differential phase fieldsdisplay = pyart.graph.RadarDisplay(radar)fig = plt.figure(figsize=(10, 10))ax1 = fig.add_subplot(221)display.plot('differential_phase', 0, ax=ax1,             title='Raw Differential Phase', colorbar_label='',             axislabels_flag=False)ax2 = fig.add_subplot(222)display.plot('corrected_differential_phase', 0, ax=ax2,             title='Processed Differential Phase', colorbar_label='',             axislabels_flag=False)ax3 = fig.add_subplot(223)display.plot('specific_differential_phase', 0, ax=ax3, vmin=0, vmax=6,             title='Raw Specific Differential Phase', colorbar_label='',             axislabels_flag=False)ax4 = fig.add_subplot(224)display.plot('corrected_specific_diff_phase', 0, ax=ax4, vmin=0, vmax=6,             title='Processed Specific Differential Phase',             colorbar_label='',             axislabels_flag=False)plt.show()# plot a fields from a single raydisplay = pyart.graph.RadarDisplay(radar)fig = plt.figure(figsize=[10, 5])ax = fig.add_subplot(111)ray_num = 191# filtered phidp and original phidpdisplay.plot_ray('corrected_differential_phase', ray_num, format_str='b-',                 axislabels_flag=False, title_flag=False, ax=ax)display.plot_ray('differential_phase', ray_num, format_str='g-',                 axislabels_flag=False, title_flag=False, ax=ax)# set labelsax.set_ylim(-200, 200)ax.set_ylabel('Differential phase shift (degrees)')ax.set_xlabel('Range (km)')# plot KDP on second axisax2 = ax.twinx()display.plot_ray('corrected_specific_diff_phase', ray_num, format_str='r-',                 axislabels_flag=False, title_flag=False, ax=ax2)ax2.set_ylim(-2, 8)ax2.yaxis.grid(color='gray', linestyle='dashed')ax.legend(display.plots,          ["Filtered phiDP", "Original phiDP", 'KDP'],          loc='upper right')plt.show()

    結果

    今回は PHASE_PROC_LP と KDP_MAESAKA を例に、 PPI 水平分布と ray 方向のグラフとをそれぞれ比較してみる。

    PPI 水平分布

    まず、PPI 水平分布を確認する。各図の説明は下記。

    PHASE_PROC_LP の処理済み ΦDP\Phi_{\rm DP} は元の ΦDP\Phi_{\rm DP} と比べてかなり平滑化されている。横軸で 0〜100 km の範囲でややまだら模様になっている部分がほぼ 0 degree になっている。一方、KDP_MAESAKA では、処理済み ΦDP\Phi_{\rm DP} は元データである程度滑らかなところはそのまま利用し、まだら模様になっている部分を落としているようだ (-50 degree 以下の色の階調が緑に変化している部分)。KDPK_{\rm DP} の分布も ΦDP\Phi_{\rm DP} の処理具合に依存し、PHASE_PROC_LP の処理済み KDPK_{\rm DP} は KDP_MAESAKA よりも PHASE_PROC_LP で分布としては滑らかに見える。

    PHASE_PROC_LP

    KDP_MAESAKA

    ray 方向のグラフ

    次に、ray 方向のグラフを確認する。各実線の説明は下記。

    先ほど PPI で確認したように、処理済み ΦDP\Phi_{\rm DP} はどちらの場合も非常に滑らかな分布となっている (青実線)。値の下限が 0 とそうでない場合になっているけれども、値域の 0〜360 と -180〜180 の違いであるため値の変動を見ることが本質。一方、処理済み KDPK_{\rm DP} (赤実線) を見ると特徴が異なっている。PHASE_PROC_LP の場合は処理済み ΦDP\Phi_{\rm DP} の勾配に応じた値の変化をしている。KDP_MAESAKA の場合は処理済み ΦDP\Phi_{\rm DP} が滑らかでも、元の ΦDP\Phi_{\rm DP} (緑実線) にある細かな勾配も値の変化として表現されている。

    PHASE_PROC_LP

    KDP_MAESAKA

    考察

    PHASE_PROC_LP では ΦDP\Phi_{\rm DP} の平滑化処理がかなり強めのようだ。その平滑化済の ΦDP\Phi_{\rm DP} から有限差分で KDPK_{\rm DP} を計算するため、水平分布・ray 方向のグラフで見ても非常に滑らか分布・値を示しているようだ。滑らかな分布が得られる一方で、局所性を全て棄却しているので KDPK_{\rm DP} の値には必然的に上限ができていそうだ。
    KDP_MAESAKA では、フィルターによる平滑化処理というよりはフィルターの重みで平滑化度合いが決まっているようだ。加えて、cost function が最小となるように平滑化した ΦDP\Phi_{\rm DP} を元の値に近づける。このため、たとえ処理済み ΦDP\Phi_{\rm DP} で平滑化されていても、元の ΦDP\Phi_{\rm DP} の変化を反映出来るという訳だ。ただ、元の ΦDP\Phi_{\rm DP} を適切に処理しておかないと、KDPK_{\rm DP} 計算時に偽の勾配を表してしまう。KDP_MAESAKA の PPI 水平分布で KDPK_{\rm DP} の値が大きい部分が散在しているのは、こういった部分が反映されているものと考えられる。このように、元の ΦDP\Phi_{\rm DP} に対するフィルター・平滑化処理のかけ方が重要となるのは、どの推定手法でも共通のようだ。

    おわりに

    Py-ART のサンプルデータ・サンプルスクリプトを用い、幾つかのKDP推定アルゴリズムについて実際に手を動かしながら確認してみた。図や式を眺めていたときとは異なり、計算時の課題やアルゴリズムの具体的な特徴を把握することが出来た。

    参考文献

    ]]>     @@ -388,7 +388,7 @@ https://takashiunuma.github.io/blog/2021/07/28/ 2021-07-27T21:47:39.000Z - 2024-01-20T10:07:05.377Z + 2024-04-15T14:13:04.504Z はじめに

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

    用語

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

    NWS OPERATIONAL

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

    HUBBERT AND BRINGI FIR FILTER

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

    CALC_KDP_BRINGI (CSU_RADAR_TOOLS)

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

    PHASE_PROC_LP (Py-ART)

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

    KDP_SCHNEEBELI (Py-ART)

    ΦDP\Phi_{\rm DP} の処理に Kalman filter (KF) を用い、その結果から KDPK_{\rm DP} を算出する手法。カルマンフィルタでは、処理前のΦDP\Phi_{\rm DP}データが観測値として用いられ、KDPK_{\rm DP}, ΦDP\Phi_{\rm DP}, δ\delta間の所定の関係が、観測値から場の変数への遷移に用いられる。カルマンフィルタはそれぞれの gate に対して 2 回行われる。1 回目はレーダーから遠ざかる方向(前方処理)、2 回目は観測範囲終端からレーダーに向う方向(後方処理)。このプロセスはそれぞれの gate に対して複数回繰り返される。各々の実行時には、固有の共分散行列の値が含まれ、個々のカルマンフィルタのアンサンブルを推定する。それぞれの gate における最終的な KDPK_{\rm DP} 値は、前方・後方処理のアンサンブル平均、もしくは 2 つのカルマンフィルタ推定値の平均が適用される。このとき、アンサンブル平均が gate 間で 0.15∘ km−10.15^{\circ}\, {\rm km}^{-1} を上回り増加する(下回り減少する)場合は、前方(後方)処理によるカルマンフィルタ推定値が用いられる。勾配がある場所で−0.15∘ km−1-0.15^{\circ}\, {\rm km}^{-1}と+0.15∘ km−1+0.15^{\circ}\, {\rm km}^{-1}の間の値をとる場合は、前方・後方処理双方によるカルマンフィルタ推定値の平均が用いられる。アンサンブルサイズはKDPK_{\rm DP}の標準偏差に依存する。KDPK_{\rm DP}の標準偏差が大きいほどアンサンブルサイズの大きな値が最終推定値に用いられる。

    KDP_VULPIANI (Py-ART)

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

    KDP_MAESAKA (Py-ART)

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

    各アルゴリズムの比較

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

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

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

    # # plot_ray.py#import numpy as npimport matplotlib.pyplot as pltimport 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 scanradar = radar.extract_sweeps([0])# Set a field for a single raydisplay = pyart.graph.RadarDisplay(radar)fig = plt.figure(figsize=[10, 5])ax = fig.add_subplot(111)# Set a ray numberray_num = 200# Plot original phidpdisplay.plot_ray('differential_phase', ray_num, format_str='b-',                 axislabels_flag=False, title_flag=False, ax=ax)# Set labelsax.set_ylim(-180, 180)ax.set_ylabel('Differential phase shift (degrees)')ax.set_xlabel('Range (km)')# Plot KDP on second axisax2 = 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)# Decorateax2.yaxis.grid(color='gray', linestyle='dashed')ax.legend(display.plots, ["phiDP", 'KDP'], loc='upper right')# Displayplt.show()

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

    参考文献

    更新履歴

    ]]>     @@ -411,7 +411,7 @@ https://takashiunuma.github.io/blog/2021/04/28/ 2021-04-28T13:47:04.000Z - 2024-01-20T10:02:53.187Z + 2024-04-15T14:13:04.503Z はじめに

    以前、Signal Processing Methods という記事を書いた。言葉だけだとどうもイメージが沸かないので、実際のデータで図を描いてみることにした。最近では Python 関連の便利でたくさんのパッケージが出回っている。いろいろと試した結果、ここでは ω\omegaradlib を用いることにした。Python や ω\omegaradlib の環境構築については、こちらを参照されたい。

    降水・非降水エコー

    ω\omegaradlib で Simple fuzzy echo classification (ファジーロジックによるエコー判別) が出来るので、一例を示す。ここに書いてあるとおりにコマンドを打つと、下記の図が出来る。

    左図は水平偏波反射強度 (dBZ)、右図は判別結果。エコー判別は、0 で降水エコー 1 で非降水エコーをそれぞれ示す。例えば、サイト中心付近の弱いエコー域は判別結果だと 1 に近いので非降水エコーと見なせる、らしい。これだとブラックボックス過ぎてよく分からないので、幾つかの要素を詳しく観てみる。

    偏波間位相差 (ψdp\psi_{dp})

    左図は上図と同じ、右図は偏波間位相差。偏波間位相差は、水平偏波の位相変化と垂直偏波の位相変化との差。受信電力から求めた場合は受信電力偏波間位相差と言うらしい。今回示した図は、受信電力偏波間位相差からノイズを取り除き、スムージングを掛け、正常値以外の領域では無効値になっている模様。このため、値域は -180度〜180度に収まっているようだ。粒経の大きい扁平な雨滴が観測された場合、垂直偏波よりも水平偏波の位相が遅れる。これを利用して降水強度推定にも利用される。判別結果と偏波間位相差とを見比べると、偏波間位相差の値が算出されている領域はほぼ降水エコーとなっているようだ。偏波間位相差の有効・無効も影響しているかもしれない。

    偏波間相関係数 (ρ\rhohv)

    右図のみ偏波間相関係数に変更。偏波間相関係数は、水平偏波と水平偏波との信号の相関係数。ある観測領域内の散乱体について一様性の度合いを表す。値域は 0〜1 をとり、1 に近いほど一様性が高い (雨粒)。雹や融解層内の雨滴は一様性がやや低い。また、生物やクラッタの場合は一様性が低い傾向にある。判別結果と偏波間相関係数とを見比べると、サイト中心付近で非降水エコーと判別された部分は、なるほど偏波間相関係数の値が低い。一方、水平偏波反射強度で比較的値の大きいエコー域でかつ偏波間相関係数の大きい領域は、降水エコーとして判別されている。

    反射因子差 (ZDRZ_{DR})

    右図のみ反射因子差に変更。反射因子差は、水平偏波と垂直偏波の反射因子の比 (或いは、強度 (dBZ) の差) の値で、散乱体の縦横比を表す。典型的な値域は -4〜10 dB。球形に近いほど値は 0 に近く、扁平な形状ほど値が大きくなる。反射因子差の値を観てみると、サイト中心の南側にあるエコー域で負の値と正の値とが混在している。この領域は、偏波間相関係数の値も低く、判別結果も非降水エコーになっている。

    以上、サンプルデータで概観してみたものの、それぞれの特徴を観る上で十分な情報量が含まれていたように思う。続いて、適切な反射強度・降水強度算出に欠かせない、降雨減衰補正について述べる。

    降雨減衰補正

    レーダービームが強雨域をレンジ方向に通過した場合、強雨域のレンジ後方で送信電力(レーダー側では受信電力)がレーダーの波長によっては著しく減衰してしまう。この減衰量を推定し、降雨による減衰としてを補正するのが降雨減衰補正である。ここに書いてあるコマンドを少し変更し、降雨減衰補正前・降雨減衰補正後の図をそれぞれ作成した。ドイツ気象局管轄の Feldberg と Tuerkheim というレーダーサイトの合成図になっている。

    降雨減衰補正後のデータは、降雨減衰補正前に比べ緑〜黄色 (例えば 100 mm/h 以上) の領域が広くなっており、適切に補正されている、らしい。パラメータは恐らくこの地域にあわせてあるのかもしれない。この例では、減衰補正の手法として Modified Kraemer algorithm が用いられているようだ。補正方法は他にもあり、オプションとして用意されている他、各国でやり方も異なるようだ。

    おわりに

    初学者にも使いやすいライブラリを提供してくださっている ω\omegaradlib プロジェクトやその関連プロジェクトに関わる方々に対し謝意を示す。

    参考

    更新履歴

    ]]>     diff --git a/categories/index.html b/categories/index.html index 7be5be2..eb72c7c 100644 --- a/categories/index.html +++ b/categories/index.html @@ -128,7 +128,7 @@ - + diff --git a/index.html b/index.html index 0713e1e..ee5d2ac 100644 --- a/index.html +++ b/index.html @@ -128,7 +128,7 @@ - + diff --git a/page/2/index.html b/page/2/index.html index d9fdc40..7f64109 100644 --- a/page/2/index.html +++ b/page/2/index.html @@ -128,7 +128,7 @@ - + diff --git a/page/3/index.html b/page/3/index.html index 2623210..1c73ff8 100644 --- a/page/3/index.html +++ b/page/3/index.html @@ -128,7 +128,7 @@ - + diff --git a/page/4/index.html b/page/4/index.html index 9968829..4517730 100644 --- a/page/4/index.html +++ b/page/4/index.html @@ -128,7 +128,7 @@ - + diff --git a/page/5/index.html b/page/5/index.html index ea9f684..df493a4 100644 --- a/page/5/index.html +++ b/page/5/index.html @@ -128,7 +128,7 @@ - + diff --git a/rss2.xml b/rss2.xml index 57f60d4..c1b1c7e 100644 --- a/rss2.xml +++ b/rss2.xml @@ -15,7 +15,7 @@ - Sat, 09 Mar 2024 11:08:09 GMT + Mon, 15 Apr 2024 14:13:04 GMT http://hexo.io/ @@ -82,7 +82,7 @@ - はじめに

    これまでに、粒径分布とレーダー変数や粒径分布パラメータのレーダー変数による推定といった記事で、主に観測データに関連した話題を扱ってきた。観測データはスナップショットではあるものの、真値に近い状態を示す。一方、時間・空間分解能に制約があるために観測結果がどうしてそうなったのか、という要因の特定には不向きな場合がある。その点を補うものとして、数値モデルによる素過程のシミュレーションや解析が挙げられる。中でも雲微物理過程は、雨や雲を扱う重要な物理過程である。計算コストの観点から、粒径分布は特定の関数で近似されたものが使われることが多い。一方、Spectral bin microphysics (SBM) model は、粒径毎の数濃度等を陽に解くものである。最近では、Hebrew University で開発された Cloud Model (Khain et al. 2004, 2008; HUCM) の spectral (bin) microphysics 部分を他の大気モデルへ移植したものが、利用可能になっている (e.g., Khain et al. 2011; Igel and van den Heever 2017)。本稿では、粒径分布の特徴を把握するためのツールの一つとして、WRF-SBM の利用と簡単な出力結果の描画を試みる。

    WRF model

    Weather Research and Frecasting Model (WRF) は、主に米国大気研究センターが開発している領域数値気象予報モデルである。WRF のバージョンは、2023年12月22日時点で 4.5.2 が最新版である (参考)。HUCM SBM は mp_physics = 32 or 30 として実装されている。これらは、それぞれ FULL-SBM、FAST-SBM と呼ばれている。前者は計算時間を要するものの物理過程の厳密さを重視したバージョン、後者は前者を特定の現象に対してやや簡略化し計算コストを減少させたバージョンである。バージョン 4.4.1 以降、FAST-SBM のみ利用可能であり、FULL-SBM は利用不可となっている (参考)。ここでは、計算速度よりも物理過程の正確さを重視することを目的とし、FULL-SBM が利用可能なバージョン 4.2.2 を用いることにした (参考)。

    動作環境

    • Rocky Linux 8.9 x86_64
    • GNU compiler version 8.5.0
    • NetCDF-C version 4.7.0
    • NetCDF-Fortran version 4.5.2
    • HDF5 version 1.10.5
    • ZLIB version 1.2.11
    • WRF version 4.2.2

    必要なライブラリは、dnf install した。

    コンパイル

    • export NETCDF=/usr && export HDF5=/usr
      • 事前に環境変数の指定が必要。
      • 場合によっては、 export NETCDF_classic=1 も必要。
    • ./configure
      • 34 番の GNU compiler smpar (OpenMP) を選択。

    configure.wrf の LIB_EXTERNAL に -L/usr/lib64 を追記し、コンパイル。

    • ./compile em_quarter_ss

    main/ 以下に ideal.exe と wrf.exe が生成されていれば、コンパイル成功。

    出力変数の追加

    今回は、粒径分布の特徴を把握することを意図し、各 bin の変数を出力できるように修正を施す。出力変数を追加するには、Registry を編集するのが一番簡単。方法は下記。

    • Registry/registry.sbm の編集
      • ff?i?? を出力する。下2桁が bin 番号になっていて、bin 毎に変数が出力される。
      • 各変数の h3rusdf=(bdy_interp:dt) の部分を hrusdf=(bdy_interp:dt) にすると良い。
        • 元の h3 の 3 は I/O stream 番号 3、つまり Fortran でいうところの出力番号に相当。
        • WRF では、主出力ファイルとは別名で変数を出力するオプションがある。
        • 詳細はこちらを参照。
    • 再コンパイル
      • Registry を変更した場合、反映させるためには再コンパイルが必要。
      • ./clean && ./compile em_quarter_ss

    定数ファイルの準備

    WRF で SBM を実行するには、追加の定数ファイルが必要となる。下記に従い、定数ファイルをダウンロード・展開・配置する。

    • こちらから3つのファイルをダウンロードし、test/em_quarter_ss/ 以下にそれぞれ展開・配置する。
      • SBM_input_33.tgz
      • SBM_input_43.tgz
      • SBM_scatter_amplit.tgz

    実行

    準備が出来たら、namelist.input を編集し、ideal.exe、wrf.exe の順に実行する。
    今回は、ideal case = 2 の supercell 理想化実験を用いる。namelist.input で mp_physics = 32 に変更し、雲微物理過程のみ用いる。計算領域は、mass point で 41x41x40 grid (水平方向 82 km、鉛直方向に 20 km)。モデル上端から 5 km までダンピング層を配置。ideal.exe を実行すると wrfinput_d01 が出力され、その後 wrf.exe を実行すると wrfout_d01_* が出力される。

    データの描画

    まず、初期プロファイルの input_sounding データを Skew-T LogP plot by MetPy の方法で描画する。ただし、read_fwf の引数に infer_nrows=222 を追加する。理由は、行数のデフォルト値が 100 に制限されているため。

    横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。Skew-T LogP plot by MetPy で示した例よりも赤色の陰影で示される面積が大きく、また、中・上層は乾燥している点も異なる。

    次に、WRF の出力結果 wrfout_d01_* を wrf-python により描画する。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import (to_np, getvar)# WRF 出力ファイル指定ncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# QRAIN 変数を取得qr = getvar(ncfile, "QRAIN")# 描画開始 (xy 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[0,15:25,15:25].west_east)+15,             to_np(qr[0,15:25,15:25].south_north)+15,             to_np(qr[0,15:25,15:25]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title("QRAIN (k = "+str(kk)+")")plt.xlabel('west_east grid number')plt.ylabel('south_north grid number')plt.grid()plt.show()# 描画開始 (xz 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[:15,22,16:23].west_east)+16,             to_np(qr[:15,22,16:23].bottom_top),             to_np(qr[:15,22,16:23]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title(qr.description + " (" + qr.units + ")")plt.xlabel('west_east grid number')plt.ylabel('bottom_top grid number')plt.grid()plt.show()

    計算開始 50 分後の雨水混合比の水平分布。グリッド位置 (19,21) で雨水混合比の極大値が計算されている。

    雨水混合比の y = 22 における x-z 断面。地上付近に計算されていた雨水混合比は、鉛直方向には一様ではなく高度とともに値が変化している。

    以後、QRAIN の極大値が計算されていたグリッド近傍の (19,22) に着目する。

    出力結果から粒径分布を計算

    方針は下記の通り。

    • 定数ファイル読み込み
      • 粒径 (bulkradii.asc_s_0_03_0_9)
      • 密度 (bulkdens.asc_s_0_03_0_9)
      • 質量 (masses.asc)
    • WRF の出力結果読み込み
      • wrf-python で ff1i01-ff1i33 を読み込み。
      • xarray の dataarray をまとめて dataset にする。
      • 粒径分布の計算し、粒径分布の描画。今回は、特定 grid の粒径分布を鉛直方向 4 grid 分を図示。

    ソースコードは下記。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import getvar# 粒径 (半径) [unit: cm] 定数データ読み込み# index とカテゴリの対応は下記。#  0: aerosol#  1: cloud + rain#  2: ice + snow#  3: graupel + hail#  4: snow#  5: graupel#  6: hailbulkradii = [[0] * 33 for i in range(7)]with open("bulkradii.asc_s_0_03_0_9") as f1:    k = 0    j = 0    for line in f1:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkradii[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 密度 [unit: g cm^-3] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。bulkdens = [[0] * 33 for i in range(7)]with open("bulkdens.asc_s_0_03_0_9") as f2:    k = 0    j = 0    for line in f2:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkdens[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 質量 [unit: g] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。masses = [[0] * 33 for i in range(7)]with open("masses.asc") as f3:    k = 0    j = 0    for line in f3:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                masses[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# WRF 出力結果を読み込みncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# 出力結果の cloud and rain bin データを xarray dataset に変換。vlist = ['ff1i01', 'ff1i02', 'ff1i03', 'ff1i04', 'ff1i05',         'ff1i06', 'ff1i07', 'ff1i08', 'ff1i09', 'ff1i10',         'ff1i11', 'ff1i12', 'ff1i13', 'ff1i14', 'ff1i15',         'ff1i16', 'ff1i17', 'ff1i18', 'ff1i19', 'ff1i20',         'ff1i21', 'ff1i22', 'ff1i23', 'ff1i24', 'ff1i25',         'ff1i26', 'ff1i27', 'ff1i28', 'ff1i29', 'ff1i30',         'ff1i31', 'ff1i32', 'ff1i33']xdsff1i = getvar(ncfile, "ff1i01")xdsff1i = xdsff1i.to_dataset()for ivar in vlist:    xdsff1i[ivar] = getvar(ncfile, ivar)# 粒径を半径から直径に、単位を cm から mm に、それぞれ変換。diameter = np.zeros((33), float)for i in range(0,33):    diameter[i] = bulkradii[0][i] * 2. * 10. # unit: diameter in mm# 特定グリッドの指定 (水平分布で極大値のあったグリッド)jj = 22ii = 19for kk in [3, 2, 1, 0]:    # 粒径分布データ用の numpy array を初期化。    ff1i = np.zeros((33), float)    # 33 bin に分かれた変数を一つの numpy array に結合    i = 0    for ivar in vlist:        ff1i[i] = xdsff1i[ivar][kk,jj,ii].to_numpy()        i = i + 1    # 粒形分布 N(D) [m^-3 mm^-1] 計算    for i in range(0,33):        # unit: (kg kg^-1) (g^-1) (g cm^-3) (mm^-1)        #    -> (kg kg^-1) (kg^-1) (kg m^-3) (mm^-1)        #    -> m^-3 mm^-1 (多分これであっているはず)        ff1i[i] = (ff1i[i] / (masses[0][i] * 0.001)) * (bulkdens[0][i] * 1000.) / (diameter[i] * 1000.)    # 粒径分布データ描画。ただし、rain bin のみ。    plt.plot(diameter[16:33], ff1i[16:33], label="k = "+str(kk))    plt.scatter(diameter[16:33], ff1i[16:33])plt.yscale('log')plt.xlim(0,6)plt.ylim(0.01, 10000)plt.title('N(D) for rain at ('+str(ii)+','+str(jj)+')')plt.xlabel('Diameter (mm)')plt.ylabel('Number concentration (mm$^{-1}$ m$^{-3}$)')plt.grid()plt.legend()plt.show()

    粒形分布を確認すると、地上に近い (i.e., k = 0) ほど粒径 2 mm 未満の頻度が低く、上空 (i.e., k = 3) ほど頻度が高くなっている。一方、粒径 2 mm よりも大きい粒径は地上ほど頻度が高い。また、頻度のピークは、地上に近いほど複数 (粒径 0.26 mm, 0.51 mm, 1.62 mm) 現れる傾向にある。

    鉛直分布を確認した時、雨水混合比の値は k = 8 を極大とし、地上に向かって減少していた。これは、粒形分布では第一近似として粒径 2 mm 未満の頻度が低くなっていたことに対応していると考えられる。雨水混合比としては減少する一方、粒径 2 mm 以上の頻度が僅かながら高くなっていた。このことは、雨滴同士の衝突併合によって雨粒が地上に向かって成長している可能性を示唆する。このように、粒形分布を確認することで雨の特徴をより多角的に捉えることが可能となる。

    おわりに

    今回は、WRF-SBM を用いた粒形分布の解析例を示した。実は、十数年前に WRF モデルの version 3 系列を利用したことがあった。その当時よりも利用可能なオプションが格段に増えており、使い勝手もとてもよくなっていると感じた。特に、描画を Python で行える点は、近年の流行を追っていて舌を巻く。当時は大雨を再現するだけでなく出力結果の解析にも時間と労力を費やしていたが、今後はそんな苦労も少なくなるだろうか。

    参考文献

    • wrf-model: https://github.com/wrf-model/WRF
    • wrf-python: https://github.com/NCAR/wrf-python
    • Khain, A. P., A. Pokrovsky, M. Pinsky, A. Seifert, and V. Phillips, 2004: Effects of atmospheric aerosols on deep convective clouds as seen from simulations using a spectral microphysics mixed-phase cumulus cloud model. Part I: Model description. J. Atmos. Sci., 61, 2963–2982.
    • Khain, A. P., N. Benmoshe, and A. Pokrovsky, 2008: Factors determining the impact of aerosols on surface precipitation from clouds: An attempt at classification. J. Atmos. Sci., 65, 1721–1748.
    • Khain, A. P., D. Rosenfeld, A. Pokrovsky, U. Blahak, and A. Ryzhkov, 2011: The role of CCN in precipitation and hail in a midlatitude storm as seen in simulations using a spectral (bin) microphysics model in a 2D dynamic frame. Atmos. Res., 99, 129–146.
    • Igel, A. L., and S. C. van den Heever, 2017: The role of the gamma function shape parameter in determining differences between condensation rates in bin and bulk microphysics schemes. Atmos. Chem. Phys., 17, 4599–4609.

    更新履歴

    • 2024-01-19: 初稿
    • 2024-01-20: 誤字を修正
    • 2024-01-23: 表記を修正
    • 2024-02-01: Registory の URL を修正
    ]]>     + はじめに

    これまでに、粒径分布とレーダー変数や粒径分布パラメータのレーダー変数による推定といった記事で、主に観測データに関連した話題を扱ってきた。観測データはスナップショットではあるものの、真値に近い状態を示す。一方、時間・空間分解能に制約があるために観測結果がどうしてそうなったのか、という要因の特定には不向きな場合がある。その点を補うものとして、数値モデルによる素過程のシミュレーションや解析が挙げられる。中でも雲微物理過程は、雨や雲を扱う重要な物理過程である。計算コストの観点から、粒径分布は特定の関数で近似されたものが使われることが多い。一方、Spectral bin microphysics (SBM) model は、粒径毎の数濃度等を陽に解くものである。最近では、Hebrew University で開発された Cloud Model (Khain et al. 2004, 2008; HUCM) の spectral (bin) microphysics 部分を他の大気モデルへ移植したものが、利用可能になっている (e.g., Khain et al. 2011; Igel and van den Heever 2017)。本稿では、粒径分布の特徴を把握するためのツールの一つとして、WRF-SBM の利用と簡単な出力結果の描画を試みる。

    WRF model

    Weather Research and Frecasting Model (WRF) は、主に米国大気研究センターが開発している領域数値気象予報モデルである。WRF のバージョンは、2023年12月22日時点で 4.5.2 が最新版である (参考)。HUCM SBM は mp_physics = 32 or 30 として実装されている。これらは、それぞれ FULL-SBM、FAST-SBM と呼ばれている。前者は計算時間を要するものの物理過程の厳密さを重視したバージョン、後者は前者を特定の現象に対してやや簡略化し計算コストを減少させたバージョンである。バージョン 4.4.1 以降、FAST-SBM のみ利用可能であり、FULL-SBM は利用不可となっている (参考)。ここでは、計算速度よりも物理過程の正確さを重視することを目的とし、FULL-SBM が利用可能なバージョン 4.2.2 を用いることにした (参考)。

    動作環境

    • Rocky Linux 8.9 x86_64
    • GNU compiler version 8.5.0
    • NetCDF-C version 4.7.0
    • NetCDF-Fortran version 4.5.2
    • HDF5 version 1.10.5
    • ZLIB version 1.2.11
    • WRF version 4.2.2

    必要なライブラリは、dnf install した。

    コンパイル

    • export NETCDF=/usr && export HDF5=/usr
      • 事前に環境変数の指定が必要。
      • 場合によっては、 export NETCDF_classic=1 も必要。
    • ./configure
      • 34 番の GNU compiler smpar (OpenMP) を選択。

    configure.wrf の LIB_EXTERNAL に -L/usr/lib64 を追記し、コンパイル。

    • ./compile em_quarter_ss

    main/ 以下に ideal.exe と wrf.exe が生成されていれば、コンパイル成功。

    出力変数の追加

    今回は、粒径分布の特徴を把握することを意図し、各 bin の変数を出力できるように修正を施す。出力変数を追加するには、Registry を編集するのが一番簡単。方法は下記。

    • Registry/registry.sbm の編集
      • ff?i?? を出力する。下2桁が bin 番号になっていて、bin 毎に変数が出力される。
      • 各変数の h3rusdf=(bdy_interp:dt) の部分を hrusdf=(bdy_interp:dt) にすると良い。
        • 元の h3 の 3 は I/O stream 番号 3、つまり Fortran でいうところの出力番号に相当。
        • WRF では、主出力ファイルとは別名で変数を出力するオプションがある。
        • 詳細はこちらを参照。
    • 再コンパイル
      • Registry を変更した場合、反映させるためには再コンパイルが必要。
      • ./clean && ./compile em_quarter_ss

    定数ファイルの準備

    WRF で SBM を実行するには、追加の定数ファイルが必要となる。下記に従い、定数ファイルをダウンロード・展開・配置する。

    • こちらから3つのファイルをダウンロードし、test/em_quarter_ss/ 以下にそれぞれ展開・配置する。
      • SBM_input_33.tgz
      • SBM_input_43.tgz
      • SBM_scatter_amplit.tgz

    実行

    準備が出来たら、namelist.input を編集し、ideal.exe、wrf.exe の順に実行する。
    今回は、ideal case = 2 の supercell 理想化実験を用いる。namelist.input で mp_physics = 32 に変更し、雲微物理過程のみ用いる。計算領域は、mass point で 41x41x40 grid (水平方向 82 km、鉛直方向に 20 km)。モデル上端から 5 km までダンピング層を配置。ideal.exe を実行すると wrfinput_d01 が出力され、その後 wrf.exe を実行すると wrfout_d01_* が出力される。

    データの描画

    まず、初期プロファイルの input_sounding データを Skew-T LogP plot by MetPy の方法で描画する。ただし、read_fwf の引数に infer_nrows=222 を追加する。理由は、行数のデフォルト値が 100 に制限されているため。

    横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。Skew-T LogP plot by MetPy で示した例よりも赤色の陰影で示される面積が大きく、また、中・上層は乾燥している点も異なる。

    次に、WRF の出力結果 wrfout_d01_* を wrf-python により描画する。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import (to_np, getvar)# WRF 出力ファイル指定ncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# QRAIN 変数を取得qr = getvar(ncfile, "QRAIN")# 描画開始 (xy 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[0,15:25,15:25].west_east)+15,             to_np(qr[0,15:25,15:25].south_north)+15,             to_np(qr[0,15:25,15:25]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title("QRAIN (k = "+str(kk)+")")plt.xlabel('west_east grid number')plt.ylabel('south_north grid number')plt.grid()plt.show()# 描画開始 (xz 断面)fig = plt.figure(figsize=(8,6))# 雨水混合比の描画plt.contourf(to_np(qr[:15,22,16:23].west_east)+16,             to_np(qr[:15,22,16:23].bottom_top),             to_np(qr[:15,22,16:23]),             10,             cmap=matplotlib.colormaps.get_cmap("jet"))# カラーバー追加plt.colorbar(shrink=.98)# タイトル描画plt.title(qr.description + " (" + qr.units + ")")plt.xlabel('west_east grid number')plt.ylabel('bottom_top grid number')plt.grid()plt.show()

    計算開始 50 分後の雨水混合比の水平分布。グリッド位置 (19,21) で雨水混合比の極大値が計算されている。

    雨水混合比の y = 22 における x-z 断面。地上付近に計算されていた雨水混合比は、鉛直方向には一様ではなく高度とともに値が変化している。

    以後、QRAIN の極大値が計算されていたグリッド近傍の (19,22) に着目する。

    出力結果から粒径分布を計算

    方針は下記の通り。

    • 定数ファイル読み込み
      • 粒径 (bulkradii.asc_s_0_03_0_9)
      • 密度 (bulkdens.asc_s_0_03_0_9)
      • 質量 (masses.asc)
    • WRF の出力結果読み込み
      • wrf-python で ff1i01-ff1i33 を読み込み。
      • xarray の dataarray をまとめて dataset にする。
      • 粒径分布の計算し、粒径分布の描画。今回は、特定 grid の粒径分布を鉛直方向 4 grid 分を図示。

    ソースコードは下記。

    import numpy as npfrom netCDF4 import Datasetimport matplotlibimport matplotlib.pyplot as pltfrom wrf import getvar# 粒径 (半径) [unit: cm] 定数データ読み込み# index とカテゴリの対応は下記。#  0: cloud + rain#  1: ice/columns#  2: ice/plates#  3: ice/dendrites#  4: snow#  5: graupel#  6: hailbulkradii = [[0] * 33 for i in range(7)]with open("bulkradii.asc_s_0_03_0_9") as f1:    k = 0    j = 0    for line in f1:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkradii[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 密度 [unit: g cm^-3] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。bulkdens = [[0] * 33 for i in range(7)]with open("bulkdens.asc_s_0_03_0_9") as f2:    k = 0    j = 0    for line in f2:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                bulkdens[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# 質量 [unit: g] 定数データ読み込み# index とカテゴリの対応は粒径と同じ。masses = [[0] * 33 for i in range(7)]with open("masses.asc") as f3:    k = 0    j = 0    for line in f3:        records = line.rstrip('\n')        columns = records.split(' ')[:]        for i in range(len(columns)):            if (columns[i] != ''):                masses[k][j] = float(columns[i])                j = j + 1                if (j%33 == 0):                    j = 0                    k = k + 1# WRF 出力結果を読み込みncfile = Dataset("wrfout_d01_0001-01-01_00:50:00")# 出力結果の cloud and rain bin データを xarray dataset に変換。vlist = ['ff1i01', 'ff1i02', 'ff1i03', 'ff1i04', 'ff1i05',         'ff1i06', 'ff1i07', 'ff1i08', 'ff1i09', 'ff1i10',         'ff1i11', 'ff1i12', 'ff1i13', 'ff1i14', 'ff1i15',         'ff1i16', 'ff1i17', 'ff1i18', 'ff1i19', 'ff1i20',         'ff1i21', 'ff1i22', 'ff1i23', 'ff1i24', 'ff1i25',         'ff1i26', 'ff1i27', 'ff1i28', 'ff1i29', 'ff1i30',         'ff1i31', 'ff1i32', 'ff1i33']xdsff1i = getvar(ncfile, "ff1i01")xdsff1i = xdsff1i.to_dataset()for ivar in vlist:    xdsff1i[ivar] = getvar(ncfile, ivar)# 粒径を半径から直径に、単位を cm から mm に、それぞれ変換。diameter = np.zeros((33), float)for i in range(0,33):    diameter[i] = bulkradii[0][i] * 2. * 10. # unit: diameter in mm# 特定グリッドの指定 (水平分布で極大値のあったグリッド)jj = 22ii = 19for kk in [3, 2, 1, 0]:    # 粒径分布データ用の numpy array を初期化。    ff1i = np.zeros((33), float)    # 33 bin に分かれた変数を一つの numpy array に結合    i = 0    for ivar in vlist:        ff1i[i] = xdsff1i[ivar][kk,jj,ii].to_numpy()        i = i + 1    # 粒形分布 N(D) [m^-3 mm^-1] 計算    for i in range(0,33):        # unit: (kg kg^-1) (g^-1) (g cm^-3) (mm^-1)        #    -> (kg kg^-1) (kg^-1) (kg m^-3) (mm^-1)        #    -> m^-3 mm^-1 (多分これであっているはず)        ff1i[i] = (ff1i[i] / (masses[0][i] * 0.001)) * (bulkdens[0][i] * 1000.) / (diameter[i] * 1000.)    # 粒径分布データ描画。ただし、rain bin のみ。    plt.plot(diameter[16:33], ff1i[16:33], label="k = "+str(kk))    plt.scatter(diameter[16:33], ff1i[16:33])plt.yscale('log')plt.xlim(0,6)plt.ylim(0.01, 10000)plt.title('N(D) for rain at ('+str(ii)+','+str(jj)+')')plt.xlabel('Diameter (mm)')plt.ylabel('Number concentration (mm$^{-1}$ m$^{-3}$)')plt.grid()plt.legend()plt.show()

    粒形分布を確認すると、地上に近い (i.e., k = 0) ほど粒径 2 mm 未満の頻度が低く、上空 (i.e., k = 3) ほど頻度が高くなっている。一方、粒径 2 mm よりも大きい粒径は地上ほど頻度が高い。また、頻度のピークは、地上に近いほど複数 (粒径 0.26 mm, 0.51 mm, 1.62 mm) 現れる傾向にある。

    鉛直分布を確認した時、雨水混合比の値は k = 8 を極大とし、地上に向かって減少していた。これは、粒形分布では第一近似として粒径 2 mm 未満の頻度が低くなっていたことに対応していると考えられる。雨水混合比としては減少する一方、粒径 2 mm 以上の頻度が僅かながら高くなっていた。このことは、雨滴同士の衝突併合によって雨粒が地上に向かって成長している可能性を示唆する。このように、粒形分布を確認することで雨の特徴をより多角的に捉えることが可能となる。

    おわりに

    今回は、WRF-SBM を用いた粒形分布の解析例を示した。実は、十数年前に WRF モデルの version 3 系列を利用したことがあった。その当時よりも利用可能なオプションが格段に増えており、使い勝手もとてもよくなっていると感じた。特に、描画を Python で行える点は、近年の流行を追っていて舌を巻く。当時は大雨を再現するだけでなく出力結果の解析にも時間と労力を費やしていたが、今後はそんな苦労も少なくなるだろうか。

    参考文献

    • wrf-model: https://github.com/wrf-model/WRF
    • wrf-python: https://github.com/NCAR/wrf-python
    • Khain, A. P., A. Pokrovsky, M. Pinsky, A. Seifert, and V. Phillips, 2004: Effects of atmospheric aerosols on deep convective clouds as seen from simulations using a spectral microphysics mixed-phase cumulus cloud model. Part I: Model description. J. Atmos. Sci., 61, 2963–2982.
    • Khain, A. P., N. Benmoshe, and A. Pokrovsky, 2008: Factors determining the impact of aerosols on surface precipitation from clouds: An attempt at classification. J. Atmos. Sci., 65, 1721–1748.
    • Khain, A. P., D. Rosenfeld, A. Pokrovsky, U. Blahak, and A. Ryzhkov, 2011: The role of CCN in precipitation and hail in a midlatitude storm as seen in simulations using a spectral (bin) microphysics model in a 2D dynamic frame. Atmos. Res., 99, 129–146.
    • Igel, A. L., and S. C. van den Heever, 2017: The role of the gamma function shape parameter in determining differences between condensation rates in bin and bulk microphysics schemes. Atmos. Chem. Phys., 17, 4599–4609.

    更新履歴

    • 2024-01-19: 初稿
    • 2024-01-20: 誤字を修正
    • 2024-01-23: 表記を修正
    • 2024-02-01: Registory の URL を修正
    • 2024-05-26: 定数ファイル読み込み時の index とカテゴリの対応表を修正
    ]]>     diff --git a/tags/Miscellaneous/index.html b/tags/Miscellaneous/index.html index 9b4c3b9..5c8c7ae 100644 --- a/tags/Miscellaneous/index.html +++ b/tags/Miscellaneous/index.html @@ -128,7 +128,7 @@ - + diff --git a/tags/Miscellaneous/page/2/index.html b/tags/Miscellaneous/page/2/index.html index 60a2e1a..35b1af5 100644 --- a/tags/Miscellaneous/page/2/index.html +++ b/tags/Miscellaneous/page/2/index.html @@ -128,7 +128,7 @@ - + diff --git a/tags/Research-paper/index.html b/tags/Research-paper/index.html index 583f188..b5b013c 100644 --- a/tags/Research-paper/index.html +++ b/tags/Research-paper/index.html @@ -128,7 +128,7 @@ - + diff --git a/tags/index.html b/tags/index.html index b5eaefd..53ed7ad 100644 --- a/tags/index.html +++ b/tags/index.html @@ -128,7 +128,7 @@ - +