土木学会論文集B2（海岸工学），Vol. 73, No. 2, I_1663─I_1668, 2017.

3次元変分法を用いた瀬戸内海流動
再解析・予報モデルの高精度化
黒澤賢太1・内山雄介2・三好建正3,4,5

1学生会員

神戸大学大学院 工学研究科市民工学専攻（〒657-8501 神戸市灘区六甲台町1-1）
E-mail:vppbjb@gmail.com
2正会員 神戸大学教授 大学院工学研究科市民工学専攻（〒657-8501 神戸市灘区六甲台町1-1）
3（国研）理化学研究所 計算科学研究機構（〒650-0047 神戸市中央区港島南町 7-1-26）
4Dept. Atmos. & Ocean. Sci., University of Maryland, College Park（College Park, MD 20742-2425, USA）
5（国研）海洋研究開発機構 アプリケーションラボ（〒236-0001 横浜市金沢区昭和町3173-25）

瀬戸内海全域を対象とした流動・水質予測精度向上を目的として，領域海洋循環モデルROMSとROMSに適合させた3次元変分法に基づくデータ同化（DA）エンジンを連成し，瀬戸内海高解像度ROMS-3DVARシステムを開発した．通常の海洋DAにおいて核心的に重要なデータである衛星海面高度やArgoデータは内湾である瀬戸内海では使えないため，試行的に船舶による水温塩分の定期観測データを同化し，観測システムシミュレーション実験（OSSE）によって性能評価を行った．3DVARの導入によって水温・塩分分布の再現性が飛躍的に向上することを定量的に示すとともに，特に塩分場の再現性向上のためには単にデータ同化だけに頼るのではなく，精緻なforwardモデルと正確な流入淡水情報が欠かせないことを示した．

Key Words : data assimilation, 3DVAR, ROMS, OSSE, NMC Method, Seto Inland Sea

1. はじめに得されている．本研究で対象とする瀬戸内海においても，国土交通省近畿地方整備局が海洋環境整備事業の一環と沿岸域における正確な海況予報は，港湾管理・海洋建して行なっている観測データ，大阪湾水質定点自動観測

設・漂流物回収・漁場予測・航路選択・海洋環境モニタシステムにより取得される多点の流速・水温・塩分の鉛リングなど，様々な分野から期待されている．しかしな

直分布データ，海洋短波レーダーによる表面流速等が存

がら，瀬戸内海のような内湾における流動や水質は，局在する．大阪湾における試行的なDA実験によって，表地気象，河川からの淡水流入，黒潮流路変動に代表され層流や密度場の推定精度が向上する可能性も指摘されてる外洋流動などの影響を強く受けて複雑に形成されるたいる3),4)．

め，ごく単純な潮流などを除くと十分な再現性が得られないことが少なくない．正確な海況予報を行うための手DAを行う際には，適用する解析領域や観測データ，使用する計算資源などを考慮し，適切な手法を選択する法として，観測データと予報モデルを融合して海洋の尤必要がある．アンサンブルカルマンフィルタや4次元変もらしい状態を推定するデータ同化（以下DAと呼称）

分法のような高度なDAでは，誤差の時間発展を考慮す

1)は非常に有効である ．DAの導入に伴い，これまでの海るために数値予報モデルを繰り返し実行する必要があり，洋モデルでは予測できなかった海洋状態の推定や重大な

その計算量は極めて膨大である．簡易な手法であるナッ発見が多くなされ，工学や科学の発展にも大きく貢献しジング法や最適内挿法では，力学的なバランスを考慮せ2)てきた ．一般的な海洋DAでは，衛星海面高度・表層水ずにモデル変数を機械的に観測値に近づけようとするた温，Argoフロート，各種現地観測データ等が使用される

め，支配方程式が満たされなくなって非現実的な循環流が，半閉鎖性内湾域や沿岸域では，陸面ノイズに伴う衛が生成されるなど，かえって再現性を損なう場合がある．

星データの精度劣化やArgoの欠測が深刻なボトルネック

変分法にもとづくDAでは，モデル予報値と観測値の間になり，外洋域や陸棚海域と比較してDAの導入が立ちのバランスを統計的に評価し，最適な解析値を求めるこ遅れていた．とがエッセンスである．このとき，予報値と観測値との一方，特に我が国の三大湾においては，海上交通，水差を評価関数として繰り返し計算により最適な解析値を産，環境アセスメント等のために様々な観測データが取求めることになる．3次元変分法（以下3DVAR）では海I_1663瀬戸内海東部流動モデル6), 7)を用いる（図-1）．ROMS-L1
モデル領域は東西640 km×南北640 km（格子数320 × 320 ×鉛直32層），L2領域は，東西480 km × 南北288 km（同800× 480 × 鉛直32層）である．領域内に流入する一級河川（L1：29本，L2：27本）を考慮し，流量には日本河川協会の雨量・流量データベースから求めた10年間の月平均気候値を与えた．その他流動モデルの計算条件については表-1に示すとおりである．なお，DA計算の際に用いる流動モデル再解析値は，特記しない限りはROMS-L2モデルの結果を使用している．
図-1 2 段ネスト高解像度瀬戸内海流動モデル（L2）領域．外
枠は親グリッド（L1）領域，カラーは水深 (m)．

3. 3次元変分データ同化
本研究では，三好らによって開発された全球気象モデ
ルSPEEDY-3DVAR8)を独自にROMSに実装した．以下で

表-1 L2 領域計算条件は，実装した3DVARの概要について説明する．

計算期間

2007 年 1 月 1 日〜2015 年 5 月 1 日

格子数

802 × 482 × 鉛直 32 層

水平空間解像度600 m

同化手法である．評価関数の最小値を繰り返し計算によ

初期・開境界条件

ROMS−L1（水平解像度 2 km）って数値的に探索し，得られた最小値を解析値とする．

外力（海上風）海面フラックスJMA GPV-MSM（1 時間値）
NOAA-COADS（月平均気候値）

つまり3DVARでは，以下のように定義される評価関数

潮汐

TPXO7.0 全球調和定数（10 分潮）

一級河川流量日本河川協会雨量・流量データベ ース（27 本，月平均気候値）

洋状態の時間変化は考慮せず，様々な観測が解析時刻に

1𝐱 − 𝐱 ( ) 𝐁 +, 𝐱 − 𝐱 (21+ 𝐇𝐱 − 𝐲0 ) 𝐑+, 𝐇𝐱 − 𝐲0(1)2
ここに，x：解析対象の状態変数，xb：状態変数の第一得られたと仮定するため，同化の際に数値予報モデルを推定値（モデル予報値），B：背景誤差共分散行列，実行せず，4次元変分法と比べると計算量が圧倒的に少なくて済むという大きな利点がある．そのため，現在運H：モデル空間から観測空間への変換行列である観測演
算子，yo：観測値からなるベクトル，R：観測誤差共分

用されている海洋再解析・予報システムでは，3DVAR散行列である．Hxはモデル空間から観測空間に変換さによるDAが実装されている1),2)．

れた状態変数を表す．BとRの対角成分は誤差分散に等DAの一つである3DVARは，最尤推定法を基礎とするJ の値を最小にする x の値を求めることになる．
𝐽 𝐱 =以上の背景を踏まえ，本研究では，瀬戸内海におけるしく，誤差共分散を表す非対角成分は誤差相関に比例す

内部流動や密度場の予報精度向上を目的として，領域海る．すなわち，Bは第一推定値の制御変数に含まれる誤5)

洋循環モデルROMS に適合する3DVARベースのDA技術差，つまりモデル内での予測のばらつきの範囲を示す．

開発を行う．また，このDAシステムを用いて，現地観また，Rの非対角成分は，観測誤差の空間的な相関を表測ネットワーク配置や，観測の時間分解能が予報精度へす．本研究では，背景誤差共分散行列Bの推定にはNMC与える影響を定量的に評価し，高精度予報を行うための法9)を用いる．NMC法とは過去のある期間の予報モデル

理想的な観測システムについての検討を行う．さらに，の誤差を蓄積し，初期時刻の異なる2つの予報について実際に瀬戸内海で行われた現地観測による水温塩分デー

同じ対象時刻での予報値の差を予報誤差とみなす方法でタを同化し，予報精度の向上を定量的に評価する．

ある．つまり，モデル予報値の誤差の情報を含む多数の事例を統計サンプルとみなすため，推定されるBは使用

2.

数値予報モデルの概要

する数値予報モデルにおける背景誤差の特性を表すこと

本研究では，水平解像度約10 kmのJCOPE21) 海洋再解
析値を最外側境界条件および初期条件として与え，領域海洋循環モデルROMSを用いた2段階の1-way offlineネス
ティングにより，ROMS−L1（水平解像度2 km）→L2（同600 m）へと順次ダウンスケーリングした高解像度になる．本研究では，2013年5月1日~5月31日の31日間
62例について，24時間予報値と48時間予報値の予報誤差を推定し，その共分散行列を求めることにより，背景誤差共分散行列を統計的に推定した．
図-2(a)に模式的に示すように，実際にROMSを数値予報モデルとする3DVAR計算では，以下の①から④の手

I_1664

図-2 ROMS を数値予報モデルとした 3DVAR データ同化サイクルの模式図．(a) 通常の同化計算の流れ，(b) OSSE の同化計算の流れを示す．「DA」はデータ同化システムを指す．

順を反復実行する．① 初期条件，境界条件，モデルパ

手順⑥において観測演算子（例えば疎らな観測網を再ラメーターなどを設定し，ROMSを用いて予報を行う．

現するためのsubsamplingなど）をNRに作用させ，観測値② 計算された予報値をデータ同化を行うための第一推に相当する物理量を算出し，その値に観測誤差を付加す

定値とみなし，第一推定値と観測データを3DVARで同化ることで，NRと整合する観測データを疑似的に作成すする．③ 同化後の結果として解析値を出力する．④ ③

る．つまりNR-OSSEでは，NRが真値である仮想世界をで計算された解析値を次の同化サイクルの初期値として構築して同化精度を評価するための手法であると解釈さ①に戻り，次の同化サイクルの計算を行う．

れる．真値が既知であるため，同化に必要な統計パラメーターの設定や同化結果の精度の評価が現実の問題を取4.

観測システムシミュレーション実験（OSSE）り扱う場合に比べてはるかに容易に行える．
本研究では，DAに伴う予報値と真値のずれをRMSE(1) OSSEの概要
本研究で開発したROMS-3DVARシステムを瀬戸内海（自乗平均誤差平方根）によって評価する．
𝑅𝑀𝑆𝐸 =海況予測へ適用し，その精度向上に対する有効性を評価するとともに，瀬戸内海における現地観測点の配置が予報精度に及ぼす影響について検討する．ここでは，観測システムシミュレーション実験（Observing System Simulation Experiment．以下OSSEと呼称）を行い，定量的な評価を行う．OSSEとは，実在しない観測システム（仮想観測システム）を計算機上に構築し，DAの振る舞いを評価する手法である．DAによる解析値を観測値と直接=(89 +:9 (;))

<>

(2)

ただし，N：比較対象とするデータ数，x：モデル予報値，yi：i地点における比較対象値（NR-OSSEではNRによる仮想真値，あるいは観測値），Hi：観測演算子である．
したがって，Hi (x)はi地点における観測空間に内挿されたモデル予報値となる．
(2)

NR-OSSEによる評価結果

比較する場合に比べ，OSSEでは観測値が含む誤差を一瀬戸内海内部における現地観測ネットワーク構造が切排除してDAの効果を計測することも可能である．あるいは，例えば将来的に配置する予定である観測システ

ROMS-3DVARシステムの解析精度に及ぼす影響を評価するために，図-3に示すように疑似観測網を瀬戸内海内ムが数値予報に与える影響や，その費用対効果を調べ，部に設定した．図-3(a)ではROMS-L2の全格子点に対し，観測システムの評価を行う際にも使われる．また，既存の観測システムにおいて，観測回数や観測地点を現状よ陸地を避けて水平方向に10グリッドごとに1点ずつ抽出した計582地点，図-3(b)では14グリッドごとの計311地点，りも増減させ，その観測データを同化した際の予報精度

図-3(c)では30グリッドごとの計69地点の観測網が仮想的への影響を評価することで，理想的な観測システムにつに存在することを仮定している．また，図-3(d)に関していての考察や提言を行う場合などにも有用である．

は，2013年の瀬戸内総合水質調査（詳細は次章）による

ここでは，NR-OSSE（図-2(b)）と呼ばれる手法を用いて開発したROMS-3DVARの効果を計測する．NR-OSSE全観測地点（229点）に対応している．

では，図-2(b)の手順⑤に示すように，3DVARサイクルの①で設定する計算条件とは異なる設定でモデル計算

OSSE解析期間は2013年1月1日から2013年12月31日までの1年間とし，表-2に示した5ケースのDA計算と，DAを行う．この結果は仮想的な真の海洋状態，ネイチャーを行わないケースの計6通りの計算を行う．ケース1〜5はそれぞれ図-3(a)〜(d)の擬似観測点で同化させ，ケーラン（Nature Run．以下NRと呼称）と呼ばれる．次に，ス4以外は，表層から鉛直下方向10層分までの鉛直分布

I_1665図-3

OSSE 解析で同化に用いた疑似観測点の水平分布．(a) L2 を 10 点ごとにサブサンプルしたもの（ケース 1），(b) 14 点ごと（ケース 2），(c) 30 点ごと（ケース 3），(d) 2013 年瀬戸内総合水質調査における観測点（ケース 4，5）

図-4 OSSE 解析結果．予報値の真値に対する RMSE の 1 日ごとの経時変化．黒線は同化なし計算を示す．

図-4は予報値の真値に対する塩分と水温のRMSEの経表 -2 OSSE 解析で行った DA あり全 5 ケースの諸条件．

時変化を示したものである．いずれのケースもDAなし

観測点の観測点の同化の場合（黒線）よりもRMSEが低下しており，DAによっ水平分布鉛直分布サイクルて予報精度が向上したことが分かる．RMSEは観測点のCASE 1(a)表層から 10 層

1 日ごと数に反比例して低下しており，最も多くの点で同化した

CASE 2(b)表層から 10 層

1 日ごと

ケース 1（マゼンダ）が最も誤差が小さい．実際の現地

CASE 3(c)表層から 10 層

1 日ごと

CASE 4(d)(d)

観測に準ずる

観測に準ずる表層から 10 層

1 日ごと

CASE 5

観測は時空間的に著しく疎らであるため，ケース4は原くく，精度向上幅も小さいが，DAによるRMSE低下は明白である．また，ケース4と5を比較すると，現行の観を同化に用いるとともに，同化サイクルは1日とする．
ケース 4 は実際の運用を想定したDAシステムの精度評価を目指しており，実際に行われた現地観測と同一座標，同一時刻（日単位）において随時同化を行う．

理的に他のDAありのケースよりもDAの効果が得られに測を使用する場合においても，観測頻度を上げることで顕著な予報精度向上が見込まれることが分かる．以上の結果から，開発した瀬戸内海ROMS-3DVARシステムはDAによって再現精度が向上しており，期待通りにDAが機能していることが確認された．

I_1666

図-6 同化したすべての観測値とそれに対応する予報値の散布図 (左: 水温, 右: 塩分)．青点, 赤点はそれぞれ同
化なしおよび同化あり計算の予報値，R は相関係数を示す．
図-5 瀬戸内総合水質調査（2013 年実施分）による月ごとの観測地点数．

5.

同化実験
本章では，実際に瀬戸内海で計測された現地観測デー

タを同化することにより，瀬戸内海モデルの予報精度が向上することを示す．同化実験に用いるデータは，
図-7 予報値の観測値に対する RMSE の 1 ヶ月ごとの計時変化 (左: 水温, 右: 塩分)．青線，赤線はそれぞれ同化なし，同化あり計算を示す．

OSSE解析にも用いた瀬戸内総合水質調査観測データセット（http://www.pa.cgr.mlit.go.jp/chiki/suishitu/）である．
このデータは，国土交通省による瀬戸内海における海洋環境整備事業の一環として，海面清掃船によって定期的に観測された各種データを取りまとめたものである．ここでは，2013年にCTDで観測された水面下2 mの水温・塩分データを同化対象とし，同化サイクルは1日とする（観測時間の日内のばらつきは無視）．2013年の月ごとの観測地点数は図-5に示すとおりである．
同化したすべての観測値と，それに対応する同化あり（ROMS-3DVAR）および同化なし（ROMS）の予報値の散布図を図-6に示す．水温に関しては，同化なし計算であっても観測値との相関は相関係数0.987と高く，そもそも同化前の瀬戸内海ROMSモデルの再現精度が極めて高く，同化によってさらに精度が向上することが分かる．
一方，塩分の相関係数は同化前は0.533，同化後は0.789であり，著しく精度向上したものの，同化後の水温の相関係数0.997には及ばない．そこで次に，同化対象となる観測値が存在するすべての地点における予報値と観測値を用いて1ヶ月ごとのRMSEを求めることで，再現性能の経時的な評価を行なった（図-7）．同化サイクルが進むにつれて水温，塩分のRMSEはともに減少していく様子が確認できるが，塩分に関しては水温ほどのRMSEの改善は見られない．特に同化なしの場合，9月と10月の塩分推定精度が低く，この時期に顕著となる降雨出水に伴う河川からの淡水影響が十分に反映できていないことが主たる原因であることが強く示唆される．本瀬戸内海モデルでは，陸水影響としては一級河川のみを考慮し，月平均気候値流量を淡水として流入させているが，出水イベントや，データ整備が立ち遅れている二級河川・地下水などからの淡水供給の影響が正確に考慮されていないことが最も大きな原因であると考えられる．また，本図-8 表層水温分布（℃）．（左）衛星データ（MODIS/Aqua）2013 年 10 月 10 日午前 3 時 48 分の瞬間値，（中）ROMS-3DVAR（同化あり），（右）ROMS（同化なし）計算結果．計算値は 2013 年 10 月 10 日午前 4 時の瞬間値である．

I_1667

研究では境界条件の信頼性は非常に高いものとして観測

謝辞：本研究は科研費補助金（15H04049，15H00977），

値との比較を行なっているが，気象による影響や境界条同国際共同研究加速基金（15KK0207）の援助を受けた．

件に含まれる不確実性に起因する誤差も少なからず含まれているとも考えられる．

参考文献

最後に，瀬戸内総合水質調査観測データを同化したこ

1)Usui, N. and coauthors: Meteorological Research Institutemultivariate ocean variational estimation (MOVE) system :Some early results, Advances in Space Research, 37, 806822, 2006.

2)

Miyazawa, Y. and coauthors: Water mass variability in theWestern North Pacific detected in 15-year eddy resolving
ocean reanalysis, J. Oceanogr., Vol. 65, pp. 737-756, 2009.

3)

入江政安，藤井裕司，岡田輝久，西田修三：海洋レーダーによって観測される表面流速の 4 次元変分法を用いたデータ同化，土木学会論文集 B2（海岸工学）,Vol. 69, No.2, pp.I_501−I_505, 2013.

4)

入江政安，西田修三，辻陽平，岡田輝久：4 次元変分データ同化法を用いた密度分布と水質の再現性向上，
土 木 学 会 論 文 集 B2 （ 海 岸 工 学 ） , Vol. 68, No.2,pp.I_446−I_450, 2012.

5)

Shchepetkin, A. F. and McWilliams, J. C.: The Regional
Oceanic Modeling System: A split-explicit, free-surface,
topography-following- coordinate oceanic model, Ocean
Modell., Vol.9, pp. 347-404, 2005.
内山雄介，小硲大地，西井達也，御手洗哲司：コネクティビティを用いた瀬戸内海全域における海洋生態系ネットワーク構造の解析，土木学会論文集 B2（海 岸工学）, Vol. 70, No. 2, pp. I_1076−I_1080, 2014.
小硲大地，内山雄介，御手洗哲司：黒潮流路変動が瀬戸内海全域の幼稚仔分散過程に及ぼす影響，土木学会論文集 B2（海岸工学）, Vol. 71, No. 2, pp. I_1231− I_1236, 2015.
Miyoshi,T. : Ensemble Kalman Filter Experiments with a
Primitive-equation global model, Ph. D. dissertation, University of Maryland, College Park, MD, 2005.
Parrish, D. F. and Derber, J. C. : The National Meteorological Center’s spectral statistical interpolation analysis system, Mon. Weather Rev., Vol. 120, pp. 1747–1763, 1992.
(2017. 3. 15 受付)とによる再現性の向上を別の角度から確認するために，3DVARに用いていない観測データとの比較を行なう．ここでは，比較対象としてMODIS/Aquaによる衛星海表面水温（SST）データを用いる．図-8は衛星データ，同化あり，同化なし計算によるSST分布を示している．DAによって衛星SSTに見られる水温分布が全体的により精緻に再現されていることが見て取れる．特に，播磨灘，広島湾，周防灘，豊後水道付近においては，明らかな再現性の向上が確認される．

6. おわりに本研究では，3次元変分法（3DVAR）を2段ネスト瀬戸内海流動モデルのROMS-L2モデルに組み込み，瀬戸内総合水質調査による水温・塩分データを同化させた
ROMS-3DVARシステムを開発した．OSSE等によってデータ同化（DA）による精度向上を定量的に評価し，瀬6)戸内海内部における水温・塩分分布の再現性がDA に よって飛躍的に向上することを明確に示した．また，瀬戸内海塩分場の再現性をより向上させるためには，正確な

7)流入淡水情報（二級河川や出水イベントを含む詳細な流量データ）が欠かせないことが示唆された．
本 ROMS-3DVAR システムの成功は，第 5 章で検討し8)たように極めて高精度な DA なしモデル（forward モデル）が下敷きになっている．つまり，高精度な海洋予報のためには単に DA 技術に頼るのではなく，使用する海洋モ

9)デルや海域の特性を深く理解し，外力・境界条件設定を含めた総合的なモデリング技術の質的向上を同時に進めていくことが本質的に重要であると考えている．

ON IMPROVEMENT OF AN ESTUARINE REANALYSIS-FORECAST MODEL
FOR THE SETO INLAND SEA BASED ON 3D VARIATIONAL ASSIMILATION
Kenta KUROSAWA, Yusuke UCHIYAMA and Takemasa MIYOSHI
Ocean data assimilation (DA) has widely been used for various engineering purposes including monitoring, forecast, and port operations. However, in semi-enclosed estuaries such as the Seto Inland Sea(SIS), Japan, DA has rarely been adapted because satellite data and Argo floats are desperatively unavailable, albeit being essential to successful oceanic predictions. In the present study, we develop a DA system toward rigorous predictions optimal for relatively small coastal and estuarine settings. The system
consists of a three-dimensional variational (3DVAR) assimilation scheme being made compatible withROMS using an in-situ dataset available for the SIS. The DA system significantly improves the reproducibility of water temperature and salinity fields from those without DA. We show an evidence that detailed
terrestrial freshwater discharge data is inevitable for more accurate salinity predictions in the estuary.

I_1668