6.4 3次元変分法（JNoVA0）の開発 1
6.4.1 はじめに初期値は、数値予報の精度を左右する重要な要素の一つである。現在主流となってきている変分法データ同化は、観測の情報をより活かすことのできる、優れた初期値解析法であり、気象庁の全球解析やメソ解析で現業的に利用されている。一方、これら既存の変分法データ同化システムは、第6.1節で既に述べられたように、静力学モデルを対象としているため、静力学平衡の仮定、微物理量の解析が十分でないなどNHMの初期値解析法としては不十分な点がある。また、変分法データ同化では、評価関数に含まれる行列の次元が大きすぎるため、いくつか仮定を置いて計算量を減らす必要があるが、この際、モデルの誤差特性などを考慮した適切な設計が必要である。以上のことから、NHMの予報の高精度化に
は、専用に設計された高度な変分法データ同化シス
テムが不可欠であり、平成14年4月より、NHM用の変分法データ同化システム（JNoVA）の開発に本格的に着手した。
これから半年以上が経過した現在、3次元変分法のプロトタイプ（JNoVA0）の動作を確認しており、今後改良を加えることとしている。本節では、現在開発中のJNoVA0の設計及び今後の展望について述べる。変分法の理論及び現業化に関しては、露木（1997）及び気象庁予報部（2002）などを参照して頂きたい。なお、JNoVA0は開発の初期段階にあり、本節の記述内容は今後の開発により変更される可能性があることに留意願いたい。
6.4.2 JNoVA0の設計方針
JNoVA0の設計の骨格には、現在現業運用中のMSM用4次元変分法（メソ4DVAR）と同様の方法（石川 2002）を用いている。すなわち、背景誤差共分
散行列の0でない成分（以下「計算成分 2 」）を減らし、最適解探索における収束を加速するため、状態変数x 3 を以下のように変数変換する。
①互いに独立とみなす変数に変数変換②鉛直共分散を独立にするため固有モードに変換③水平相関行列をCholesky分解した行列で変換
この際、いくつか仮定を置くことで、背景誤差共分
散行列の計算成分を減らしている。ここで重要なのは、解析にとって重要な情報を減らさずに、いかに計算量を減らすか、ということである。JNoVA0で1

三好建正

は、この①から③までの変数変換を行った変数を制御変数として用いている。
この変数変換により、JNoVA0の評価関数がどのように変形されるか見てみよう。
まず、3次元変分法の評価関数 J は、次式で与えられる。
J = 12 (x − x b ) T B −1 (x − x b )+ 12 ( Hx − y o ) T R −1 ( Hx − y o )ここで、 x は解析しようとする状態変数、 x b は状態変数の第一推定値、B は背景誤差共分散行列、H は観測演算子（一般には非線形）、 y o は観測値からなるベクトル、 R は観測誤差共分散行列である。状態
変数のインクリメント δx = x − x b を定義し、(6.4.1)
式の評価関数をインクリメント法で書き直すと、第1次近似で次式のようになる。
J = 12 δx T B −1δx + 12 (Hδx − d) T R −1 (Hδx − d)(6.4.2)

ここで用いた第1次近似は、次式による。
Hx = H (x b + δx) = H (x b ) + ∂∂Hx δx + O(| δx | 2 )≈ H (x b ) + Hδx(6.4.3)

ここで、線形化された観測演算子 H = ∂∂Hx 及びD値
(departure) 4 d = y o − Hx b を用いた。
次に、①及び②の変数変換 C （ C : δx a δu ）により、 B u = δu(δu) T の変数間共分散及び鉛直共分散の非対角成分が0となる。なお、この変数変換は一
般には非線形だが、インクリメント法では δx を微少な摂動とみなすことにより線形化できる。線形化した変数変換を C と書き、 Cδx = δu とおいて(6.4.2)式
を変形すると、次式を得る。
J = 12 δu T B u−1δu + 12 (HC −1δu − d) T R −1 (HC −1δu − d)(6.4.4)

最後に、これまでの手続きでは対角化されずに残っている水平共分散をCholesky分解する③の手続
き（ B u = LLT ）により、評価関数は次式のように変形される。
J = 12 (L−1δu) T (L−1δu)+ 12 (HC −1 L(L−1δu) − d) T R −1 (HC −1 L(L−1δu) − d)(6.4.5)2 0 となっている部分の計算をスキップすることにより、実際上の計算量を減らすことができる。0 でない成分のみを計算することから、ここでは計算成分と呼んでいる。
3 全格子点の全予報変数からなるベクトルをいう。

(6.4.1)4 (観測値)－(第一推定値)

ここで、 δv = L−1δu とおけば、この δv が制御変数となり、評価関数は最終的に次式のようになる。
J = 12 δv 2 + 12 (HC −1 Lδv − d) T R −1 (HC −1 Lδv − d)(6.4.6)(6.4.6)式には B が陽には現れておらず、①～③の変数変換の手続きで表現されている。①及び②の変数変換は、背景誤差共分散行列をブロック対角化 5 している。実際は、①の変換後の変数間相関は完全には0にならないが、0でない成分を無視し、互いに独立とみなすことにより、計算成分を減らしている。
また、①～③のすべての手続きで水平一様性を仮定
し、計算成分をさらに減らしている。③の手続きは、
右辺第1項から行列変換を除いており、準ニュートン法による最適解探索の際の収束を加速させる効果がある 6 。共分散行列は対称行列なので、正定値であれば完全にCholesky分解可能だが、ここでは水平誤差相関の広がりをGauss分布で近似し、正定値性を確保するとともに、情報を東西及び南北方向の半値距離のみとしている。こうしても計算量は減少しな
いが、必要とされる統計情報が激減する。背景誤差
共分散行列の計算成分を減らすための仮定についてまとめると、以下の通りである。
・ 水平一様性の仮定
・ ①の変数変換で変数間相関を無視
これらの仮定により B の計算成分が激減する。近似を行わずに背景誤差共分散行列を対角化することも可能であり、最適解探索の際の収束は加速する。
しかし、この場合、対角化するための変数変換の方に情報量がすべてシフトし、事実上計算量を減らす
ことにはならない点に注意が必要である。
②及び③の変数変換の手続きは機械的に行われるが、①の変数変換については任意性がある。背景誤
差共分散のうち変数間の共分散は、①の変数変換で表現され、解析結果に大きな影響を与えるため、モデルの特性を考慮した慎重な設計が求められる。このため、実際の予報実験によりNHMの背景誤差共分散行列を推定し、どのような変数を選べば互いの相関を無視できるか調査しつつ変数変換を設計する必要がある。

5

水平誤差共分散を除く成分を対角化している。
通常、観測ベクトルの次元は状態変数の次元よりかなり小さく、評価関数の第 1 項（背景項）が第 2 項（観測項）に比べて大きな割合を占める。このため、背景項を対角化すれば、評価関数の形状が球形に近くなり、最適解探索の際の収束が加速する。正確に言えば、Hessianの条件数が1 に近づく（露木(1997)のp.109 を参照）。
6

6.4.3 背景誤差共分散行列の統計的推定
背景誤差共分散行列の推定には、NMC法（Parrishand Derber 1992）を用いた。NMC法は、初期時刻の異なる2つの予報について、同じ予報対象時刻（valid time）での予報値の差を予報誤差と見なす方法であり、通常は24又は12時間程度を初期時刻の差とする。今回は、局所域モデル（600km×600km、
水平解像度は5km）を用いたため、通常のNMC法のように24時間や12時間の差では境界から流入する
外部モデルの影響を受けて有意な情報が得られない。
このため、6時間の差とした。具体的には、1日2回の同じvalid timeについて12時間予報値と6時間予報値の差をとることにより、予報誤差とした。1日2回としたのは、昼と夜の差が統計に現れないようにするためである。
NMC法は、状態変数のすべての要素について誤差が推定できるという点で優れており、2メンバーのアンサンブル予報による誤差の推定と見なせる。多数のメンバーを用いたアンサンブル予報による誤差推定法（解析アンサンブル法）もあるが、NMC法は簡便であるため、よく使われている 7 。
今回は、2002年4月24日から2002年5月30日の間の31日間62例について、以上の方法により予報誤差を推定し、その共分散行列を求めることにより、背景誤差共分散行列を統計的に推定した。なお、鉛直層数は45層とした。
6.4.4 変数間誤差相関を表現する変数変換
モデルは現実に近い気象現象を表現するため、気象学的に妥当なバランスが予報誤差統計に現れ、これにより変数間誤差相関が説明されると考えることは、変数間相互のバランスを考える上でごく自然な発想である。
具体的には、地衡風の関係等から推測される水平
風と気圧の相関、及び、質量保存の関係から推測される水平風と鉛直風の相関について考える。前者は静力学モデルでも考慮されるものであるが、後者については鉛直風を予報変数として持つ非静力学モデ
ル特有のものである。
(1) 気圧と水平風の相関
まず、気圧と水平風の相関が、大気の運動方程式の定常解として表される力学バランスにより表現されることを考える。
大気の運動方程式の水平成分は、一般鉛直座標系
ζ において次式のように与えられる（Haltiner andWilliams 1980; Kalnay 2002）。

7 筆者の知る限り、ECMWFを除く各数値予報センターではNMC法が用いられている。なお、Fisher (2001)によると、ECMWFは解析アンサンブル法を用いている。

1∂v & ∂v+ζ= − ∇ z p − v ⋅ ∇ ζ v − fk × v + F∂t∂ζρ(6.4.7)

ここで、この方程式は2次元で書かれており、 ∇ は2次元の微分演算子、v は ζ 面上での水平風ベクトル、
ρ は密度、 p は気圧、 f はコリオリパラメータ、 k
は鉛直方向の単位ベクトル、 F は摩擦などの外力である。これ以降本節では、z面上での微分を表す場合
のみ微分記号に添字を記し、 ζ 面上での微分を表す
場合は省略する。(6.4.7)式の定常解は、左辺の時間微分項を0とおくことにより得られる。左辺第2項（鉛直移流項）は他の項に比べて十分小さいため無視し、右辺第4項（外力項）を除くと、定常解の満たすべき方程式として次式を得る。
1

ρ

∇ z p B = − v ⋅ ∇ v − fk × v(6.4.8)

この式における気圧は、水平風とバランスしているものであることを明示するため、添字 B をつけた。
ここで、次式で定義される質量変数 Φ を導入する。
∇Φ =

1

ρ

∇z p

(6.4.9)

これを用いて、(6.4.8)式は次式のようになる。
∇ Φ B = − v ⋅ ∇ v − fk × v(6.4.10)

厳密にはこの方程式を扱う必要があるが、これを扱うにはPoisson方程式を解く必要があり、計算コストが高い。このため、次のように簡単化する。総観
規模のバランスを考える場合、移流項はコリオリ項に 比 べ て 十 分 小 さ い た め 無 視 し 、 (6.4.8) 式 を
v = −∇ ×ψk で定義される流線関数 ψ を用いて変形すると、次式を得る。
1

ρ

∇ z p B = f∇ψ

デル最下層の気圧 p surf と各層の温位 θ を用いて静力学平衡の関係により計算される気圧（以下「バランス気圧」）だと考える。
以上の議論から、気圧とのバランスを差し引いた
流線関数の非バランス成分（非バランス流線関数
ψ U ）を変換後の独立変数と選ぶことにより、水平
風と気圧の地衡風関係が表現される。なお、気圧は
R/c
Exner関数 π = pp
を用いる。

( )

ψ U = ψ (u , v) − r1 (π B ( p surf , θ ) )c / R
p(6.4.13)

ここで、r1 は第6.4.3項の統計から決まる回帰係数であり、考慮したバランスが統計的にどの程度成り立
つかを示すような係数である。地衡風関係は(6.4.12)式により表現されるため、線形回帰関係で説明がつc /R= 0 の条件から、次式のく。正確には、 ψ U , π Bp
ように決まる。なお、括弧<,>は共分散を表す。

ψ ,π Bcp / Rr1 =πB

cp / R,π B

cp / R

(6.4.14)以上のように、水平風2成分のうち1成分は流線関数にした。水平風のもう一つの成分として、v = −∇χ
で定義される速度ポテンシャル χ を考える。流線関数と速度ポテンシャルの水平相関を調べたところ、
図6.4.1のようにおおむね同心円状となったため、線形回帰関係でほぼ説明がつくものと考えられる 8 。この相関を考慮することにより、地衡風以外の外力の効果等による影響がある程度説明される。特に、地表面摩擦などの影響が大きい下層においては、(6.4.11)

風の鉛直シアがそれほど強くない状況を考えれば∇ψ ≈ ∇ zψ とできる。また、密度は完全圧縮系のNHMでは細かい振動を含みうるが、ここで表現しているバランスはこのような細かい振動を含まない
総観規模のものであると考えられ、z面上でほぼ一定と扱える（ ρ ( x, y, z ) ≈ ρ ( z ) ）。コリオリパラメータは緯度の関数だが、モデルの対象領域が狭ければほぼ一定と考えられる（ f ( x, y ) ≈ f ）。これらの近似により、次式が成り立つ。

p B ≈ ρ fψp0

図6.4.1 下層における流線関数と速度ポテン
シャルの水平相関。中心が同じ点同士の相関を表す。上が北の方向で、一辺の長さは1200km。

(6.4.12)

総観規模のバランスを表現すると仮定したため、風とバランスする左辺の気圧は動圧成分を含まず、モ

8

同じ点での流線関数と速度ポテンシャルの相関（線形回帰関係）に加え、水平自己相関を考慮することで、同心円状の相関が説明されうる。

重要な役割を果たす。このほか、気圧との相関も考慮できる形として、次式の非バランス成分（非バラ
ンス速度ポテンシャル χ U ）を独立変数として選ぶ。

χ U = χ (u, v) − r2ψ (u, v) − r3 (π B ( p surf , θ ) )c / R
p(6.4.15)
ここで、 r2 及び r3 は回帰係数である。右辺第3項については、水平風と気圧の相関として地衡風関係のみを考える場合は r3 = 0 であるが、低気圧と水平収束の関係や圧縮の効果等も考慮できるようにするため、導入した。
最後にメソ4DVARで扱っている非バランス水平風について触れておきたい。(6.4.8)式から気圧の勾配を直接計算すれば、気圧の勾配と水平風との相関を考えることができる。すなわち、気圧の勾配と地衡風の関係で釣り合う成分をバランス水平風 v B とし、次式で定義される非バランス水平風 v U を独立変数として選ぶ。
⎛rv U = v − ⎜⎜ uu
⎝ rvu

ruv ⎞
⎟ v B ( p)rvv ⎟⎠(6.4.16)

ここで、 r は回帰係数である。地衡風関係は ruu 及びrvv で表されている。NHMの場合は、気圧の摂動
P = ( p − p )G 1 / 2 を用いて、v B として次式を用いることができる（Saito et al. 2001）。なお、 p はz面上で一定となる背景場の気圧であり、G 1 / 2 等は計量テンソルである（吉崎 1999）。
⎛ ∂p ⎞⎜−⎟
1 ⎜ ∂y z ⎟v B ( p) =ρf ⎜⎜ ∂p ⎟⎟⎜ ∂x ⎟z ⎠⎝⎛ ∂P ∂G 1 / 2 G 23 P ⎞
⎜−
⎟
−
1/ 2
1
∂
y
∂
*
G
z
⎜
⎟
=
ρ G 1 / 2 f ⎜ ∂P ∂G 1 / 2 G 13 P ⎟+⎜⎟G 1 / 2 ∂z * ⎠⎝ ∂x∂ρ ∂ρu
∂ρv
∂ρw
+
+
+
=0
∂t
∂x z ∂y z∂z(6.4.18)

静的バランスは(6.4.18)式の定常解により与えられるため、時間微分項を0とすると、運動量の全発散が0という関係となる。z*座標系の微分にすると、この関係は次式のようになる（Saito et al. 2001）。
∂ρG 1 / 2 u ∂ρG 1 / 2 v ∂ρG 1 / 2 w *
+
+
=0
∂x
∂y
∂z *

(6.4.19)

ここで、 w* = dzdt* とした。厳密にはこの式を扱う必要があるが、水平風の運動量の鉛直シアがそれほど強くない状況を考えると、(6.4.18)式の微分をz*座標系の微分と近似的に見なすことができる。これを考慮した上でこの関係を変数変換として表現するため、次式のような非バランス鉛直風を独立変数とする。
⎛ ∂ρu ∂ρv ⎞
⎟dz
+
∂y ⎟⎠
⎝ ∂x

ρwU = ρw − r4 ∫ ⎜⎜(6.4.20)

ここで、r4 は回帰係数であり、第6.4.3項の統計からその値を見積もると、図6.4.2を得た。
(3) その他の変数間相関
ここまでで、風及び気圧という力学的な変数の相関を表現した。残る変数は、温位、水蒸気量及びその他の微物理量である。
温位は、静力学平衡の関係で、気圧と強い相関があると考えられる。静力学平衡の場合、気圧は各層の気塊の重さとの釣り合いにより求められるため、(6.4.17)(6.4.17)式を(6.4.16)式に代入し、第6.4.3項による統計から回帰係数を計算すると、気圧の予報誤差に空間的な細かい振動があり、最大0.01程度 9 と極めて小さくなって、この方法はNHMの場合は適用できないことを確認した。
一方、流線関数を用いると、(6.4.12)式から分かる
ように、気圧の勾配は必要ない。このため、気圧の予報誤差に空間的な細かい振動があっても、その影響は大きくならないことから、JNoVA0では(6.4.13)
式で定義した非バランス流線関数を用いることとし
メソ 4DVARでは最大となる中層で 0.3 程度の値を持つ（石川 2002）。

9

たのである。
(2) 水平風と鉛直風の相関
次に、水平風と鉛直風の相関が、次式で与えられる質量保存の法則（連続の式）により説明されることを考える。

図6.4.2 非バランス鉛直風の回帰係数 r4 の値。縦軸は鉛直層番号、横軸は回帰係数の値を表す。

平面の自由度しか持たない。これを考慮し、メソ
4DVARと同様に、モデル面最下層気圧と各層の温位をあわせて一つの変数と見なし、その鉛直相関を考慮することにより、温位と気圧の相関を表現する。
本来鉛直相関は1変数の鉛直層間の相関のことであり、第6.4.2項に示した②の変換により、固有モードに展開することで表現される。ここでは、モデル面
最下層の気圧を、温位の新たな1層とみなすことにより、1変数として扱う。当然この変数は、通常の変数よりも鉛直層数が1層多い変数となる。すなわ
ち、鉛直層数を n z とすると、次式で定義される θPs を独立変数に選ぶ。なお、 θ (n) はモデル面第 n 層の温位を示す。
⎧θ (n) n = 1, L , n z
θPs (n) ≡ ⎨n = nz +1⎩ p surf

π B ( z ) = π surf + ∫
surf

dz

θ

変

(6.4.22)

水蒸気量 q v については、力学変数や温位との気象学的な関係が明確でないため、変数間相関を0とする。ただし、以下の点に注意が必要である。水蒸気
量はNHMでは混合比として扱われるが、水蒸気の混合比は下層から上層に向かって指数関数的に減少する量である。これをこのまま扱うと、上層の小さなインクリメントが鉛直相関により下層の水蒸気場を不自然に大きく変えるなどの問題が生じる可能性がある。これを未然に防ぐため、全層同程度の量として扱えるように対数をとって扱う方法も考えられるが、ここでは、背景場の気温から決まる飽和水蒸
気量 q vs 0 との比をとって扱う。すなわち、次式で定義される q vU を独立変数に選ぶ。
qv =
U

qv
q vs 0

表 6.4.1(6.4.21)

なお、これから静力学平衡の関係によりバランス気圧が次式のように計算される。
z

量に関する相関を統計的に表現することは一般に困難である。このため、直接推定できるような観測がある場合にそれをその微物理量に解析することは視野に入れうるが、微物理量に関する変数間相関については、ここでは考慮しないこととする。微物理量の同化については、今後の展望も含め、第6.4.6項でもう一度触れる。
(4) まとめ
以上に述べたJNoVA0の独立変数を、表6.4.1にまとめる。なお、NHMでは風の変数として運動量を扱っていることを考慮した。また、気圧は他の変数との静的バランスがないと考えられる動圧成分としている。 D は水平風の運動量の発散である。

(6.4.23)対数関数は非線形関数だが、(6.4.23)式は線形であり、変数変換の際に非線形の影響がないという利点がある。
その他の微物理量については、これまで考慮したような明確な関係を気象学的に与えることはできない。一方、地表付近に強い水平風の収束があればその上空に雲や雨等があることが多く、また、地表で雨が観測されればその上空に雲や霰等があることが多いだろう。これらを考慮すると、微物理量に関する相関は統計的に表現されうるのではないかと考えるかもしれない。しかし、雲底や雲頂の高度、雨や霰の割合等多くの要素が雲の種類などによって異なり、個々の場の状態による影響が強いため、微物理

JNoVA0 の独立変数

数 名

定義式

非バランス流線関数

ψ U = ψ − r1π B

非バランス速度ポテン
シャルχ U = χ − r2ψ − r3π B

非バランス鉛直風

WU = W − r4 ∫ Ddz

cp / R

cp / R

温位－地表面気圧

π U = π −π B
θPs

水蒸気量

qv =

非バランス気圧U

qv
q vs 06.4.5 JNoVA0の構築及び動作結果
これまでの設計で、(6.4.6)式に示した評価関数を計算することができる。また、変数変換のアジョイ
ントを用いて評価関数の勾配が計算でき、更にメソ
4DVARで使われている準ニュートン法による非線形最適化アルゴリズム（L-BFGS）を用いることで、最適解が求まる。
このようにして構築したJNoVA0のプロトタイプ
が、設計通りの動作を示すことを確認した。
まず、モデル第20層（全部で45層）の領域の中心に+1m/sの西風をD値として与えて解析した結果を
図6.4.3に示す。図6.4.3(a)には東西風の解析インク
リメントを示しており、東西方向に長い同心円状に西風が広がり、南北方向に離れた点で東風となっていることが分かる。なお、中心で最大+0.65m/s程度の西風となっている。図6.4.3(b)には南北風の解析イ
ンクリメントを示しており、南北風は第1,3象限に正、第2,4象限に負となっている。東西風と南北風を総合すると、水平風が南北に渦を作るように解析されている。斜め方向の歪みは、流線関数と速度ポテンシ
ャルの相関を考慮した影響である。ここでは示さな
いが、この相関を表す回帰係数 r2 を0とすると、この歪みが現れないことが確認できる。図6.4.3(c)には
水平風ベクトル及び気圧のインクリメントを示している。北側は低気圧性の渦及び低気圧、南側はその逆となっており、地衡風の関係が表現されていることが分かる。図6.4.3(d)には中心での東西風速の鉛直
プロファイルを示しており、統計的に求めた背景誤差の鉛直相関を反映している。図6.4.3(e)には鉛直風のインクリメントを示しており、風下側の上層に上昇流、下層に下降流、風上側ではその逆となっている。第20層から少し下にずれているのは、鉛直風がスタガードに配置され（Lorenz格子）、半層ずれているためである。

(a)

次に、第20層の領域の中心に+1Paの気圧をD値として与えて解析した結果を図6.4.4に示す。同心円状に広がった気圧の回りに、地衡風の関係を満たすよ
うな渦が解析されていることが分かる。

(b)
(c)

図6.4.4 第20層（全45層）の領域の中心に+1Paの気圧をD値として与えた場合のJNoVA0による解析結果。観測
誤差標準偏差は1Paとした。領域は一辺600kmの正方形領域で、格子間隔は一辺5km。

L
H

(d)

(e)

図6.4.3 第20層（全45層）の領域の中心に+1m/sの西風をD値として与えた場合のJNoVA0による解析結果。観測誤差標準偏差は1m/sとした。領域は一辺600kmの正方形領域で、グリッドサイズは一辺5kmとした。 (a)は東西風の第20層における解析結果、(b)は南北風の第20層における解析結果、(c)は第20層における水平風ベ
クトル及び気圧の解析結果をそれぞれ示す。(d)は領域の中心での東西風の解析の鉛直分布を示す。横軸は風速[m/s]、縦軸は鉛直モデル層番号である。(e)は鉛直風の領域の中心を通る東西断面における解析結果を示す。

6.4.6 今後の課題及び展望
気象庁におけるNHM用のデータ同化はJNoVA0が開発の黎明期にあり、今後様々な問題に直面する可能性があるとともに、大きな期待もある。最後に、今後の課題及び展望について、現時点で考えられる
いくつかの点を指摘したい。
(1) 変数変換の設計の妥当性について
これまで述べたように、互いに独立と見なす変数への変数変換①には任意性がある（第6.4.2項）ため、予報誤差の統計調査を行い（第6.4.3項）、NHM用の変数変換を設計した（第6.4.4項）。第6.4.4項(1)で述べた地衡風関係の考慮の際、Poisson方程式の計算を避け計算コストを下げるため、z座標系とz*座標系を同一と見なし、密度がz面上で一定であるといった仮定を置いた。これらの仮定の影響がどの程度のもの
か、調べる必要がある。なるべく仮定の少ない変数変換を設計し、結果を比較することにより、本設計の妥当性を調査する必要がある。
(2) 初期値化の問題（ペナルティ項の設計）
ここまでの議論では、初期値化の問題は考えていない。ここまでの方法で解析された初期値には、不自然な重力波モード等を発達させるような成分が含まれている可能性がある。完全圧縮系のNHMでは音波も表現されるため、不自然な音波モードが現れる可能性もある。これらについては実際のモデルの挙動を見て確認する必要があるが、一般に、これらの不自然なモードの発達を防ぐため、初期値化の手続きが必要となる。
一方、変分法では評価関数の最小化問題としているため、適当なペナルティ項を評価関数に加えることができる。不自然なモードを発達させる成分を初期値から消散させるようなペナルティ項を加えることにより、解析後の初期値化の手続きを省略できる場合がある。ただし、ペナルティ項の追加は、評価関数の形状を変化させ、最適解探索時の収束に悪影響を与えるという欠点がある。ペナルティ項の設計については、実際の事例について不自然なモードがどの程度現れるか調査し、これらのモードを消散させるように設計する必要がある。
最終的にどの手法を用いるかという点については、
解析にかかる時間（計算量）とその予報の結果とを比較検討した上で、初期値化を行うか、どのペナル
ティ項を追加するか、最適な方法を選択する必要があるだろう。
(3) 現象スケールに合わせた解析
非静力学モデルでは、従来の静力学モデルでは表現できないような小さなスケールの現象を表現できる。こういった細かい現象に対する解析を行うための観測データとしては、時空間的に密なドップラー
レーダーデータをはじめ、航空機観測データや、今後重要性が増す衛星観測データなどが考えられる。
一方、小さなスケールの現象は予測可能性が低く、誤差の成長が速い。このため、現象のスケールに合わせた適切な解析には多くの問題が想像される。
典型的な雷雲の時間スケールは1時間程度である。
総観規模の場の状態からは、小さなスケールのシビ
アな現象の潜在的な起こりやすさを知ることはでき
るが、発生の場所や時間、強さは正確にはわからない。高解像度の非静力学モデルではこのような現象も表現しうるが、現実的に意味のある予測をするにはそのスケールに合った適切な初期値が必要である。
現在のように6時間予報値を第一推定値とすることは、この場合適切とは言えない。というのも、1時間程度の時間スケールの現象は、6時間後にはほぼランダムに発生すると考えられるからである（Lorenz 1993）。それどころか、ランダムに発生した現実にはないシビアな現象が第一推定値に含まれ、解析に悪影響を与えることさえありうる。
予報誤差統計情報の取り出し方や4次元変分法における同化タイムウィンドウの長さも問題となる。
現象の規模に応じて誤差相関スケールや予測可能な時間スケールが異なるが、これを考慮した取り扱いが必要だろう。あるスケールの現象に着目した統計
情報をどのように取り出すか、4次元変分法においてどのようなスケールに着目してタイムウィンドウの長さを設定するか、問題である。
観測データの取り扱いについても、単純ではない。
確かに、時空間的に細かい観測データは、細かい現象を表現しうるNHMにおいて、より多くの情報を提供する。しかし、観測の時空間的な代表性の違い
やデータの間引き方、観測誤差相関が無視できない点など、新たに考慮すべき問題がある。代表性の違いにより観測誤差を推定し直す必要があるし、デー
タの間引き方も新たに開発する必要がある。観測誤
差共分散の非対角成分が無視できない場合は、背景
項と同様、観測項についても非対角成分の扱い方に工夫が必要かもしれない。
以上に述べたように、スケールに応じて解析の取り扱いが異なるため、まず総観規模の解析を行い、その後それにメソスケールの解析を埋め込む、といった方法も考えられる。モデルで表現しうるすべてのスケールについて、解析を一度に行うことがふさわしいのか、その可能性も含め、細かい現象を表現しうる非静力学モデル特有の問題であり、今後の研究を待つ必要があるだろう。
(4) 微物理量の扱いモデルで微物理過程を詳細に表現するようになり、より現実に近い雲や降水の形成過程を表現できるよ
うになる一方、これに合わせた雲水や雨、霰といった微物理量の解析についてはまだ実用段階にはなく、今後の課題となる。
それぞれの雲の微物理構造は様々であると考えられるため、微物理量の解析の際、統計的な背景誤差共分散を用いることはふさわしくない。3次元変分法ではモデルのバランスを動的に扱うことができな
いため、微物理量の同化には特別な扱いが必要となる。これについては、背景誤差共分散では独立して扱い、観測演算子で積雲のプロファイル等を仮定して扱うといった方法や、観測前処理でPI（物理的初期値化）のような考え方で擬似観測データを生成するといった方法が考えられるが、現時点では良い方法が知られていない。
微物理量を直接推定できる観測については、レー
ダー反射強度、衛星センサー等があるが、これらの扱いも含め、今後の大きな課題の一つである。
(5) 3次元変分法の様々な運用形態
3次元変分法は、4次元変分法と比較して計算コス
トが格段に低いため、この点を利用した様々な運用形態が考えられる。毎時解析サイクルの実現や、PC上で動作するものの構築も可能だろう。前者に関しては、次のような利点がある。3次元変分法では原理的に解析時刻に近い時刻の観測しか有効に使うことができないが、頻繁に解析を行うことにより、より多くの観測データを有効に生かすことができる。
また、実況監視を目的とした現在の毎時下層風解析

を3次元変分法により高度化する可能性も検討されており、この場合、特にドップラーレーダーの動径
風を直接同化できる点で、現行の最適内挿法よりも優れている。
後者に関しては、領域を限定することでNHMをPC上で動作させることが可能であり、JNoVA0についても同様にPC上で動作可能であることから、
NHMのようにソースコードを提供することで、気象台等のより多くのユーザからのフィードバックにより、問題点の抽出とプログラムの安定化や改良に資することが期待される。

参考文献
Fisher, M., 2001: Assimilation techniques (3):3DVAR. Meteorological Training Course
Lecture Series, ECMWF, 11pp.
Haltiner, G. and R. Williams, 1980: Numerical
Prediction and Dynamic Meteorology. 2nd ed.,John Wiley and Sons, p.15.
石川宜広, 2002: メソ4次元変分法. 数値予報課報
告・別冊第48号, 気象庁予報部, 37-47.
Kalnay, E., 2002: Atmospheric Modeling, Data
Assimilation and Predictability. Cambridge
University Press, p.62.
気象庁予報部, 2002: 変分法データ同化システムの現業化. 数値予報課報告・別冊第48号, 87pp.
Lorenz, E., 1993: The Essence of Chaos.
University of Washington Press, 240pp.
Parrish, D. and J. Derber, 1992: The National
Meteorological Center’s spectral statisticalinterpolation analysis system. Mon. Wea. Rev.,
120, 1747-1763.
Saito, K., T. Kato, H. Eito, C. Muroi, 2001:
Documentation of the Meteorological ResearchInstitute / Numerical Prediction Division
Unified Nonhydrostatic Model. Tec. Rep. MRI,42, p.10.
露木義, 1997: 変分法によるデータ同化. 数値予報
課報告・別冊第43号, 気象庁予報部, 102-165.
吉崎正憲, 1999: 地形に沿う座標系における支配方
程 式 . 気 象 研 究 ノ ー ト , 196, 日 本 気 象 学 会 ,37-43.