統計数理（2006）
第 54 巻 第 2 号 247–264
c 2006 統計数理研究所


特集「予測と発見」［総合報告］

日本海循環データ同化モデル
1

広瀬 直毅 ・万田 敦昌

2（受付 2005 年 11 月 1 日；改訂 2006 年 8 月 21 日）

要旨

大規模な地球流体推定問題の先駆例として，日本海のデータ同化・予報モデルを構築した．
高解像度の海洋循環モデルに各種観測データを適切に反映させるため，現実的な計算量で精度を確保できる同化手法を開発した．海面の風応力や熱収支の時空間変化だけでなく，海底地形も制御した結果，日本海表層から亜表層の水温や流動場をかなり現実的に再現することができた．予報結果は，流出重油やエチゼンクラゲの漂流計算にも利用されている．
キーワード： 日本海，海洋循環モデル，衛星観測データ，データ同化，海況予報．
1.

はじめに

海洋・気象といった地球規模の物理現象を研究対象とする場合，必ずといっていいほど問題となるのが，観測データの不足である．相手が巨大であり，またその時空間変動スケールも多岐に渡るため，ある事象を解析するために必要十分な量の観測値を得ることは，多くの場合，困難である．
一方，コンピュータの高速化により，かなり現実的な地球流体の変動を再現することが可能になってきたが，まだその解像度も不十分である．例えば，新世紀初頭に世界最大・最速であった「地球シミュレータ」をもってしても，大気海洋結合モデルの格子は約 100 km，海洋単独の循環モデルでようやく 10 km 程度である．特に沿岸域では，漁業活動や防災の観点から高い解像度が要求される．非静水圧の乱流混合現象など，微小な変動の影響をパラメータ化する手法も確立しておらず，将来や過去はおろか，現在の大気・海洋の状態をどれだけ正確に把握できているか，その精度も疑わしい．
このような観測・モデルの問題点を解決する有力な手段と期待されているのが，データ同化手法である．データ同化とは，地球物理学を中心とした用語で，より一般的には逆問題を解く
手法 inverse method と認識される．断片的な観測データを数値モデルへ適切に融合（同化）させることにより，より現実的な地球流体の 4 次元的な変動を推定することが可能となる．
本研究では，大規模な海洋流体推定問題の先行モデルとして，日本海を取り上げる．日本海は太平洋北西部に位置する縁辺海の一つにすぎないが，太古の昔は湖であったといわれる通り，
現在でもかなり閉鎖的な系を成す．例えば，図 1 に示すように，日本海南西部から流入した対
馬暖流の大半は津軽海峡・宗谷海峡から流出するが，北部には独立した亜寒帯循環が形成されている．中層・深層の水塊は，ほぼ 100％この亜寒帯循環域で冬季に形成されていると考えられており，太平洋の同じ深度と比較してかなり冷たい．
1
2

九州大学 応用力学研究所：〒816–8580 福岡県春日市春日公園 6–1
長崎大学 生産科学研究科：〒852–8521 長崎市文教町 1–14

248

統計数理 第 54 巻

図 1.

第 2 号 2006

長沼（1977）と Yarichin（1980）による夏季の日本海表層循環の概念図．

暖水渦や前線波動などの中規模変動も激しく，海氷まで生成するこの海域は，水平スケール
さえ無視すれば，大洋の特性と大差なく，しばしば「ミニ大洋」と称される．領域が小さく観測活動を行いやすいため，外洋と比較するとかなり現場観測の密度が高い．
領域が小さいおかげで，先駆的な数値シミュレーションや同化モデリングも試行しやすい．
特に，大規模な地球流体モデルに対して，高精度な非線形フィルタリング手法を適用するのは非現実的であり，精度と負荷が両立する適切な同化手法を探求する意味でも，日本海循環は格好の実験場となる．本報告では，データ同化手法を利用した日本海循環の再現・予報モデルを構築する中で，過去数年間に得られた知見を紹介していきたいと思う．
2.

日本海循環モデル

本研究で使用する海洋循環モデルは，九州大学応用力学研究所で開発された RIAM Ocean
Model である
（Lee et al., 2003）．このモデルは球面座標系において，ブシネスク近似，静水圧平衡，自由表面を仮定したナビエ ストークスの方程式，連続の式，移流拡散方程式を有限差分法
で解く地球流体モデル，一般的に海洋大循環モデル（OGCM; ocean general circulation model）と呼ばれるものである．主な予報変数は 3 次元流速ベクトル，海面高度，温度，塩分（濃度），乱流運動エネルギーなどである．運動量の移流項は一般化された荒川スキームと斜め移流で差分化され，エネルギーとエンストロフィーが同時に保存される（Ishizaki and Motoi, 1999）．斜め
移流は，海底地形における水平運動量の鉛直移流をより正確に再現する．格子配置は Arakawa
B-grid である．
水温・塩分の移流項には 3 次精度の Modiﬁed Split-QUICK scheme（Webb et al., 1998）を使用し，粘性項は 4 次精度（biharmonic）で差分化している．また，海洋中の物質は等密度面方向に拡散しやすい傾向があるため，拡散項はテンソル変換によって等密度方向の拡散が卓越するように工夫した（Gent and McWilliams, 1990）．さらに鉛直混合過程は Noh and Kim（1999）の混合層スキームで計算される．

249

日本海循環データ同化モデル

モデルの計算領域は，126.5◦ ∼ 142.5◦ E，33◦ ∼ 52◦ N の日本海全域で，水平の解像度は 1/12◦とした．浅海域の流動場の再現と予報に興味があるので，鉛直解像度は，0∼290 m 深までを
5∼20 m 間隔と密にし，290 以深を粗く，25∼600 m 間隔に設定した．鉛直の総層数は 36 層である．海洋流動場に対して海底地形は決定的に重要な影響を及ぼすので，後述するように，様々な海底地形を試行し，地形データを最適化する．陸地部分を除いた本モデルの実質的な自由度
は 107 程度となる．
風応力によって海面から与えられる運動量と海面を通じた熱交換・蒸発・降水量は，基本的にヨーロッパ中規模気象センター（European Centre for Medium-range Weather Forecast，
ECMWF）の長期月平均を使用した．2002 年以降の予報実験に対しては，より解像度の高い気象庁データ（GPV MSM）を用いる．しかし，境界条件や表層混合過程の不正確さに起因して，海洋水温や塩分を必ずしも正確に再現することができない．表層の水温・塩分分布は内部の海
洋循環に強く影響するため，しばしば海面水温・塩分を観測値へ緩和（漸近）する措置が取られる．この緩和法は最も単純な部類のデータ同化（ナッジング法）と理解することもでき，その妥当性について，後に検証する．
対馬・津軽・宗谷海峡を通じた流入・流出として気候月平均の流速・水温・塩分観測値を与えている．水温・塩分は主に JODC，KODC（Japan or Korea Oceanographic Data Center）で蓄積された船舶観測データから，流速はフェリー観測値に基づいている．
計算時間は Fujitsu VPP-5000（ピーク性能 9.6GFlops/cpu のベクトル並列計算機）により，1
年間の時間積分に必要な計算時間は約 16 時間であり，MPI を利用して 16 cpu 並列計算を行う

(a)

(b)

図 2. （a）衛星高度計観測と（b）数値シミュレーションによって得られた海面高度の標準偏差（cm）．

250

統計数理 第 54 巻

第 2 号 2006

ことにより，かなり線形的な並列化効率を得ることができた（1 年間の時間積分に計算約 1 時間）．外洋のシミュレーションと比較して，かなり計算時間は短い．同じ解像度であれば，当然，領域が小さいほど計算コストも節約される．
海面高度変化の標準偏差を図 2 に示す．まだデータ同化していない段階で，変動強度の分布
が衛星観測とよく対応していることがわかる．日本海内部の流動場や渦変動がよく再現されていると類推される．衛星観測の標準偏差がより大きいのは，観測誤差のためか，日本海全体の海面変化のためではないかと考えている．即時的予報のために，高度計データとしては誤差補
正の不十分な速報値 I-GDR を使用している．また，この日本海モデルにおいては領域全体の体積保存を課しているため，日本海の平均水位は診断的に得るしかなく，過小評価になりやすい．例えば，気圧変化に伴う水位の非平衡変化も大きいと報告されている（Inazu et al., 2006）．
3.

海底地形の最適化

データ同化の最初の例として，海底地形の推定問題を取り上げる．既知のデータとして利用される機会の多い海洋水深であるが，海面高度計の解析から南大洋に新たな海山が発見される（Smith and Sandwell, 1997）など，未だその精度は不確かである．
日本近海においても，各種水深データセット間の差は顕著である．Hirose（2005）は，海流の分布に対して直接的に海底地形が影響すると考えられる浅い対馬海峡を例に，水深データの精度を検証し，より正確な推定値を提案した．
まず，代表的な海底地形データとして，DBDB-V（Sandy, 1996），ETOPO2（Smith and Sandwell, 1997），SKKU（SungKyunKwan University, Choi et al., 2002），JTOPO1（海洋情報研究センター, 2003）の 4 種類を比較した．ETOPO5 は，海岸線付近を除いてほぼ DBDB-V と同一であったので，使用していない．
図 3 に示すように，対馬海峡のような沿岸浅海域でも各データセット間の差が大きい．特に，ETOPO-2 では他と異なる小スケールの海盆・海山が目立つ．逆に，DBDB-V はかなり滑らかな分布傾向である．33-36◦ N，128-131.5◦ E の 200 m 以浅における，各データセット間の rms 差を取ると，10 m∼20 m にも達する．対馬海峡の平均水深が 100 m 程度であることを考えると，地形データの誤差は 10∼20％程度と推察できる．
次に，それぞれの地形データを用いて前出の 1/12◦ 日本海循環モデル（RIAM Ocean Model）
を実行し，計算された流速を，日韓往復フェリーに設置した ADCP（acoustic Doppler current
proﬁler）で長期間観測された海流の平均値（Takikawa et al., 2005）と比較する．
図 4 に，観測データとモデル計算結果による冬季平均流速分布を示す．この海域では，対馬
暖流が平均的に北東方向へ通過している．特に冬季は，海面冷却によって水質や流速場も鉛直方向に均質化しているので，より直接的に海底地形の影響を受ける．
DBDB-V を用いた場合は計算結果と観測値との差が大きいが，他の 3 データを用いればある程度観測値に近く，それぞれの水平速度分布も大局的には類似していることがわかる．だが，沿岸近くの速度ベクトルの方向や強さなど，細部では計算結果の間に差が認められる．RMS 差にして 3 ∼ 5 cm/s となり，対馬海峡における冬季平均流速（10 ∼ 20 cm/s）と比較して，無視できない大きさである．
Losch and Wunsch（2003）は，擬似観測された流速データを満足するような海底地形を，逆問題の解として推定する問題に取り組んだ．本研究では，冬季の平均流速データとの rms 差を最小にする最小二乗問題を定義し，各データセットの線形和として適切な海底地形を求めてみた．
(3.1)

y ∼ Ax

ここで，ベクトル y は観測された流速データ（図 4（a）），行列 A は図 4（b）
（
，c）
（
，d）
（
，e）に示

日本海循環データ同化モデル

図 3. （a）DBDB-V，（b）ETOPO-2，（c）SKKU，（d）JTOPO1 の各データセットによる対馬
海峡海底地形の分布．

図 4.

対馬海峡におけるフェリー観測と海洋循環モデルによって得られた長期間冬季平均流速
分布．モデル計算には 4 種類の海底地形を試行した．

251

252

統計数理 第 54 巻

第 2 号 2006

表 1.

フェリー観測とシミュレーションの流速差を最小にする最適重み係数．

図 5.

最適加重平均（a）E1，（b）E3 によって求めた対馬海峡海底地形の分布．

すモデルの流速分布を並べた応答関数，ベクトル x が求めるべき重み係数である．加重平均として解を得るために，
(3.2)



xi = 1(i = 1 ∼ 4)iの条件も付加した．
本研究においては，特異値分解によって行列 A の一般化逆行列を求め，最小二乗となるべ
き解 x を得た．類似した地形データを使用するため，多重共線性も疑われるが，式（3.2）の拘束条件のため自由度が一つ下がり，相関係数行列の行列式は 0.40，最小固有値は 0.25 となり（4データセットすべてを用いた場合; E1），とりあえず安定な解を得ることができる．
得られた最適重み係数を表 1 に示す．結果として，日韓のデータセットが重視され，DBDB-V
や ETOPO2 には負の重みが推定された（E1，E2，E4，E5）．特に DBDB-V は他の 3 地形データ
より誤差が大きいと予想されており，DBDB-V に対する強い負加重（E1 の −0.24，E5 の −0.32）
は信頼しがたい．正の加重平均だけに限定すると，JTOPO1 に 3/4，SKKU に 1/4 の重みで作成した水深データ（E3）が最も正確で，観測値との rms 差は 3.28 cm/s である．
観測値との rms では，SKKU より ETOPO2 の方がわずかに正確であったにもかかわらず，
線形結合で ETOPO2 にあまり加重されないのは（例えば，E1，E2），JTOPO1 と ETOPO2 に相関（共線性）があるためだろう．実際に，JTOPO1 と ETOPO2 を結合しても，有意な加重は得られない（E4）．両者とも同じインバース法を使用して，海面高度分布から海底地形を推定しているためと考えられる．従って，説明変数として JTOPO1 あるいは ETOPO2 のどちらか一方を採用すれば十分であり，結局，JTOPO1 と SKKU（E3），または ETOPO2 と SKKU（E6）のような線形結合が推奨される．

日本海循環データ同化モデル

253

最適化された海底地形 E1，E3 を図 5 に示す．JTOPO1 データへの強い重みを反映して，
JTOPO1 水深データとよく似た分布をしている．だが，JTOPO1 と比較して，やや平滑な分布を示しているのは，SKKU データなどの影響である．
海洋流動の非線形性のために，上記の線形解析は必ずしも保証されるものではないが，実際に 3/4 JTOPO1 + 1/4 SKKU の加重平均地形でモデル計算すると，期待通りの流速分布が再現されたので，冬季平均流速と海底地形の間にはほぼ線形関係が成立している，といえる．
一般的に，浅海域で再現される流動場は，海底地形に強く依存する．この物理関係を逆に利用し，流速データから海底地形を推定することは，まさに逆問題を解く作業である．データ同
化とは，地球物理学で広く認知された用語であるが，その本質は逆問題である，と理解するのが妥当であろう．
4.

海面水温データの同化

フェリー観測では，その往復線に沿って海洋環境諸量を断面的にモニタリングすることができるが，水平的な代表性に乏しい．地球上の 2 次元平面を最も広く観測することができるのは，
人工衛星による遠隔観測であろう．特に，海洋変動に関しては，海面水温と海面高度の観測が長期間に渡って行われている．本節と次節では，その両者の適切なデータ同化手法を検討する．
まず，海面水温の変化は，表層混合の力学を強く反映する．衛星海面水温データを数値モデ
ルに同化させることによって，海洋混合層の再現性が向上すると期待できる．しかし，対流混合過程においては静水圧平衡がほとんど成立しないなど，地球流体力学の基本的な近似が通用せず，同化手法の線形化も困難である．正確にデータ同化を行うためには，非線形・非正規変化を考慮した同化手法が要求されるが，OGCM では常にその膨大な計算量が問題となる．そこで，
Manda et al.（2003, 2005）では，鉛直 1 次元の乱流混合層モデル（Noh and Kim, 1999）に対して，海面水温データを Ensemble Kalman ﬁlter（EnKF; Evensen, 1994）や sampling importanceresampling ﬁlter（SIR; Gordon et al., 1993）などで同化する実験を行い，簡便なナッジング法との比較を行った．
ナッジング法とは，モデル値を観測値へ徐々に漸近させる同化手法（緩和法）である．海面水図 6.

海洋観測定点 Papa（50◦ N，145◦ E）における鉛直 1 次元データ同化モデル結果と観測値との rms 差（◦ C）．〇実線，△破線，▽点線，□一点鎖線はそれぞれ Run 1（シミュレーション），Run 2（EnKF），Run 3（SIR），Run 4（ナッジング法）の実験結果を示す．

254

統計数理 第 54 巻

第 2 号 2006

温の場合は，海面の熱境界条件へ反映させて，
Q = Q∗ + ρCp ∆z1 (Ts∗ − Ts )/τ

(4.1)のように海面を通して大気と交換される熱量データ Q∗ を，海面水温の計算値 Ts を観測デー
タ Ts∗ へ，徐々に（時間スケール τ で）漸近させる．ここで，ρ，Cp ，∆z1 は，海水の密度，比熱，モデル鉛直 1 層目の厚さである．緩和項の係数部分 ρCp ∆z1 /τ が大きい（小さい）ほどより強く（弱く），モデル値が観測値へ漸近することになる．EnKF や SIR では，この緩和項の代わりに，全支配変数（水温・流速・乱流混合エネルギー等の鉛直分布）の誤差変化の非線形性（非正規性）まで評価して適切な修正を行うため，高精度なデータ同化が期待できる半面，計算負荷が桁違いに大きくなる欠点がある．
Burchard et al.（1999）
を参考にして，海洋観測定点 Papa（50◦ N，145◦ E）
における 1961/3/25∼
1962/3/24（1 年間）の風応力，熱交換データを用いて，Noh and Kim（1999）の混合層モデルを駆動した．定点 Papa における現場海面水温データに標準偏差 0.6◦ C の白色誤差を加え，これを擬似的な衛星海面水温データとして同化に使用した．
海面水温データを同化せず，単純に風応力と熱交換でモデルを駆動した場合（Run 1），水温の海面から水深 80 m までの平均自乗誤差（rms error）は 0.75◦ C となった（図 6）．EnKF と SIR
によって海面水温データをモデルに同化した場合，それぞれ rms 差は 0.42，0.43◦ C となり，シ

図 7.

気象庁日本海ブイにおける（a）1m，（b）50m，（c）100 m 深における水温の時間変化，赤線は水温の観測値，黄，青，緑線はそれぞれ 0.1，1，10 日の時定数を用いた場合のモデル水温を示す．

日本海循環データ同化モデル

図 8.

255

気象庁ブイにおける水温のモデルと観測値の rms 差（◦ C）．

ミュレーションに比べ約半分に減少した．ナッジング法を用いた場合でも，重み係数（時定数）を最適化することで rms 差は 0.45◦ C まで減少し（時定数約 0.5 日の場合），EnKF や SIR と比べてそれほど劣らない結果となった．計算コストを考えると，ナッジング法は実用的で適当な精度を持った衛星水温データ同化手法であるといえる．
ナッジング係数に対する依存性をさらに詳しく調査するため，同じ Noh 型の混合層モデルを組み込んだ日本海海洋循環モデル（RIAM Ocean Model）に対して，衛星水温データを同化する．
実験期間は 1997∼2000 年で，ECMWF の予報値を日平均して海面の運動量や水・熱交換の境界条件とする．さらに，A-HIGHERS 海面水温データ（Sakaida et al., 2000）を日平均し，3 種類の時定数（0.1 日，1 日，10 日）でナッジング法により同化した．なお，EnKF や SIR によって，大量の衛星水温データを 3 次元循環モデルに同化することは，計算負荷の面から困難である．
同化の効果を確認するため，計算結果を気象庁ブイの水温観測データと比較する（図 7）．海面付近では，短い時定数（0.1 日，1 日）で強く海面水温を衛星観測値に緩和させる方が，現場観測にも近い結果が得られる．しかし亜表層（50, 100 m）においては，時定数が短すぎる（0.1 日）と，逆に差異が広がってしまう．
現場観測による水温データとの rms 差をとると，その傾向がはっきりする（図 8）．衛星海面
水温データを短い時定数によって強く同化することにより，海面付近での温度の再現性は向上
するが，逆に亜表層（50, 100 m）での rms 差は広がってしまう．全体として適切なナッジング係
数は，約 1 日程度であり，定点 Papa の 1 次元モデルによる実験結果と矛盾しない結果である．
鉛直 1 次元の実験結果がそのまま 3 次元計算でも保証されるわけではないが，少なくとも海面を通じた熱交換条件が表層混合の主因である限り，ナッジング法は計算の精度と負荷を両立する現実的な海面水温データ同化手法であるといえよう．ただし，最適なナッジング係数は解像度や同化間隔などによって異なるため，問題に応じて注意深く選択する必要がある．
5.

海面高度データの同化

海面水温データは，ごく表層の温度を直接代表する観測値であるが，海洋内部の変化をより広く反映する観測値として，海面高度データが挙げられる．海洋内部の密度変化や流動変化によって海面が上下する関係は，比較的簡単な熱膨張や地衡流の関係によって説明することができ，逆に海面高度データから内部変化を推定することも期待できる．
近年では，Jason-1，Envisat，GFO など，複数の人工衛星によって海面高度計データが提供
されており，かなり広範囲に渡って，しかも 10 日程度の等時間間隔で海面の凹凸を知ることができる．本研究では，アメリカ海軍研究所（U. S. Naval Research Laboratory）で統合化されている海面高度計データを取得し，200 m 以浅のデータを除去したり，標準偏差の 3 倍を超え256

統計数理 第 54 巻

第 2 号 2006

る異常値を除去するなどの品質検査を施した．また，モデルの格子間隔に合わせ，1/12◦ ，1 日
毎に平均値を求めてデータ同化に使用する．
同化手法として，拡張カルマンフィルターを採用する．しかし，107 の自由度を持つモデルに対して，厳密にカルマンフィルターを適用すると，膨大な計算量を消費する．そこで本研究
では，Fukumori et al.（1999）などを参考に下記のような近似を行って，カルマンフィルターを現実的な計算量で実行する．
（1）
衛星高度データの平均分布（ジオイド面）は特に大洋縁辺部で信頼できないので，観測デー
タの長時間平均値（バイアス）は同化に使用せず，モデルの長時間平均海面高度に修正は加えない．さらに，モデルの時間平均値を基準として各システム行列（状態遷移行列 A，外力行列 G，
観測行列 H ）を作成する．つまり，拡張カルマンフィルターは，時間的に変化しない基準値 xからの摂動場のみを修正することになる．
（2）漸近定常な誤差共分散行列を使用する．人工衛星による高度観測はかなり規則的である
ので，実質的に観測行列 H を一定値として定義できる．遷移行列 A や外力行列 G も時間不変とみなせば，しばらくリカッチ方程式を時間積分すると（倍化法を利用），多くの場合，漸近定
常値を得ることができる．結局，カルマンフィルターより原始的なウィーナーフィルターに近い状態になる．
（3）水平格子を間引きし，1/3◦ 間隔で行列を作成する．モデルの水平格子間隔は 1/12◦ であるが，有限差分モデルにおいては短波長の再現性が悪く，実質的な解像度は水平格子間隔の数倍以上である．1/3◦ 毎に計算しても，ほとんどのモデル変化に対応した誤差を推定することができると考えられる．なお，低解像度のシステム誤差を元の解像度へ変換するために，最適内挿法（Bretherton et al., 1976）を用いる．
（4）鉛直方向の層数も圧縮する．大規模な海洋風成循環は数少ない鉛直モードに縮退していることが多いので，本研究でも，鉛直密度分布を固有値・固有ベクトル分解し，水平速度の上
位 2 モード（順圧モードと傾圧第一モード）と密度分布の最上位モード（傾圧第一モード）のみを，
カルマンフィルターによって制御される縮小モデル xBT ，xBC として定義する．近似（3）と併せて，順圧モードと傾圧第一モードそれぞれ 3000 以下の自由度で縮小モデルを構築することができる．
（1）∼（4）の近似をまとめて表現すると，(5.1)

x ≈ B BT xBT + B BC xBC + x(5.2)

y ≈ Hx + rとなる．x が海洋モデルの水温，塩分（濃度）や流速を表す状態ベクトル，B BT ，B BC がそれ
ぞれ順圧・傾圧モードの偏差を，鉛直モード積分や最適内挿によって元の状態ベクトルへ射影（変換）する 行列，y は観測値，r が観測誤差である．なお，観測誤差は特に近似しない．同様にして，システム誤差共分散行列も(5.3)

P ≈ B BT P BT B TBT + B BC P BC B TBCのように近似される．
さらに，縮小場における状態遷移行列を A ≈ B −I AB ，外力行列を G ≈ B −I G，観測行
列を H  ≈ HB と定義することにより，カルマンフィルターによって記述される一連の定式（リカッチ方程式など）を，ほぼそのまま縮小場でも利用することが可能である（Fukumori andMalanotte-Rizzoli, 1995）．

日本海循環データ同化モデル

図 9. プロセスノイズ（海面風応力）の相関スケールと同化結果（1 日予報値）の検出分散量（海面高度）の関係．検出分散はデータ分散と残差分散の差で定義される．

図 10.

257

プロセスノイズ（海面風応力）の平均値と同化結果（1 日予報値）の検出分散量（海面高度）の関係．

縮小場における各システム行列 A ，G ，H  は，仮想変位の原理を利用して数値的に作成した（Fukumori et al., 1993）．しかし，初期摂動に対して多様な時間スケールの波動が生じ，特に同化間隔以内の短周期変化がエイリアシングするため，摂動を一定期間与え続け，応答も同じ期間だけ平均化する．特に今回問題となったのが外部重力波である．1 日間のモデル時間発展に対応する状態遷移行列を作成するため，単位摂動を 1 日間与え続け，次の 1 日間のモデル応答を平均して，対応する各行列の列ベクトルとした．この措置により，1 日以内の高周波変動（そのほとんどは慣性重力波）のエイリアシングを防ぐことができ，制御可能な遷移行列 A を作成することができた．
海洋循環を駆動する主要因は海上風であるので，縮小モデル xBT ，xBC を制御するプロセス
ノイズも風応力の形で定義した．従って，外力行列を作成する際は，海面風応力に微小な摂動を加え，遷移行列と同様に 1 日平均して短周期変動を除去した．
順圧モードには海面高度が陽に含まれるため，観測行列 H を単純に定義できるが，傾圧第一モードに対しては，密度変化を鉛直方向に積算し，海面高度変化相当量を求めて，観測行列の各列を定義した．やはり，各モードに単位量の変化を加えて，基準値との差を観測行列の各列とする．
非同化のシミュレーション結果と海面高度計データとを比較し，covariance matching（Fu etal., 1993）に基づいてプロセスノイズと観測誤差の共分散行列を定量化したが，さらに同化結果を改善するため，予測値（prediction）の残差分散を最小化するようにチューニングした．例えば，図 9 に示すように，プロセスノイズの空間スケールは，約 100 km 程度で最適値を示すことがわかった．ここで，検出分散（accounted variance）あるいは説明分散（explained variance）とは，データの分散（data variance）からモデルとデータの残差分散（residual variance）を差し引いた量 < d2 > − < (h − d)2 > である．残差分散は unexplained variance といえるので，残りは説明（検出）できた分散，とみなすこともできるだろう．日本海の表層循環や渦構造の典型的な大きさから判断して，妥当な相関スケールである．もちろん，プロセスノイズの大きさ自身も最適化した（図 10）．その結果，約 0.6 dyn2 /cm4 で極大値を得た．海面変化に換算して 3 ∼ 4 cm2程度である．一方，観測誤差共分散行列は 50 cm2 程度の大きさなので，データ同化による修正

258

統計数理 第 54 巻
表 2.

第 2 号 2006

シミュレーション（simulated），同化予測値（predicted），同化修正値（filtered）．それぞれの海面高度変化の分散（model variance）とデータとの残差分散（residual variance），
説明分散（explained variance），相関係数（correlation）．分散の単位は cm2 ．

図 11.

海面高度における検出分散の水平分布（cm2 ）．1/3◦ 格子毎に平均した．

はかなり小さいと期待される．
表 2 に，プロセスノイズを最適化して得た同化結果の統計量を示す．モデルの海面変動量（状態推定値の標本分散）は，シミュレーションでも同化結果でも約 45 cm2 とほとんど変わらず，
観測データの分散 83 cm2 と比較して，積極的な修正を行っていないことがわかる．しかし，変動量は同等にも関わらず，海面高度観測データとの残差分散が小さくなり，検出分散と相関係数がそれぞれ向上しているので，データ同化の効果を定量的に認めることができる．また，1
日毎にデータ同化を実施した際の予測値と修正値（ﬁltered estimates）の残差分散が 4 cm2 程度と小さく，データ同化による修正量は決して大きくない（海面変動の分散と比較して 1 割弱）．
予測値が修正値と大差なく，一方シミュレーションより大幅な改善が見られるので，同化効果にかなり持続性があり，長時間予報が可能と期待される．
図 11 にモデルによって説明することができた海面変動量を示す．非同化のシミュレーションの場合は，日本海に流入した対馬暖流の影響が強い日本沿岸部や朝鮮半島東岸で説明分散が大きく，内部領域では相対的に説明量が小さいが，高度計データ同化によって中規模渦や流動構造が修正され，特に南部中央海域で観測値と一致する成分が増加している．北部の亜寒帯循環域でも修正の効果が認められるが，修正値（ﬁltered estimates）と予測値（predicted estimates）の差（図 11（b）と（c）の差）がある程度認められることから，その持続性は短いようである．北方
冷水域では成層が弱く，変化の早い順圧成分が卓越するため，風応力などの強制力に対する応
答時間が早く，相対的に同化効果が早く消えてしまうようだ．逆に，南部の対馬暖流域においては，傾圧変化が主役であるので，その表層循環や渦活動において，データ同化の効果は数ヶ月間持続可能という解析結果を得ている（Hirose et al., 2005）．北海道から樺太沿岸で検出分散が小さい．津軽・宗谷海峡における流出条件に問題があるのかもしれない．

日本海循環データ同化モデル

図 12.

図 13.

259

船舶観測による 2005 年 7 月初旬の 200 m 深水温分布（◦ C）．水産総合研究センター
日本海区水産研究所による日本海漁場海況速報．

本研究のシミュレーションと高度計同化による 2005 年 7 月 4 日の 200 m 深水温推定値（◦ C）．

データ同化の効果は，海面高度だけでなく，内部でも検証すべきである．一例として，図 12・
13 に 2005 年夏季の水温分布を示す（水産総合研究センター日本海区水産研究所，2005）．非同化のシミュレーションでは，5◦ C 以上の暖水域がかなり広く分布していたが，高度計データの同化により，ある程度，現場観測値に近い推定値が得られている．特に，能登半島北部の 10◦ Cを超える暖水渦や隠岐諸島東部の冷水塊などでその差が顕著である．このように，検出分散が高い海域（図 11）において，亜表層（100∼200 m）の水温分布が修正されていると期待できるので，日本機関の観測が及ばない韓国水域でも同様により現実的な状態が推定できていると考えてよいだろう．日本の経済水域内においても，燃料代や人的コストなど，現場船舶観測には限界があり，時空間的な空白を埋める手段として，データ同化に寄せる期待は大きい．
以上，フェリー観測データや衛星観測データを数値モデルに同化し，推定精度を向上してきたが，まだ他にも利用すべきデータが多く残されている．例えば，船舶や中層フロートなどによる鉛直観測データをモデルへ同化することにより，海洋深層の変化をより現実的に再現・予260

統計数理 第 54 巻

第 2 号 2006

報できるはずである．例えば，本節ではあまり推定精度が向上しなかった北部の亜寒帯循環も，中層フロートの移動速度を利用すれば，かなり改善されることが示されている（笹島 他, 2003）．
6.

海況予報

結局，現在の日本海データ同化モデルでは，上述してきたように最適化された海底地形を用い，衛星海面水温データ NGSST（東北大学）を時定数 1 日で緩和している．さらに衛星海面高度計データ（T/P，Jason-1，ERS，Envisat，GFO）を近似カルマンフィルターで同化し，現実的な中規模渦や循環流をモデルに再現し，予報初期値とする．
予報計算を行うにあたり，将来の気象条件を知り得ないので，予報期間（現在から 10 週間後まで）の海面境界条件として 1 年前の値を使用する．気候値（過去の長期平均値）に切り替えるよりも，高低気圧の通過に伴った総観規模の気象条件を陽に与えることができ，時間的不整合も小さいと考えられる．なお，有意な予報期間（persistency）は現行モデルでは未調査であるが，水平 1/6◦ の準渦解像モデルの海面変化では約 2ヶ月であった（Hirose et al., 2005）．渦解像モデ
ルを用いた他の調査結果でも，同様に約 2ヶ月の予報可能期間を得ている（Komori et al., 2003;
Miyazawa and Yamagata, 2003）．ただし，海洋内部の水温・塩分や流速においては，表層ほど予報可能期間が短く，深層ほど長い．短周期の気象変化の影響を受ける表層ほど短く，逆に時間変動スケールの長い深層循環では予報期間も長い，と素直に理解してよいだろう．
同化計算値と予報結果は，研究所内の www サーバー（http://jes.riam.kyushu-u.ac.jp/）で公開している．本報告で紹介したデータ同化に基づく 2ヶ月予報は，毎週更新している．代表的な海面の水温や流動場の分布をトップページで紹介しているが，任意の時間における水温・塩分・海面高度・流速などの水平・鉛直断面分布を，ユーザーが対話的に選択することもできる．
海況予報は様々な分野に有益な情報となりうる．応用例を挙げれば，原油や汚染物質流出事
故の拡散予測，漁業資源の管理，船舶の省エネ航行，海難救助，陸上気象の長期予報など，暇ない．
例えば，予報された流速値を利用して，エチゼンクラゲの漂流シミュレーションを行った（日
本海区水産研究所と共同研究）．2005 年 8 月中旬を初期値として，1ヵ月後までの漂流計算を行い，9 月 15 日頃にエチゼンクラゲが津軽海峡に達する予測を発表したところ（図 14），実際に

図 14.

日本海データ同化モデルによって予報されたエチゼンクラゲの分布域（2005/8/25
独立行政法人水産総合研究センター報道資料）．

日本海循環データ同化モデル

261

9/11∼12 に下北半島で目撃情報があった．本モデルにより，日本海表層の流動場をかなり正確に再現することに成功し，一方，クラゲ自身の遊泳力は無視できる程度であったと考えられる．
7.

結言

日本海データ同化・予報システムを構築するために，長年九州大学応用力学研究所で培ってきた海洋循環モデルを基に，フェリー観測データを利用して海底地形データを修正し，現実的な計算量で精度を確保できる同化手法を選択し，衛星海面水温と高度計データを同化してきた．その結果，非同化のシミュレーションに倍するデータ検出量を得るなど，一定の推定成果があった．
しかし，本報告で紹介した統計手法は，最先端の統計学や制御工学における精密な解法と比較すると，何世代も前の簡単な手法ばかりである．
「ミニ大洋」日本海とはいえ，単純な制御問
題などと比較して自由度が桁違いに大きいため，適切な推定結果を得るには様々な困難が伴う．
さらに，今回の同化計算によって得た検出量は，観測データの変動強度と比べると半分程度で，なお説明できない（未発見の）変動が残っている．原因としては，観測誤差の影響もあるが，潮汐成分の不正確さや，全域的な気圧変化の影響（Inazu et al., 2006）など，1 日間隔のデータ同化では十分に処理できなかった慣性重力波の影響も考えられる．
大規模で長時間に渡る地球物理現象を前にして，「観測データが不足している」「数値モデルの解像度が低い」「統計推定手法のコストが大きすぎる」等々，言い訳をするのは簡単である．
だが，地球環境の変化は，明日を待たない．地球環境のモニタリングやモデリングの進展と比較して，両者を精度よく統合して地球を知る手段（データ同化）はまだ未熟である，と著者は感じる．推定手法を工夫することによって大きな発見や応用（例えば海底地形の推定，トレーサーの源流追跡，新たなパラメータ依存傾向など）が期待でき，またその積み重ねが我々の明日へと繋がって行くであろう．

参 考 文 献Bretherton, F. P., Davis, R. E. and Fandry, C. B.（1976）. A technique for objective analysis anddesign of oceanographic experiments applied to MODE-73, Deep-Sea Research, 23, 559–582.
Burchard, H., Bolding, K. and Villarreal, M. R.（1999）. GOTM, A general ocean turbulence model,
theory, implementation and test cases, Report EUR 18745, European Commission.
Choi, B. H., Kim, K. O. and Eum, H. M.（2002）. Digital bathymetric and topographic data forneighboring seas of Korea, Journal of Korean Society of Coastal Ocean Engineering, 14, 41–50（in Korean with English abstract）．
Evensen, G.（1994）. Sequential data assimilation with a nonlinear quasi-geostrophic model using
Monte Carlo methods to forecast error statistics, Journal of Geophysical Research, 99, 10143–
10162.
Fu, L.-L., Fukumori, I. and Miller, N.（1993）. Fitting dynamic models to the Geosat sea level observations in the Tropical Pacific Ocean, II, A linear, wind-driven model, Journal of Physical
Oceanography, 23, 2162–2181.
Fukumori, I. and Malanotte-Rizzoli, P.（1995）. An approximate Kalman filter for ocean data assimilation: An example with an idealized Gulf Stream model, Journal of Geophysical Research,
100, 6777–6793.
Fukumori, I., Benveniste, J., Wunsch, C. and Haidvogel, D. B.（1993）. Assimilation of sea surface
topography into an ocean circulation model using a steady-state smoother, Journal of Physical

262

統計数理 第 54 巻

第 2 号 2006

Oceanography, 23, 1831–1855.
Fukumori, I., Raghunath, R., Fu, L.-L. and Chao, Y.（1999）. Assimilation of TOPEX/Poseidon
altimeter data into a global ocean circulation model: How good are the results?, Journal ofGeophysical Research, 104, 25647–25665.
Gordon, N. J., Salmond, D. J. and Smith, A. F. M.（1993）. Novel approach to nonlinear/non-Gaussian
Bayesian state estimation, IEEE Proceedings-F, 140, 107–113.
Gent, P. R. and McWilliams, J. C.（1990）. Isopycnal mixing in ocean circulation models, Journal of
Physical Oceanography, 20, 150–155.
Hirose, N.（2005）. Least-squares estimation of bottom topography using horizontal velocity measurements in the Tsushima/Korea Straits, Journal of Oceanography, 61, 789–794.
Hirose, N., Fukumori, I., Kim, C.-H. and Yoon, J.-H.（2005）. Numerical simulation and satellite
altimeter data assimilation of the Japan Sea circulation, Deep-Sea Research II, 52, 1443–1463.
Inazu, D., Hirose, N., Kizu, S. and Hanawa, K.（2006）. Zonally asymmetric response of the Japan
Sea to synoptic pressure forcing, Journal of Oceanography (accepted).
Ishizaki, H. and Motoi, T.（1999）. Reevaluation of the Takano-Oonishi scheme for momentum advection on bottom relief in ocean models, Journal of Atmospheric and Oceanic Technology,
16, 1994–2010.

海洋情報研究センター（2003）. JTOPO1 北西太平洋 1 分グリッド水深データ, 日本水路協会，CDROM.
Komori, N., Awaji, T. and Ishikawa, Y.（2003）. Short-range forecast experiments of the Kuroshio pathvariabilities south of Japan using TOPEX/Poseidon altimetric data, Journal of GeophysicalResearch, 108, Art. No. 3010.
Lee, H. J., Yoon, J. H., Kawamura, H. and Kang, H.-W.（2003）. Comparison of RIAMOM and MOM
in modeling the East Sea/Japan Sea circulation, Ocean Polar Research, 25, 287–302.
Losch, M. and Wunsch, C.（2003）. Bottom topography as a control variable in an ocean model,
Journal of Atmospheric and Oceanic Technology, 20, 1685–1696.
Manda, A., Hirose, N. and Yanagi, T.（2003）. Application of a nonlinear and non-gaussian sequential
estimation method for an ocean mixed layer model, Engineering Sciences Reports, Kyushu
University, 25, 285–289.
Manda, A., Hirose, N. and Yanagi, T.（2005）. Feasible method for the assimilation of the satellitederived SST with an ocean circulation model, Journal of Atmospheric and Oceanic Technology,
22, 746–756.
Miyazawa, Y. and Yamagata, T.（2003）. The JCOPE ocean forecast system, First Argo Science
Workshop, Tokyo, Japan.
長沼光亮（1977）. 日本海の海況変動，海洋科学, 9, 137–141.
Noh, Y. and Kim, H.-J.（1999）. Simulations of temperature and turbulence structure of the oceanic
boundary layer with the improved near-surface process, Journal of Geophysical Research, 104,
15621–15634.
Sakaida, F., Kudoh, J. and Kawamura, H.（2000）. A-HIGHERS The system to produce the highspatial resolution sea surface temperature maps of the Western North Pacific using the
AVHRR/NOAA, Journal of Oceanography, 56, 707–716.
Sandy, R. J.（1996）. The navy’s bathymetric databases... From the sea, Sea Technology, 37, 53–56.
笹島雄一郎，広瀬直毅，尹 宗煥（2003）. 近似カルマンフィルターを用いた日本海順圧モデルに対する中層フロートのデータ同化，九州大学総合理工学報告，25, 369–378.
．Global sea floor topography from satellite altimetry andSmith, W. H. F. and Sandwell, D. T.（1997）
ship depth soundings, Science, 277, 1956–1962.
水産総合研究センター日本海区水産研究所（2005）. 日本海漁場海況速報，604.
Takikawa, T., Yoon, J.-H. and Cho, K.-D.（2005）. The Tsushima warm current through Tsushima日本海循環データ同化モデル

263

straits estimated from ferryboat ADCP data, Journal of Physical Oceanography, 35, 1154–
1168.
Webb, D. J., de Cuevas, B. A., Richmond, C. S.（1998）. Improved advection schemes for ocean
models, Journal of Atmospheric and Oceanic Technology, 15, 1171–1187.
Yarichin, V. G.（1980）. Study state of the Japan Sea circulation, Problems of Oceanography（ed.
V. Pokudoov），46–61, Hydrometeoizdat, Leningrad（in Russian）.

264

Proceedings of the Institute of Statistical Mathematics Vol. 54, No. 2, 247–264(2006)

Data Assimilation Model for Japan Sea Circulation
Naoki Hirose1 and Atsuyoshi Manda2
1 Research Institute for Applied Mechanics, Kyushu University
2 Graduate School of Science and Technology, Nagasaki University

A data assimilation model has been developed for the Japan Sea circulation as apioneering work for estimation problems in geophysical fluid dynamics. Feasible assimilation methods are selected to reflect various measurement data in a high-resolution ocean
circulation model considering computational cost and accuracy. Realistic surface and subsurface states are estimated and predicted by controlling not only the surface boundary
conditions but also the bottom topography in the Japan Sea. Ocean current forecasting
benefits many scientific and social applications including prediction of drifting oil spills
and movement of giant jellyfish.

Key words: Japan Sea, ocean circulation model, remote-sensing observation data, data assimilation,
ocean forecasting.