フェーズドアレイレーダにおける隣接仰角間の類似性を利用した気象パラメータ推定

Weather Parameter Estimation Using the Similarity
in Adjacent Elevation Angles for Phased Array Radar
北原 大地 †
森川 侑奈 †
平林 晃 †
吉川 栄一 ‡
菊池 博史 ♯
牛尾 知雄 ∗
† 立命館大学 情報理工学部‡ 宇宙航空研究開発機構 航空技術部門♯ 電気通信大学 宇宙・電磁環境研究センター∗ 首都大学東京 航空宇宙システム工学域

Daichi KITAHARA†
Eiichi YOSHIKAWA‡

Yuna MORIKAWA†
Hiroshi KIKUCHI♯

Akira HIRABAYASHI†
Tomoo USHIO∗† College of Information Science and Engineering, Ritsumeikan University‡ Aeronautical Technology Directorate, Japan Aerospace Exploration Agency♯ Center for Space Science and Radio Engineering, The University of Electro-Communications∗ Department of Aeronautics and Astronautics, Tokyo Metropolitan University
はじめに

各仰角の散乱信号から推定する気象パラメータは局所的集中豪雨の発生頻度は最近 20 年で有意に 平均電力密度，平均ドップラー周波数，周波数幅の 3増加傾向にあり，被害も甚大なものとなっている [1]． 種類であり，それぞれ降水量，平均風速，風速の分散
この深刻な気象現象の原因は積乱雲であり，高さが と対応する．PAWR により各パラメータが時間的・
10 キロメートルを超えることもある．一方で，発生 空間的に高分解能で得られるようになったが，推定から消滅までは 10 分から 30 分程度と，短時間で状態 精度に関してはドップラー気象レーダよりも劣ってが大きく移り変わる現象である．従来のドップラー しまう．これは，ビームフォーミングで抽出された気象レーダでは，特定の約 10 から 20 の仰角方向に 散乱信号が誤差を多く含むためである．また，送信のみペンシルビームを送信し，対象仰角の気象状態 パルス数を制限することも精度劣化の原因となる．
を観測する．全方位を観測するためには，アンテナの本論文では PAWR でもドップラーレーダと同精度
向きを変えながら観測仰角数と同じ回数アンテナを の気象観測を可能にするために，ドップラーレーダに回転させる必要があり，5 分から 10 分程度の時間が おける高精度なパラメータ推定法 [8], [9] を PAWR にかかってしまう [2]．したがって，局所的集中豪雨を 応用する．文献 [8], [9] では散乱信号の離散フーリエ観測するには，時間・空間分解能が不十分であった． 係数の振幅の二乗 (離散ピリオドグラムと呼ばれる)
時間・空間分解能向上のためにフェーズドアレイ を混合密度関数でモデル化し，最尤推定法により各
気象レーダ (Phased Array Weather Radar: PAWR) が パラメータを高精度に推定した．ドップラーレーダの開発された [3]．PAWR ではフェーズドアレイ技術と 場合と違い，PAWR では隣接仰角それぞれに対して
ビームフォーミング技術を併用している．フェーズド 低精度のピリオドグラムを得ることができる．提案
アレイ技術によりアンテナの向きを変えることなく 法ではこれら隣接仰角の観測と，隣接仰角間の気象
瞬時にビームの送信方向を変えることができ，上空 パラメータの類似性を利用して，最大事後確率推定の散乱信号をまとめて一度に受信する．そして，受信 により全仰角のパラメータをまとめて求める．事後した混合信号からビームフォーミング技術 [4]–[7] に 確率の最大化は EM アルゴリズムを用いて行われる．
より各仰角の散乱信号を抽出する．これによりアン 実際の PAWR データを模倣したシミュレーションで，テナを 1 回転させるだけで全天の観測が可能となり， 非線形ビームフォーミング [7] と提案法の組み合わせ
観測時間が 10 秒から 30 秒程度と大幅に短縮された． が，非常に高精度な推定結果をもたらすことを示す．

1

2

ドップラーレーダによる気象観測

2.1 観測モデル
従来のドップラー気象レーダは，特定の仰角方向にペンシルビームを送信し，送信信号と対象物体 (雲や雨滴) の衝突によって生じた散乱信号をx
bl := xl + εl := x(lT + 2rc ) + εlと表される．最右辺の近似には，式 (2) を用いている．
ドップラーレーダによる気象観測とは，式 (1) の観測
値xbl から，式 (3) の P , µ, σ を推定することであり，それぞれ降水量，平均風速，風速の分散と対応する．
x(t) は離散的に観測されるため，離散時間信号の電力スペクトル密度を考える．xl の電力スペクトル
密度は，自己相関関数の標本値 RX (lT ) を用いて
(
)
∞
∞
∑
1 ∑
k
RX (lT )e−j2πf lT =Sd (f ) =S f−
TT(l = 0, 1, . . . , L − 1)(1)のように受信する．ここで，T [s] はパルス繰り返し
l=−∞
k=−∞
時間，r [m] は対象物体までの距離，c [m/s] は光速でκcos(2πf
T
−ν)
e
あり，x
bl ∈ C は第 l 番目の観測値，xl ∈ C は時刻 lT≈P= 2πP V (2πf T | ν, κ) (4)I0 (κ)に送信したパルスによって生じた散乱信号 x(t) の標
本値，εl ∈ C は散乱信号と独立な ∀l E[|εl |2 ] = N を と表される．ここで，1 行目から 2 行目の近似には，満たす白色雑音である．また，E[·] は期待値を表す． 式 (3) と巻き込みガウス分布がフォン・ミーゼス分布
κ cos(ω−ν)
散乱信号 x(t) は時間的にゆらいでいる不規則信号 V (ω | ν, κ) := e 2πI0 (κ) で近似できることを用いた．
(確率過程 X(t) の 1 つの見本過程) であり，xl の値も ν := 2πµT [rad] は正規化された平均ドップラー角周
確率変数の実現値となる．一般に，散乱信号 x(t) は弱 波数，κ > 0 は 1/σ 2 に相当する変数，I0 (κ) は 0 次
定常1 かつエルゴード的2 であることが知られている． の第一種変形ベッセル関数である．電力スペクトル
更に，散乱信号の平均値は E[X(t)] = 0 で，自己相関 密度 Sd (f ) の推定には，離散ピリオドグラム
2
関数 RX (τ ) = E[X(t + τ )X(t)] の振幅は τ が大きく
L−1
∫∞
1 ∑
−j 2πklLxl e(k = 0, 1, . . . , L − 1) (5)なると急激に減少するため， −∞ |τ RX (τ )| dτ < ∞ hk :=L
l=0が成り立つ．厳密には，x(t) はパルス幅とビーム幅で決まる散乱体積内の独立な散乱信号の総和であり， を用いる．実際に，離散ピリオドグラムの期待値は[( L−1)( L−1
)]
∑
∑
∑
2πkl
2πkl
1
x(t) =as (t)ej(2πfs t+ϕs ) ≈ a(t)ej(2πµt+ϕ) (2)−j L 1j L2
E[hk ] = E
xl 1 ex̄l2 e
s
L
l1 =0
l2 =0
と表される [10]–[13]．ここで，as (t) ≥ 0 は第 s 番目の
L−1
L−1
2πk(l1 −l2 )1 ∑∑散乱信号の振幅，fs [Hz] はドップラー周波数，ϕs [rad]L
E[xl1 x̄l2 ]e−j=L
l1 =0 l2 =0
は見本過程ごとに [0, 2π] の一様分布から生成される
2 σ2 τ 2−2π
L−1 L−1
初期位相であり，a(t) ≥ 0 は RA (τ ) = P e
を() −j 2πk(l1 −l2 )
1 ∑∑
L
=
R
(l
−
l
)T
e
1
2
X
満たす不規則振幅，µ [Hz] は平均ドップラー周波数，L
l1 =0 l2 =0
ϕ [rad] はランダムな初期位相である．最右辺の近似
L−1
L−1
∑
∑ |l|RX (lT )
2πkl
2πkl
では，形状の等しい散乱体がランダムに分布し，fs=
RX (lT )e−j L −e−j LLの確率密度関数が平均 µ のガウス分布であると仮定
l=−(L−1)
l=−(L−1)(6)
している．また，j ∈ C は虚数単位，P := RX (0) は
∑∞
|lR
(lT
)|
<
∞
であることから，となり，散乱信号の平均電力密度，σ [Hz] は周波数幅である．
X
l=−∞
(
)
∞
x(t) の電力スペクトル密度は，自己相関関数を用いて∑k−j 2πkl∫ ∞L
lim E[hk ] =RX (lT )e= Sd(7)(f −µ)2
PL→∞LT
RX (τ )e−j2πf τ dτ ≈ √S(f ) =
e− 2σ2l=−∞2πσ−∞(3) が成り立つ．また，任意の L に対して，∫[ L−1 ]
[ L−1]∞
1
確率過程 X(t) が，∀t E[X(t)] = −∞ xp(x | t) dx = µX かつ
∫∫ ∞
1 ∑
1 ∑
2
∀t ∀τ E[X(t+τ )X(t)] = −∞ x1 x̄2 p(x1 , x2 | t+τ, t) dx1 dx2 =E
hk = E|xl | = P
(8)LL
RX (τ ) である (期待値が t に依らない) とき，X(t) は弱定常であるという．ここで，p(x | t) は x(t) の確率密度関数，p(x1 , x2 | t+τ, t)
は x1 = x(t + τ ) と x2 = x(t) の同時確率密度関数である．
2
弱定常過程 X(t) がエルゴード的であるならば，任意の見本∫T
過程 x(t) に対して，µX = limT →∞ T1 0 x(t) dt かつ RX (τ ) =∫TlimT →∞ T1 0 x(t + τ )x(t) dt が成り立つ．ここで，確率変数の極限は二乗平均の意味で収束することに注意されたい．

k=0

l=0

k
も成り立つ．本論文では，式 (6) の右辺が Sd ( LT)に似ると考え，離散ピリオドグラムのモデルを()
k
V (ωk | ν, κ)E[hk ] = CL Sd= P L ∑L−1
(9)LT
k′ =0 V (ωk′ | ν, κ)と簡略化する．ここで，CL :=2π∑L−1 L，V (ωk′ | ν,κ)k′ =0

ωk := 2πk
L であり，limL→∞ CL = 1 となる．式 (9) は，
式 (4), (7), (8) とは矛盾しないことに注意されたい．

2
1.8
1.6
1.4

2.2 EM アルゴリズムによるパラメータ推定1.2

式 (5) の離散ピリオドグラムは雑音 εl を含まない
理想的な値であり，実際は，εl を含む観測値 x
bl から

1
b
hk :=
L

L−1
∑

x
bl e

0.8
0.6

2
−j 2πklL1(k = 0, 1, . . . , L − 1) (10)

0.4
0.2

l=0

のように離散ピリオドグラムを計算する．xl と εl は互いに独立で，εl は白色雑音であることから

E[ bhk ] = E[hk ] + N = 2πCL P V (ωk | ν, κ) + Nが成り立つ．したがって，以下の混合密度モデル

E[ bhk ] ∝ ρV (ωk | ν, κ) + (1 − ρ)U (ωk ) =: g(ωk | Θ)(11)が得られる (図 1)．ここで，ρ ∈ (0, 1) は
ρ :=

CL P
P
≈
CL P + N
P +N

0

-3

-2

-1

0

1

2

3

図 1: 離散ピリオドグラム b
hk と式 (9) の混合密度関数
は，ベイズの定理を用いて
ρ(i) V (ωk | ν (i) , κ(i) )

(i)


αk =
g(ωk | Θ(i) )(1 − ρ(i) )U (ωk )(i) βk =g(ωk | Θ(i) )

(14)

(12)

と表される．これらの値は，ωk が気象エコーと雑音と定義され，最右辺の近似は L の値が大きいときに それぞれの分布に属する確率を表しており，負担率1
は [0, 2π] 上の一様分布， とも呼ばれる．完全データ対数尤度の期待値は成立する．また，U (ω) := 2π
L−1Θ := (ρ, ν, κ) は混合密度関数 g のパラメータを表す．
∑ [ (i)
(
)
(i)
b
α
Q
(Θ)
=
h
log
ρV
(ω
|
ν,
κ)
k
k
k
式 (12) から，ρ は観測信号 x
bl の電力 P + N のうち，
k=0
気象エコー xl の電力 P が占める割合を表している．
(
)]
(i) b
+
β
h
log
(1
−
ρ)U
(ω
)
(15)
k
k
k
文献 [8], [9] では，式 (10) の離散ピリオドグラム
から，EM アルゴリズムを利用して混合密度関数 gのパラメータ Θ を推定している3 ．式 (11) の関数 gはあくまで E[ bhk ] のモデルにすぎないが，bhk の値を正規化ドップラー角周波数 wk の観測回数と見なし，
g は ω ∈ [0, 2π] の確率密度関数であると近似する．
これによって，混合確率分布のパラメータ推定問題と同様に，EM アルゴリズムで Θ が推定可能となる．
EM アルゴリズムは，対数尤度関数の最大化問題

maximize
Θ

L−1
∑

(
)
b
hk log g(ωk | Θ)

(13)

k=0

を解く手法であり，ある初期値 Θ(1) から E ステップと M ステップの 2 つを繰り返すことで解を求める．
E ステップでは，現在の推定値 Θ(i) のもとで得られる潜在変数の事後確率に関して，完全データ対数
尤度関数の期待値を計算する．潜在変数の事後確率のように計算される．
M ステップでは，式 (13) の対数尤度関数の下界でもある式 (15) の期待値 Q(i) (Θ) を最大化する．具体的には，Q(i) (Θ) を各変数 ρ, ν, κ で偏微分し，各変数に関して最大値を求めることで Θ を更新する．まず，∂Q(i)∂ρ = 0 から

∑L−1

ρ

(i+1)

(i) b
k=0 αk hk(i) b
L−1 (i)k=0 (αk + βk ) hk

=∑

(16)

(i)

と更新する．次に， ∂Q∂ν = 0 から∑L−1 (i) b
αk hk sin ωk
(i+1)
tan ν
= ∑ k=0 (i)L−1α bhk cos ωk
k=0

(17)

k

が成り立つ．この式を満たす ν ∈ [0, 2π) は 2 つあり，
一方は式 (15) を最小に，もう一方は最大にするため，(i)
最大値の方を ν (i+1) とする．最後に， ∂Q∂κ = 0 から3
文献 [8], [9] では，散乱信号には，山岳，森林，建築物など∑L−1 (i) bからの信号であるグランドクラッタも含まれているとしていた．
(i+1) )
I1 (κ(i+1) )k=0 αk hk cos(ωk − ν
本論文では，グランドクラッタが存在しない式 (11) のモデルに=(18)∑
L−1 (i) b
I0 (κ(i+1) )対して，文献 [8], [9] のアルゴリズムを適用している．
α hk
k=0

k

が成り立つ．この方程式の解は陽に求まらないため， 3.2 ビームフォーミング
ニュートン法やグリッドサーチなどで解を計算する． 3.2.1 線形ビームフォーミング
式 (14)–(18) を繰り返すことで，式 (13) の問題を解く．
ビームフォーミングの中でも，最も基本的な手法で3 フェーズドアレイレーダによる気象観測
ある Fourier 法 (以下 FR 法) [4]，サイドローブの寄与を低減する Capon 法 (以下 CP 法) [5]，少数の観測値
3.1 観測モデル
からでもロバストに推定できる MMSE 法 [6] の 3 つ
従来のドップラーレーダでは，10 方向程度の仰角を説明する．これらの手法はいずれも，観測信号 yl
に対して，散乱信号をそれぞれ L = 80 回程度観測と複素重みベクトル wm ∈ CK の複素内積の計算することで，気象パラメータ P , µ, σ を高精度に推定HHTしていた．しかしながら，レーダを中心とする 3 次元
b l := (wHx(20)
1 y l , w2 y l , . . . , wM y l ) = W y lの上半球のスキャンに 5 分から 10 分の時間がかかり，により xl の推定を行う，y l に関して線形なビーム
仰角標本数も少ないため，近年増加しているゲリラ
フォーミングである．FR 法の複素重みベクトルは豪雨等の局所的大雨を観測しきれない．時間的・空間
sm
wFR,m :=
(21)
的に高分解能な観測を実現するために，フェーズドK
アレイ気象レーダが開発された [3]．このレーダは，と表される．FR 法は高速ではあるが，サイドローブ
観測対象の仰角区間全体 [θmin , θmax ] にファンビーム
で他方向の信号も同時に受信してしまう欠点がある．
を送信し，K 個あるアンテナ素子で散乱信号をそれCP 法では，対象仰角にメインローブを向けながらぞれ yn,l ∈ C (n = 1, 2, . . . , K) と観測する．全仰角
受信電力を最小化する複素重みベクトル∆θ
区間を M 個の部分区間 [θm − ∆θ
2 , θm + 2 ] に分割
b −1 smRyし，それぞれの部分区間における散乱信号を xm,l ∈ C
wCP,m :=
(22)H
θmax −θminb y−1 smsR(m = 1, 2, . . . , M ) とする．ここで，∆θ :=m
M
かつ θm := θmin + (m − 21 )∆θ である．アンテナ素子 を適応的に求める．ここで，Rb y := 1 ∑L−1 yl y H で
l=0l
Lは直線状に間隔 d [m] で配置されており，送信信号の ある．CP 法は，観測数 L が十分大きいとき，他方向波長を λ [m]，各素子の観測雑音を εn,l ∈ C とすれば， の信号を良好に抑圧できる．しかし，L が小さくなるM
b y の逆行列も定義∑
につれて電力が過少に推定され，Ryl =xm,l sm +εl = Sxl +εl (l = 0, 1, . . . , L−1)
できなくなる．MMSE 法では，複素重みベクトルを
m=1

( L−1
)
(19)

∑
1

(i)
(i)
(i)H
b MMSE,l x
b MMSE,l ⊙ IM
Rx =x
のように観測信号 y l := (y1,l , y2,l , . . . , yK,l )T ∈ CK

L

l=0
を表現できる．ここで，ステアリングベクトル sm は

(i)
(i)

(
2(K−1)πd sin θm )T
2πd sin θm
Ry = SRx S H + N IKλλsm := 1, e−j, . . . , e−j
∈ CK

(i)−1

Ry sm
 (i+1)T
Mで表され，xl := (x1,l , x2,l , . . . , xM,l ) ∈ C ，εl :=(m = 1, 2, . . . , M )wMMSE,m =
(i)−1

sH
sm
m Ry
(ε1,l , ε2,l , . . . , εK,l )T ∈ CK ，S := (s1 , s2 , . . . , sM ) ∈

(i+1)
(i+1)
b MMSE,l = WMMSE y l (l = 0, 1, . . . , L − 1)x
CK×M である．また，Rε := E[εl εHl ] = N IK とする．
(23)
式 (19) における xm,l は，式 (1) における xl に相当する．したがって，フェーズドアレイレーダによる のように求める．ここで，⊙ はアダマール積，IM は
b y を計算
気象観測では，最初にビームフォーミング (詳細は第 M ×M の単位行列である．CP 法と異なり R
3.2 節) を用いて y l から xl を推定する．その後，各 しないため，小さな L でもロバストに推定を行える．
仰角方向の平均電力密度 Pm ，平均ドップラー周波数
fm ，周波数幅 σm を推定することで，時間的・空間的に高分解能な気象観測が可能となる．一方で，得られる気象パラメータ Pm , fm , σm の精度に関しては，
ドップラーレーダの方が優れている．これは，ビーム
フォーミングによる推定結果 x
bm,l が誤差を多く含むためである．また，時間分解能向上のために L = 20回程度しか観測しないことも精度劣化の原因となる．

3.2.2

非線形ビームフォーミング

筆者らは仰角分解能向上のために，凸最適化問題

1
T G2
1
∥Y −SX∥2F + ξ1 ∥X∥G
1 + ξ2 ∥BF X ∥1
2
(24)
K×L
を解くことで Y := (y0 , y1 , . . . , yL−1 ) ∈ Cから
X := (x0 , x1 , . . . , xL−1 ) ∈ CM ×L を推定するビーム
フォーミングを提案した [7]．ここで，ξ1 , ξ2 > 0 は各
minimizeX項の重み，F ∈ CL×L はフーリエ変換，B ∈ RbL×L
はサイズ b のブロック抽出行列，∥·∥G
1 はグループ ℓ1c
ノルムである．推定結果 X は ADMM [14] を用いて計算され，式 (20) では表せない非線形な手法 (以下NL 法) となり，散乱信号 xm,l を高精度に推定できる．

4

隣接仰角間の類似性を利用したパラメータ推定

ビームフォーミングによって得られた散乱信号の推定値 xbm,l (l = 0, 1, . . . , L − 1) から，式 (10) と同様に離散ピリオドグラム b
hm,k (k = 0, 1, . . . , L − 1) を計算できる．したがって，仰角毎に第 2.2 節の手法を適用すれば，各仰角の混合密度関数 g のパラメータΘm := (ρm , νm , κm ) が推定される．しかしながら，
ビームフォーミング結果 x
bm,l に含まれる残存誤差と観測回数 L の少なさに起因して，パラメータの推定
精度はドップラーレーダと比べて低下してしまう．
本論文では，フェーズドアレイレーダでドップラー
レーダと同精度の気象観測を行うために，全仰角のパラメータ Θ := (Θ1 , Θ2 , . . . , ΘM ) をまとめて推定する手法を提案する．ドップラーレーダでの観測と異なり，フェーズドアレイレーダでは連続する仰角区間∆θ[θm − ∆θ2 , θm + 2 ] それぞれに対してピリオドグラム
bhm,k を取得できる．隣接仰角間では，平均電力密度
Pm ，平均ドップラー周波数 µm ，周波数幅 σm は類似すると考えられ，この性質を Θ の事前分布によって表現する．提案法は，最大事後確率 (MAP) 推定

maximize
Θ

M L−1
∑
∑

(i)

(i)

M L−1
∑
∑ [ (i)()
Q (Θ) =
αm,k bhm,k log ρm V (ωk | νm , κm )(i)

m=1 k=0()](i) b+ βm,k hm,k log (1 − ρm )U (ωk ) − Ψ(Θ)

(27)

(i)

(i)
(i)
2 (i+1)3
2ζ1 Pbmρm− ζ1 Pbm (2Pbm + Pem+1 + Pem−1 )ρ(i+1)2
m
[
]
L−1
∑ (i)
(i)
(i)
(i)
+ Pe
)− (α +β )b+ ζ1 Pbm (Pehm,k ρ(i+1)
m+1

m−1

m,k

m,k

m

k=0

+

L−1
∑

(i)
hm,k = 0
αm,k b

(28)

k=0
(i)
(i)
を満たす．ただし，Pem±1 := ρm±1 Pbm±1 である．これは ρm ∈ (0, 1) に関する 3 次方程式であるので，解は最大 3 つ存在する．その中で式 (27) の関数値を最も(i)大きくする解に更新する．次に， ∂Q∂νm = 0 から
(i+1)
tan νm
∑L−1

(i) b
(i)
(i)
k=0 αm,k hm,k sin ωk − ζ2 (sin νm+1 + sin νm−1 )
=∑
(i)
(i)
L−1 (i) b
k=0 αm,k hm,k cos ωk − ζ2 (cos νm+1 + cos νm−1 )(29)が成り立つ．この式を満たす νm ∈ [0, 2π) は 2 つあり，
一方は式 (27) を最小に，もう一方は最大にするため，
(i)
(i+1)
最大値の方を νmとする．最後に， ∂Q∂κm = 0 から
(i+1)

(i+1)

(i)

(i)

ζ3 (2κm − κm+1 − κm−1 )+∑L−1 (i) b(i+1)
I0 (κm )k=0 αm,k hm,k∑L−1 (i) b(i+1))
k=0 αm,k hm,k cos(ωk − νm=
∑L−1 (i) bk=0 αm,k hm,k
I1 (κm

m=1 k=0と定義される．ただし，ζ1 , ζ2 , ζ3 > 0 は各項の重み，∑
Pbm := L1 L−1xm,l |2 は雑音を含む電力 Pm + N の
l=0 |b推定値，ρm Pbm は式 (12) から Pm の推定値である．
式 (25) の最大化問題を EM アルゴリズムによって(i)
解く．E ステップでは，現在の推定値 Θm のもとで完全データ対数尤度関数の期待値を計算し，下界関数(i)ν (i) , κ(i) を ρm , νm , κm に置き換えて得られる負担率である．M ステップでは，Q(i) (Θ) を最大にする
(i)
(i+1)
は
パラメータを求める．まず ∂Q∂ρm = 0 から，ρm()
bhm,k log g(ωk | Θm ) − Ψ(Θ)(25)
により Θ を高精度かつロバストに計算する．ここで，Ψ は Θ の対数事前分布を −1 倍した関数に相当し，
M
−1[
∑
ζ1
Ψ(Θ) :=(ρm+1 Pbm+1 − ρm Pbm )2
2
m=1
]
ζ3
2
(26)
− ζ2 cos(νm+1 − νm ) + (κm+1 − κm )2(i)

を設計する．ここで，αm,k と βm,k は式 (14) の ρ(i) ,

)

(30)

が成り立つ．この方程式の解は陽に求まらないため，
ニュートン法やグリッドサーチなどで解を計算する．
更新する仰角 θm をランダムに選びながら，式 (28)–(30) を繰り返すことで，式 (25) の問題を解く．また，
ζ1 = ζ2 = ζ3 = 0 とすれば，式 (28), (29), (30) はそれ
ぞれ，従来法の更新式 (16), (17), (18) と一致する．

5

数値シミュレーション

実データを模倣した数値シミュレーションにより提案法の有効性を示す．大阪大学吹田キャンパスに設置された PAWR で 2014 年 3 月 30 日に観測されたデータを用いた．このデータの第 30 番目の方位角の
距離 r = 7.5 [km] 地点において，θmin = −15 [deg]から θmax = 30 [deg] の 55 個のレーダ反射因子を 3次スプライン補間して小さな雑音を加えることで，
M = 110 個の真の平均信号電力密度 Pm を作成した．
真の νm と κm は正弦波に定数と小さな雑音を加える

ことで作成した．散乱信号 xm,l は周波数領域で作成され，式 (9) を満足する複素ガウス分布から各離散フーリエ係数を独立に生成した．PAWR の各設定は，アンテナ素子数 K = 128，パルス数 L = 20，波長λ = 0.0318 [m]，素子間 d = 0.0165 [m]，パルス繰り
返し時間 T = 0.0004 [s] とした．雑音電力 N は信号
電力 Pm の最小値の半分とした．ビームフォーミング
には，L < K で CP 法 (式 (22)) は使用不可能なため，
FR 法 (式 (21))，MMSE 法 (式 (23))，NL 法 (式 (24))の 3 種類を用いた．NL
法の各パラメータは，q = 5，√
K qL
b = 3，ξ1 = 0.0025 M ，ξ2 = 0.25 MK√b とした [7]．
ドップラーレーダで L = 20 回ペンシルビームを送信した後に従来法を用いた結果，PAWR でビーム
フォーミングした後に従来法を用いた結果，PAWR
でビームフォーミングした後に提案法を用いた結果を比較する．提案法では式 (26) の重みを，ζ1 = 2.0 ×−2 と設定した．
10−10 ，ζ2 = 3.5×10−3 ，ζ3 = 4.5×10∑[4] E. Kudeki and F. Sürücü, “Radar interferometric imaging of field-aligned plasma irregularities in the equatorial electrojet,” Geophys. Res. Lett., vol. 18, no. 1,
pp. 41–44, 1991.

∑
∑
100 M
κm −κm100 M
νm −νm )2
m=1 (bm=1 (b
∑
∑
，
M
M
2
2
ν
m=1 m
m=1 κm[5] J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” Proc. IEEE, vol. 57, no. 8, pp. 1408–
1418, 1969.

表 1 は，各推定値の正規化誤差

100

M
(b
ρm Pbm −Pm )2
m=1∑M
2
m=1 Pm)2，を表している．

パルス数が少ないためドップラーレーダの場合でも高精度な推定は難しく，特に周波数幅の逆数に相当
する κm の推定精度は極めて低い．PAWR でビーム
フォーミング後に従来法を用いた場合には，各仰角の離散ピリオドグラムの精度が下がるため，Pm , νm ,
κm の推定精度も下がる．一方，ビームフォーミング
後に提案法を用いた場合には，νm と κm の推定精度が大きく向上し，特にビームフォーミングに NL 法を用いたときが最も高精度にパラメータを推定できた．

6

おわりに

表 1: パルス数 L = 20 における正規化誤差
Pm
νm
κm
ペンシルビーム

FR + 従来法
FR + 提案法
MMSE + 従来法
MMSE + 提案法
NL + 従来法
NL + 提案法

43.18
62.53
64.04
47.69
47.55
44.35
43.58

19.02
33.86
30.63
17.80
15.85
20.42
12.21

1249.501565.5444.92
1363.2039.50
1364.05
24.00[3] F. Mizutani, T. Ushio, E. Yoshikawa, S. Shimamura,
H. Kikuchi, M. Wada, S. Satoh, and T. Iguchi, “Fastscanning phased-array weather radar with angular imaging technique,” IEEE Trans. Geosci. Remote Sens.,
vol. 56, no. 5, pp. 2664–2673, May 2018.

[6] E. Yoshikawa, T. Ushio, Z. Kawasaki, S. Yoshida,
T. Morimoto, F. Mizutani, and M. Wada, “MMSE beam
forming on fast-scanning phased array weather radar,”IEEE Trans. Geosci. Remote Sens., vol. 51, no. 5,
pp. 3077–3088, 2013.
[7] D.Kitahara, M.Nakahara, A.Hirabayashi, E.Yoshikawa,
H. Kikuchi, and T. Ushio, “Nonlinear beamforming via
convex optimization for phased array weather radar,”in Proc. Asia-Pac. Signal Inf. Process. Assoc. Annual
Summit Conf. (APSIPA ASC), 2018, pp. 1831–1835.
[8] S. Kon, T. Tanaka, H. Mizutani, and M. Wada, “A machine learning based approach to weather parameter estimation in Doppler weather radar,” in Proc. IEEE Int.
Conf. Acoust. Speech Signal Process. (ICASSP), 2011,
pp. 2152–2155.

本論文では，フェーズドアレイレーダによる気象観測に対して，隣接仰角間の類似性に基づく高精度なパラメータ推定法を提案した．ビームフォーミング
[9] 根 智志, 田中 聡久, “ロバストな混合モデルによるドップラー気象レーダの気象データ推定,” 電子情報通信
で得られる離散ピリオドグラムは低精度であるため，
学会技術研究報告, 2011, vol. 111, no. 104, pp. 79-84.
従来法をそのまま適用すると気象パラメータの推定[10] V. N. Bringi and V. Chandrasekar, Polarimetric Doppler
精度も劣化してしまう．提案法では，隣接仰角間のWeather Radar: Principles and Applications. Cam類似性を事前分布で表現し MAP 推定を行うことで，
bridge University Press, 2001.
少ないパルスでも非常に高精度な推定を実現した． [11] 関根 松夫, レーダ信号処理技術. 電子情報通信学会,
参考文献[1] Committee on Extreme Weather Events and Climate
Change Attribution, Attribution of Extreme Weather
Events in the Context of Climate Change. Washington D.C., The National Academy Press, 2016.
[2] M. Wada, J. Horikomi, and F. Mizutani, “Developmentof solid-state weather radar,” in Proc. IEEE Radar Conf.,
2009, 4 pages.

1991.
[12] 吉野 文雄, レーダ水文学. 森北出版, 2002.
[13] 深尾 昌一郎, 浜津 享助, 気象と大気のレーダーリモー
トセンシング. 京都大学学術出版会, 2005.
[14] D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear variational problems via finite element approximation,” Comput. Math. Appl., vol. 2,
no. 1, pp. 17–40, 1976.