Freak wave estimation by WAVEWATCH III and HOSM

Size: px

Start display at page:

Download "Freak wave estimation by WAVEWATCH III and HOSM"

ありさえいさか
4 years ago
Views:

1 2017 年 9 月 25 日第 2 回理研データ同化ワークショップアンサンブル 4 次元変分法を活用した巨大波浪の再現藤本航早稲田卓爾東京大学新領域創成科学研究科海洋技術環境学専攻 1

導入 : フリーク波とは海洋に突発的に生じる巨大波波浪場の平均的な振幅 (H s : 有義波高 ) を大きく超える波 H > 2 or η crest > 1.3 H s H s (Dysthe et al. http://folk.uio.

2 導入 : フリーク波とは海洋に突発的に生じる巨大波波浪場の平均的な振幅 (H s : 有義波高 ) を大きく超える波 H > 2 or η crest > 1.3 H s H s (Dysthe et al. 工学的 : 船舶海洋構造物の安全 Taylor, P The Draupner wave of 1 st January 1995 ICMS 北海石油プラットホームでの観測日本近海でも観測 2

3 観測データからの波形再構成フリーク波の推定物理の研究船舶海洋構造物発電機の事故解析逆問題データ同化再解析津波インバージョンの深海波 ver. フリーク波インバージョン想定する観測データブイ波高計レーダー LIDAR など津波インバージョン藤井佐竹 (2011) 建築研究所 HP Nouguier and Guerin (2014) 船舶からの LIDAR 波浪計測 3

4 本研究の概要 1. フリーク波を再構成するには非線形性を考慮する必要がある 2. 波浪の非線形時間発達を解ける手法があるが通常の 4 次元変分法 ( アジョイント法 ) を適用するには課題が多い 3. アンサンブル 4 次元変分法を用いてフリーク波を再構成する 4

波浪の非線形性とフリーク波の関係連続の式水底境界条件運動学的自由境界条件力学的自由境界条件基礎方程式系 (Euler 方程式 ) 2 φ z 2 = 2 φ Φ z = 0 d z η at z = η t = Φ η + 1 + η 2 W at z = η Φ t = 1 2 Φ 2 + gη + 1 2 1 + η 2 W 2 at z = η O z 自由表面 z =

5 波浪の非線形性とフリーク波の関係連続の式水底境界条件運動学的自由境界条件力学的自由境界条件基礎方程式系 (Euler 方程式 ) 2 φ z 2 = 2 φ Φ z = 0 d z η at z = η t = Φ η η 2 W at z = η Φ t = 1 2 Φ 2 + gη η 2 W 2 at z = η O z 自由表面 z = η(x, y) z = x, y 表面上の速度ポテンシャル Φ = φ(x, y, z = η, t) 表面上の鉛直速度 W = z φ(x, y, z = η, t) 水平方向のナブラ = ( x, y ) 非線形性によって振幅が局所的に増大 = 変調不安定フリーク波のモデル (Benjamin & Feir 1967) 水槽実験 Chabchoub et al. (2012) 5

6 Higher Order Spectral Method (HOSM) Dommermuth & Yue (1987), West et al (1987) Euler 方程式を解く海洋波浪の DNS ラプラス方程式の解から空間微分を FFT で評価 ( スペクトル法 ) 高速高精度任意の次数の非線形性広いスペクトル幅を考慮自由境界条件について Φ, W は z=0 周りにテイラー展開 ( 摂動展開 ) 表面上の速度ポテンシャル Φ = φ(x, y, z = η, t) 表面上の鉛直速度 W = z φ(x, y, z = η, t) 水平方向のナブラ = ( x, y ) η t = w 1 Φ η + w 2 +w 3 +w 1 η 2 Φ t = gη 1 2 Φ w 1 w 1 + w 1 w 2 φ 0 1 = Φ φ 2 0 = η φ 0 1 z φ 3 0 = η φ 0 2 z η2 1 φ 0 2! z w (1) = φ 0 1 z w (2) = φ 2 0 z + η 2 1 φ 0 z 2 w (3) = φ 3 0 z + η 2 2 φ 0 z 2 + η2 3 1 φ 0 2! z 3 HOSM で計算されたフリーク波 6

7 観測値と波浪予報モデルを用いた波浪場の再構築 WAVEWATCH III などの波浪予報モデル ( パワースペクトルのみ出力 ) 初期波浪場 ηƹ 0 k = A k exp(2πi ψ(k)) パワースペクトルの背景値 S k 観測データ y HOSM データ同化振幅 A k 位相 ψ(k) 波浪予報モデルはパワースペクトルしか出力しない波浪の DNS と組み合わせることで観測データから位相も復元する ( 振幅も調整 ) 4 次元変分法を適用 ( 空間的に小量な時系列データを想定 )

8 4 次元変分法の適用とコスト関数の定式化 WAVEWATCH III などの波浪予報モデルによるパワースペクトル S k y: 観測時系列 H η 0 : モデル内時系列パワースペクトルの背景値 S k 観測データ y 制御変数 η 0 : 初期波形のフーリエ係数 L η 0 = λ 2 η 0 D 1 η H η 0 y 2 HOSM HOSMで力学を考慮することで観測演算子を時間方向に拡張 Dは対角行列とすることができ逆行列 D 1 を容易に計算メモリの節約初期波形がフーリエ係数で表されるため ( 異なる成分波は無相関 ) diag D = S k S k を背景誤差共分散とみなしている λはチューニングパラメタ L2 正則化 (Tikhonov 正則化 ) 8

9 Higher Order Spectral Method (HOSM) Dommermuth & Yue (1987), West et al (1987) Euler 方程式 (NLS と同じ基礎方程式系 ) を解く海洋波浪の DNS ラプラス方程式の解から空間微分を FFT で評価 ( スペクトル法 ) 高速高精度任意の次数の非線形性広いスペクトル幅を考慮自由境界条件について Φ, W は z=0 周りにテイラー展開 ( 摂動展開 ) η t = w 1 Φ η + w 2 +w 3 +w 1 η 2 Φ t = gη 1 2 Φ w 1 w 1 + w 1 w 2 φ 0 1 = Φ φ 2 0 = η φ 0 1 z φ 3 0 = η φ 0 2 z η2 1 φ 0 2! z w (1) = φ 0 1 z w (2) = φ 2 0 z + η 2 1 φ 0 z 2 w (3) = φ 3 0 z + η 2 2 φ 0 z 2 アジョイント法を適用するのは困難摂動展開を多用する (Aragh & Nwogu 2008) フォワード計算に時間がかかるのでアジョイント方程式にも時間がかかるアンサンブル 4 次元変分法 + η2 3 1 φ 0 2! z 3 HOSM で計算されたフリーク波 9

10 アンサンブル 4 次元変分法 (1) η [m] H η 0 は初期波形から観測時系列を出力する ( 力学を考慮し時間方向に拡張された ) 観測演算子初期波形 η 0 の摂動ベクトル P = p m が張る空間 ( モデル空間 ) δy = HP H の線形化 H ( ヤコビアン ) 観測値 y の摂動ベクトル δy = δy m が張る空間 ( 観測空間 ) δy m = H η 0 + p m H η 0 Hp m 摂動 p m cos 波 sin 波モデル予測値の摂動 δy m = H η 0 + p m H η 0 摂動を与えられたアンサンブルラン η 0 m = η 0 + p m 時間発展モデル予測値 H η 0 m 10

11 アンサンブル 4 次元変分法 (2) 初期波形 η 0 の摂動ベクトル P = p m が張る空間 ( モデル空間 ) コスト関数の勾配 L = H T H η 0 y P が張る空間上で勾配の残差を最小化する P T L = P T H T 初期値の更新 Ps H η 0 + Ps y = 0 δy = HP ミスフィット y H η 0 観測値 y の摂動ベクトル δy = δy m が張る空間 ( 観測空間 ) H T を解析的に求めるのが通常の 4 次元変分法 ( アジョイント法 ) 次の初期値を摂動 p m の線形和で表す s は重み係数 η 0 + Ps η 0 P T H T H η 0 + HPs y + O( s 2 ) = 0 δy T H η 0 + δys y = 0 δy T δys = δy T y H η 0 各イタレーションで観測演算子 H η 0 のヤコビアンH をアンサンブル摂動計算によって η 0 の周りで線形近似アジョイント方程式が不要並列化が容易 11

12 アンサンブル 4 次元変分法の種類 Maximum Likelihood Ensemble Filter (MLEF, Zupanski 2005) 摂動はヘッシアン行列のコレスキー分解から与えるヘッシアン行列の固有値が大きいモードを調整する 4DVAR とアンサンブルカルマンフィルター (ETKF) の折衷摂動から勾配とヘッシアン行列を近似 En4DVAR (Liu et al. 2008, 2009) 摂動を固定する気象では発達する擾乱モードを常に追いかけているのでそれを摂動に与えるヘッシアン行列を直接求めない ( 最急降下法 ) 予報誤差共分散行列に局所化を用いる Adjoint-free 4DVAR (a4dvar, Yaremchuk et al. 2009, 2016) 摂動はモデルデータや観測値の EOF から与える固有値が大きいモードを与える発達する擾乱モードについての事前情報が無く観測がスパースな海洋の問題に適している摂動で得られたデータを再利用する全てのモードを調整する保存系に近くスペクトルが広い問題に適している GMRES と関係摂動から勾配とヘッシアン行列を近似 12

13 1 ミスフィットのパワースペクトルから摂動を選択本手法ではミスフィット y H η 0 のパワースペクトルを計算しそのピーク付近に相当するフーリエモードを摂動 P として加えるミスフィットパワースペクトル S(ω) 青 : 現在のイタレーション黒 : ひとつ前のイタレーションミスフィットパワースペクトルS(k) 線型分散関係で変換 ω 2 = gk 13

14 具体的な計算設定 HOSM 3 次非線形時間 : 50 Tp ( ピーク周期 ) JOSNWAP γ=3.3 波形勾配 : Hs/2*kp=0.11 比較的強い非線形性周期境界条件の影響 128 λ p の領域で真値を生成 32 λ p 領域内で推定真値の水位の時空間プロファイル 128 λ p の領域のうち 64λ p の範囲を図示 η crest = 1.51 H m0 ノイズ : 正規乱数標準偏差は真の値の 10% 32 λ p 領域内で同化青 : 線形推定値赤 : 真値特にフリーク波につながる波群の周辺で精度が不足 14

15 λ=0.005, Nens=50, e TOL = 0.2 の場合の解析値初期波形のプロファイル真値解析値 (200 イタレーション後 ) 一定の範囲内において解析値は真値とよく一致推定可能範囲は分散関係による群速度によって決まる時系列の長さがN T p であった場合推定可能範囲はおよそ N/2 λ p 正則化によって境界付近は減衰される特に高波数成分周期境界条件の影響を緩和 15

16 並列化とスタッキングの効果 ( λ=0.005 ) 赤 :e TOL = 0 ( 常にスタックしない ) 青 :e TOL = 0.2 黒 :e TOL = ( 常にすべてのアンサンブルをスタックする ) ー実線 :Nens=50 -- 破線 :Nens=10 -- 破線 :Nens=10 ー実線 :Nens=50 e TOL が大きいほどスタックされる摂動が多くなり収束が早いが不安定になるアンサンブル数を増やすことによってイタレーションを減らせる ( 並列化 ) 初期値の変化 s n が小さいときに摂動がスタックされ s n が大きいときに精度の悪い摂動が消去される ( スタックされた摂動の数がノコギリ状に変化する理由 ) 16

17 結論以下の要素により水位時系列からの波形再構成が可能であり推定を並列化できることを示した 1 HOSM アンサンブル計算による 4 次元変分法 2 非線形スピンアップによる力学的整合性の確保 3 事前予測パワースペクトルによる正則化安定性向上 4 ミスフィットのパワースペクトルから摂動を選択 5 過去のイタレーションで得られた摂動の近似精度の確認と再利用 ( スタッキング ) 今後の課題多方向のケースについてより広い領域設定の中でインバージョンを行うフリーク波再構成に適した観測システムを提案する 17

パソコンシミュレータの現状

パソコンシミュレータの現状第 2 章微分偏微分, 写像豊橋技術科学大学森謙一郎 2. 連続関数と微分工学において物理現象を支配する方程式は微分方程式で表されていることが多く, 有限要素法も微分方程式を解く数値解析法であり, 定式化においては微分積分が一般的に用いられており. 数学の基礎知識が必要になる. 図 2. に示すように, 微分は連続な関数 f() の傾きを求めることであり, 微小なに対して傾きを表し, を無限に