海洋物理におけるデータ同化システムの発展と今後の展望

石川洋一 (海洋研究開発機構・地球情報研究センター)

はじめに• 「データ同化」は地球科学において、一般的なものとなってきている• さらに、データ同化はこれまで使われてこなかった分野にも広がっており、ますます重要な技術となっている• 惑星科学分野では観測データの量が問題となる▫ 海洋分野では気象分野に比べ、データが少ない中でデータ同化研究が進められてきた

データ同化とは• 観測データと数値モデルを組み合わせてデータセットを作成すること

観測データからみると…
数値モデルを利用した4次元的な内挿／外挿／グリッド化
異なる種類の観測データの統合

数値モデルからみると…
モデルの入力パラメータ(初期条件、境界条件など)のチューニングの自動化

海洋観測データ• 人工衛星リモートセンシングデータ▫ 海面水温、海面高度、海色など▫ 面的にカバーされ、最も多いが海面の情報のみ• 現場観測データ▫ 水温、塩分プロファイル▫ 船舶からの観測▫ ARGOフロート:自動的に昇降して水温、塩分プロファイルを観測し、衛星通信によってデータを送る

海洋観測データ• 人工衛星リモートセンシングデータ▫ 海面水温、海面高度、海色など▫ 面的にカバーされ、最も多いが海面の情報のみ• 現場観測データ▫ 水温、塩分プロファイル▫ 船舶からの観測▫ ARGOフロート:自動的に昇降して水温、塩分プロファイルを観測し、衛星通信によってデータを送る

観測データ数はどれくらいあるか“同化”「大辞林 第二版」
1. 本来異なる性質や考え方が同じものになること。
2. 外から得た知識などを理解して自分のものとすること。
3. 生物体が外界から摂取した物質に特定の化学変化を加え、その生物に固有あるいは必要な物質を作り出すこと。同化作用。アナボリズム。
4. マグマが周囲の岩石や外来物質を取り込み、一つのものに混合すること。同化作用。
5. ある音素が隣接する音素に影響されてそれと同じ、または似た性質のものに変化すること。エビス（ebisu）の i が先行する e と等しくなってエベス（ebesu）となる類。

データ同化とは？
• 観測データと数値モデルを組み合わせて、現実的なデータセットを作るための手法▫ 複数の観測データを統合した時系列データセットを作成する▫ 数値モデルで用いているパラメータなどの最適化を観測データをもとにして行う▫ 数値予報のための初期条件、境界条件を作成する▫ 観測システムの設計のための評価、改善にむけた情報提供

データ同化結果の例

７月

同化結果

(亜表層水温場)

東北海区水温図(東北水研より)

100m

200m

データ同化と関連した数学• 微分方程式• 行列、線形代数• 確率論(確率微分方程式)• 統計学(誤差論)• 最適化理論
初期値、パラメータなど• 最適制御
数値モデル
• 逆問題

データ同化

時系列データ、数値予報

観測データ
比較

データ同化とは？
• 観測データと数値モデルを組み合わせて、現実的なデータセットを作るための手法▫ 複数の観測データを統合した時系列データセットを作成する▫ 数値モデルで用いているパラメータなどの最適化を観測データをもとにして行う▫ 数値予報のための初期条件、境界条件を作成する▫ 観測システムの設計のための評価、改善にむけた
情報提供

数値モデルに観測データを“同化”する。

数値予報とデータ同化
予報サイクル

修正
予報

予報

予報

修正
データ同化
データ同化: 観測データを用いた予報値の修正
予報値が得られる
観測値が得られる
両者を組み合わせて解析値を作成する

データ同化による初期値化データ同化手法の基礎
線形不偏最適推定理論
Best Linear Unbiased Estimation (BLUE)
ある変数x (未知の真値:xtrue)に対して異なる2つの推定値x1, x2およびそれぞれの誤差σ1,σ2が与えられたとき、xに対する最適推定値は次式で与えられる。
ただし、x1、x2の期待値はxtrueと等しいとする。

不偏最適推定の簡単な例題
ある変数ｘ(未知の真の値=0)についての推定値としてx1:期待値 0, 標準偏差1
x2:期待値 0, 標準偏差2があたえらた時の、最適推定値を求める。

不変最適推定理論から最適解はで与えられる。

不偏最適推定理論における注意点
求められた最適解は統計的(確率的)に最適であり、個々のケースでは必ずしも真の値に近づくとは限らない。
最適解の精度はあらかじめ推定された誤差の精度に依存している。実際に推定誤差の比が重要である。
現実には、バイアスがない、誤差に相関がないという仮定が成り立つ保証がない。
実際のデータ同化手法はこの不偏最適推定理論を多変量、多次元に拡張したものである。

データ同化手法の分類
Sequential法

ここで修正最適内挿法
カルマンフィルタ
3D-VARなど

前回の同化時刻

Iterative法
4D-VAR
RTSスムーザなど

今回の同化時刻

前回の予報=今回の初期推定同化結果

期間で修正カルマンフィルター(x = x + K y - Hxaff)(

K = PH HPH + RTT)

-1

行列Hは観測行列と呼ばれ、モデルの予報値を観測点に射影するオペレータである。
カルマンフィルタは予報誤差共分散行列の時間発展を予報モデルを用いて求めるものである。
予報誤差共分散行列Pはあらかじめ与える必要があり、最適内挿法などの簡単な同化手法では、空間方向にガウス関数を仮定する。

flow dependent covariance
流れや波によって情報が運ばれるので、データ同化による修正は非等方的であるべき

Kuragano & Kamachi 2000

カルマンフィルター(x = x + K y - Hxaff)(K = PH HPH + RTT)

-1

カルマンフィルタは予報誤差共分散行列の時間発展を予報モデルを用いて求めるものである。
誤差の空間スケール、変数間の関係がモデルの方程式に対応して求められる
流速場などの変動にあわせて時々刻々誤差共分散行列が得られる

アンサンブルカルマンフィルター(x = x + K y - Hxaff)(

K = PH HPH + RTT)

-1

カルマンフィルタは予報誤差共分散行列の時間発展をリカッチ方程式を用いて求める。
アンサンブルカルマンフィルターでは予報誤差共分散を予報モデルの多数回試行から求める。
リカッチ方程式における線形化の問題を解決
少ないアンサンブル数で誤差共分散行列を求める手法を開発することにより、計算負荷の軽減

3D-VAR(3次元変分法)(J−= x - xaf) (T-1aB x -xf)() (T+ Hx - y R -1 Hx a - ya)

誤差共分散行列が同じ場合、最適内挿法と同じ最適解が得られる。
Ｈとして非線形の観測オペレータを用いることができるため、観測データの有効利用ができる。
⇒リモートセンシングデータのリトリーバル
容易に付加的な制約条件を課すことができる。
例えば、空間的な滑らかさ、力学的バランスなど

4次元変分法(adjoint法)
J = (c - c ) B (c - c )+ (Hx - y ) R -1 (Hx - y )af T

L=J+λT・Ax-1afT(制御変数cは初期・境界条件、パラメータなど)

数値モデルの前方積分とadjoint方程式(λn=A*λn+1)の後方積分を繰り返して、制御変数の最適解を求める。
数値モデルと随伴方程式

線形システムにおいてadjoint法の最適解は•同化期間の最後でカルマンフィルターと一致し•期間にわたってRTSスムーザーと一致する
非線形システムでは線形化の方法に依存
アンサンブルカルマンフィルターなどと一致する保証はないadjointモデルによる後方積分データ同化システムの設計• 観測データ、数値モデル、同化手法のバランスが重要• 観測データが豊富にあれば、簡単な同化手法でも良い推定値が得られる• 数値モデルが複雑になると、表現される現象が多様になるので、制御するのが難しい▫ どのような変数／パラメータを修正して、どのような現象を再現するか

結合データ同化システムの概要• Assimilation with a fully coupled GCM• Assimilation method is 4D-VAR▫ Long assimilation window (9 month)▫ Correction of model climatology by parameter
estimation▫ Correction of seasonal to interannual trajectory byinitialization▫ Atmospheric data are also assimilated▫ Weather mode is treated as noise季節から数年の時間スケールをもつ気候変動現象をターゲットにする。Weatherモードの変動についてはモデルの挙動にゆだねるが現象によっては平均場が改善されることで、implicitに物理場が修正されることが期待される。

結合同化システムの3つのタイプ:(Conventional independent/forced A & O, Quasi CDAS, & CDAS)(a) Independent Systems (Conventional DAS)(b) Quasi CDAS（MOVE-C）
Atmos. Obs.

Atm. Obs.

AGCM
大気モデル
結合

Atmos. Forcing

予測△ Prediction
Coupled

Coupling

AGCM
大気モデル

Coupling

SST

CoupledPrediction
予測
OGCM
海洋モデル

OGCM
海洋モデル
Ocean Obs.

Initial Shock

Ocean Obs.

3DVAR(c) Fully coupled DAS (CDAS: K7）
Atmos. Obs.
大気観測データ
Coupling
結合AGCM
大気モデル
OGCM
海洋モデル

CoupledPrediction
予測？
esp. for seasonal to interannual variation
Different Scales

Ocean Obs.
海洋観測データ->filtering consistently

1. Straight-forward concept
2. No initial shock
3. Higher forecast potential

結合同化システムの概要•••

Coupled Model (CFES):▫ T42L24 AFES for AGCM
▫ 1x1deg L45 MOM3 for OGCM
▫ IARC SeaIce model▫ MATSIRO Model for Land
Observational Data
▫ Atmosphere:
 NCEP’s BUFR data U,V,T,Q (10daily) SSM/I sea wind scalar x ERA40 wind direction (10daily)▫ Ocean: T/P altimeter data(10daily)
 Reynolds SST (10daily)
 WOD data T,S (monthly ) Ocean Data Assimilation Product T,S(monthly)
Adjoint Code▫ Adjoint OGCM and adjoint AGCM are coupled [Line by linetransformation by TAMC,TAF]▫ Temporal averaging of forward field for the adjoint integration is applied
to smooth the basic field▫ Adjoint AGCM contains damping terms to suppress the strong adjoint
sensitivity from weather fluctuations.

(∂λ  ∂ M T−1− = λ − Γλ + H R H x − y∂t  ∂x  x = xT)

λ : adjoint variables, x : temporal average, − Γλ : damping.

気候変動現象indexの再現性• CDA Reanalysis experiment for the period of
1980-2000

４D-VAR同化手法を用いたENSO予報研究(初期値：Case of 3-month assimilation window;
1997/98)
NINO3.4 Anomaly Correlation Score
Prediction by CDA

Prediction by ODA &IAU

Persistent

imply
capability of
1.5-year lead
prediction ofENSO

結合同化システムによるインドモンスーン予測• インド熱帯気象研究所（IITM）との協働により、気候災害の軽減と水産資源管理に向けたインド洋域での海洋・大気間の熱・水リンクのルート解明及びインドモンスーン予測研究を開始した。今後も継続的に協働研究を進めていく。
2000年5月

SeaWiFS観測によるchl-a

2000年9月

Dynamical monsoon index

DMI=U*850-U*200
結合同化システムによる
インドモンスーン循環の再現Correlation coefficients（vs. CDA）
ERA40 → 0.720
JRA25 → 0.765
NCEP1 → 0.674
NCEP2 → 0.676

CDA forecast for 2010

Assimilation window
2010Jan.3 month
2010Mar. 2010Apr.

同化データ (10-day mean)
・NOAA OISSTSSTアノマ・SSM/I NINO3
wind speed* JRA25 wind direction
リT,Q,U,V・NCEP PREPBUFR

forecast

2010Dec.

初期条件
・50年海洋長期再解析から2010年に似たケース(Jan.1995) を使用

インドモンスーン指数

結合4DVARによる4月初旬の初期値化予測は、
６・７・８月のモンスーンの観測結果を7割程度予測できた。

Parameter Estimation using Green’s Function method

Toyoda et al. (in revision, IMBER/IMBIZO special issue)
We divide the model domain (75oS-80oN) into 40 provinces based onLonghurst (1995).
Coastal provinces in Longhurst (1995) are excluded regarding the coarse
resolution of the ocean reanalysis.
Japanese Sea, Okhotsk Sea, Costa Rica Dome region, and region off thecoast of Peru are added taking into consideration of the specificcharacteristics of the physical environment.
In each province, the optimal set of the NEMURO parameters areestimated. Note that the region within 2 degree from the coast line (shaded)is eliminated from the optimization.

Optimization

We select 32 independent control variables from all 70 NEMURO parameters.
those which differ largely between sites (Fujii et al., 2007; Yoshie et al., 2007)
those which have large sensitivity (Yoshie et al., 2007)
coefficients of photosynthetic rate, respiration, and mortality for each
phytoplankton and of grazing and mortality for each zooplankton are estimateddependently.
Cost function to be minimized is defined as follows:interannual trend during the experiment (Jtr)differences between the analyzed and the observations of nutrients (Jnt)same as Jnt but of surface chlorophyll concentration (Jcl)difference between the analyzed control variables and their first guesses (Jbg)

Examples of improved time series of surface chlorophyll concentration:(left) central ant (right) eastern tropical Pacific (180oE-160oW, 5oS-5oN; 90-100oW,
20oS-5oN)
blue: control experiment, red: optimized experiment, black: SeaWIFS data
•In these regions (ones of the most improved regions), the mean and interannualvariability of the surface chlorophyll concentration are quantitatively reproduced in theoptimized experiment.
•By analyzing the NEMURO variables and the physical fields in the 4DVAR reanalysis,
the interannual variability in phytoplankton in these regions is generated by thevariability in nutrients, which is attributed mainly by the zonal advection change, not by
the vertical one.

1997

2006

1997

2006

海洋５０年再解析を用いた北太平洋最深層部の昇温メカニズム解明• 4D-VAR海洋長期再解析データセット(1957-2006年)を用いると共に、アジョイントモデルによる感度実験を行うことにより、北太平洋最深層部の昇温メカニズムの解明研究を行いScienceに掲載された。
• 北太平洋（北緯47度）の底層でJAMSTECが発見した水温上昇は、時間を遡って解析するとおよそ４０年前の南極アデリーコースト沖周辺での海面フラックス変動の影響であることがわかった。

47Nに沿った船舶観測による1999年と1985年の水温差(Fukasawa et al.,2004)

海洋５０年再解析データセットを基にした４D-VAR全球全層海洋同化システムによる感度

How about the Physical Mechanism?

?

Difference of the heat storage between
WOCE-WOCE revisit observationalperiods.

The physical mechanisms governing bottomwater warming are poorly understood since insitu observations are spatially and temporallysporadic.
The changes in heat storage between WOCEWOCE revisit imply northward running of thewarming signal, but…

Validation with reanalysis dataset;Term balance of the temperature equation

海洋５０年再解析を用いた北太平洋最深層部の昇温メカニズム解明• 4D-VAR海洋長期再解析データセット(1957-2006年)を用いると共に、アジョイントモデルによる感度実験を行うことにより、北太平洋最深層部の昇温メカニズムの解明研究を行いScienceに掲載された。
• 北太平洋（北緯47度）の底層でJAMSTECが発見した水温上昇は、時間を遡って解析するとおよそ４０年前の南極アデリーコースト沖周辺での海面フラックス変動の影響であることがわかった。

47Nに沿った船舶観測による1999年と1985年の水温差(Fukasawa et al.,2004)

海洋５０年再解析データセットを基にした４D-VAR全球全層海洋同化システムによる感度

より詳しく勉強したい人へ• データ同化についての教科書があります。
• 「データ同化」観測・実験とモデルを融合するイノベーション (京都大学出版会)• データ同化入門 ―次世代のシミュレーション技術―(朝倉書店)