JP2017533514A - 順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法 - Google Patents

Abstract

本発明は、順序再帰的フィルタリング３次元変分に基づく実測海洋環境データの同化方法を提供する。海洋環境実測データに対して前処理を行うステップと、目的関数値を計算するステップと、目的関数の勾配値を計算するステップと、目的関数の最小値を計算するステップと、実測データにおける空間マチルスケール情報を抽出するステップと、背景場データを更新し最終的なデータ同化解析場を形成するステップとを含む。本発明は、従来の再帰的フィルタリング３次元変分方法を改良し、異なるスケールの情報を順に同化することで、従来の３次元変分方法がマルチスケール情報を効率的に抽出できないという問題を効果的に克服し、且つ、高次再帰的ガウスフィルターを用い、高次再帰的フィルターのカスケード形式を並列接続構造に変換することで、再帰的ガウスフィルタリングの再帰的フィルタリングプロセスを並列実行でき、カスケードフィルターによる多くの問題を克服する。【選択図】図１

Description

本発明は、実測海洋環境データ同化方法に関し、具体的には、順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法に関する。

海洋は潜水機の活動空間であり、海洋内部の躍層や内部波などの小・中規模の海洋過程は、直接潜水機の航行安全を脅かし、そのタスクの順調な実行に影響を与える。従って、潜水機の周囲の航行環境場のリアルタイムな構築及び小・中規模の海洋過程の迅速な診断の実現は重要な意義がある。従来、海洋環境データの取得は主に２つの方法があり、１つの方法として、数値モデルを用いて計算して取得する方法があるが、該データは海洋環境場の規則的な特徴を近似的に反映しかできず、小・中規模の海洋過程を抽出しにくい。もう１つの方法として、センサを用いて直接測定する方法があるが、当該データは、海洋環境場に対して最も直接的で、精度の高い観測であり、重要な海洋情報を含んでいる。しかし、観測装置の制約と観測点の物理量のランダムな変動のため、観測結果には不可避な系統誤差とランダム誤差が生じる。これから分かるように、２種類の方法で取得されるデータはそれぞれの利点や欠点がある。

データ同化方法は、数値計算と数値予報サービスの発展に伴って発展してきた、観測データと理論モデルを組み合わせる方法である。データ同化方法は、一定の最適化基準と方法に基づき、異なる空間、異なる時間で、異なる観測手段を用いて取得した観測データを数学モデルと有機的に組み合わせ、モデルとデータが互いに協調する最適化関係を構築する。データ同化技術により、観測データに含まれる有効な情報を最大限に抽出し、分析・予報システムの特性を高めて改善することができる。従って、データ同化によってリアルタイムな観測データと理論モデル結果とを組み合わせ、両方の利点を取り入れることで、実際に近い海洋環境状態場を構築することができる。３次元変分は、常用の海洋資料同化方法である。３次元変分同化において、背景場誤差共分散行列は重要な役割を有し、背景場に対する解析場の修正程度を決める。従来の３次元変分同化方法において研究している主な問題の１つは、背景場誤差共分散行列をどのように良好に構築するのかということである。従来、３次元変分同化方法には主に２種類があり、１つの方法は相関長法であり、もう１つの方法は再帰的フィルタリング法である。相関長法は、相関スケールで背景場誤差共分散行列を構築するものである。しかし、当該方法は、実際の応用時に大きなランダム性があり、相関スケールの正確な推定値を提供することが難しい。それと共に、当該方法は巨大なメモリと計算オーバーヘッドが必要である。再帰的フィルタリング法は、誤差共分散行列を明示的に構築する必要がなく、相関長法に比べて、計算速度が速く、メモリを節約するという利点を有する。しかし、相関長法と再帰的フィルタリング法は同じ欠点があり、即ちマルチスケールの観測情報を効果的に抽出しにくい。従って、潜水機の水中航行の周囲海洋環境場に対する要求に対して、高効率の実測海洋環境データ同化方法を設計して、水中環境状態場のリアルタイムな構築とマルチスケール海洋過程の迅速な抽出を実現することは重要な実用的価値を有する。

本発明は、水中環境状態場のリアルタイムな構築とマルチスケール海洋過程の迅速な抽出を実現可能な順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法を提供することを目的とする。

本発明の目的は以下のような技術案により達成され、
当該技術案は、
ステップ１：先ず実測データにおける無効データと異常データを取り除き、次に同化に必要なフォーマットに変換する、海洋環境実測データの品質管理及び前処理ステップと、
ステップ２：前のステップで前処理した実測海洋環境データと背景場データを用い、並列再帰的フィルターで、且つフィルタリングパラメータを与えることによって、目的関数値を計算する、目的関数値の計算ステップと、
ステップ３：目的関数値と、正接線形・随伴モデルコンパイラにより生成される目的関数勾配のソースコードに基づき、目的関数の勾配値を計算する、目的関数の勾配値の計算ステップと、
ステップ４：目的関数値と目的関数の勾配値をＬＢＦＧＳ最小化アルゴリズムに入力し、目的関数の最小値を計算し、最小値に対応する状態変数が観測データの長波推定値である、目的関数の最小値の計算ステップと、
ステップ５：実測海洋環境データからステップ４で得られた長波推定値を減算して新たな実測データとし、フィルタリングパラメータを小さく調整し、同時にステップ２における実測海洋環境データの代わりに、新たな実測データを使用し、ステップ２〜ステップ４を繰り返し、実測データにおける空間マルチスケール情報を順に抽出する、実測データにおける空間マルチスケール情報の抽出ステップと、
ステップ６：背景場データの更新を行い、最終的なデータ同化解析場を形成する、背景場データの更新ステップと、を含む。

本発明は、順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化技術を提供する。従来の３次元変分方法による実測海洋環境データ同化技術に比べて、本発明は、フィルタリングパラメータσを逐次に変更することによって、異なるスケールの観測情報を順に抽出することを顕著な特徴とする。即ち、先ず、σを大きな値とするにより、フィルターの周波数窓の幅は狭くなり、長波（低周波）情報のみ通過でき、再帰的フィルタリング法で同化して長波の推定値を得ることができ、次に、実測海洋環境データから得た長波推定値を減算して新たな実測データとして、σを小さくするように適宜調整することにより、フィルターの周波数窓の幅は広くなり、中波（Ｍｅｄｉｕｍ Ｆｒｅｑｕｅｎｃｙ）情報が通過できるようになり、同化して中波の推定値を得て、このようにして、異なる波帯の情報を順次に抽出する。本発明に係る方法は、従来の３次元変分方法がマルチスケール情報を効率的に抽出できないという問題を克服することができる。

本発明は、
（１）従来の再帰的フィルタリング３次元変分方法を改良し、異なるスケールの情報を順に同化することで、従来の３次元変分方法がマルチスケール情報を効率的に抽出できないという問題を効果的に克服し、
（２）高次再帰的ガウスフィルターを用い、高次再帰的フィルターのカスケード形式を並列接続構造に変換することで、再帰的ガウスフィルタリングの再帰的フィルタリングプロセスを並列実行でき、カスケードフィルターによる多くの問題を克服する、といった利点を有する。

海洋環境実測データ同化のフローチャートである。 順序再帰的フィルタリング３次元変分実行のフローチャートである。

以下、例を挙げて本発明をさらに詳細に説明する。

図１に示すように、本発明は主に以下のようなステップを含む。

ステップ１：海洋環境実測データの前処理
実測海洋環境データを背景場データベースに同化するために、先ず実測データの品質管理を行い、実測データにおける無効データと異常データを取り除き、次に品質管理後の実測データを同化に必要なフォーマット、例えば常用のＡＳＣＩＩフォーマット又はバイナリのストレージフォーマットに変換する。

ステップ２：目的関数値の計算
前のステップで前処理した実測海洋環境データと背景場データを用い、並列再帰的フィルター、且つフィルタリングパラメータを与えることによって、目的関数値を計算する。

ステップ３：目的関数の勾配値の計算
目的関数値と正接線形及び随伴モデルコンパイラ（Ｔａｎｇｅｎｔ ａｎｄ ａｄｊｏｉｎｔ ｍｏｄｅｌ ｃｏｍｐｉｌｅｒ）により生成される目的関数勾配のソースコードに基づき、目的関数の勾配値を計算する。

ステップ４：目的関数の最小値の計算
目的関数値と目的関数の勾配値をＬＢＦＧＳ（ｌｉｍｉｔｅｄ ｍｅｍｏｒｙ Ｂｒｏｙｄｅｎ−Ｆｌｅｔｃｈｅｒ−Ｇｌｏｄｆａｒｂ−Ｓｈａｎｎｏ）最小化アルゴリズムに入力し、目的関数の最小値を計算し、最小値に対応する状態変数が観測データの長波推定値である。

ステップ５：実測データにおける空間マルチスケール情報の抽出
実測海洋環境データからステップ４で得られた長波推定値を減算して新たな実測データとし、フィルタリングパラメータを小さく調整し、同時にステップ２における実測海洋環境データの代わりに、新たな実測データを使用し、ステップ２〜ステップ４を繰り返し、実測データにおける中波と短波情報を順に抽出する。

ステップ６：背景場データの更新
背景場データの更新を行い、最終的なデータ同化解析場を形成する。

以下、図２を参照しながら本発明をさらに詳細に説明する。

ステップ１、実測海洋環境データの前処理及び品質管理
観測データに含まれるエラー及び非標準データをできるだけ除去するために、先ず実測海洋環境データに対して前処理及び品質管理を行う。本発明は、地域検査、重複深さ検査、深さ逆転検査、温度データ範囲検査、温度・塩度勾配検査、及び密度安定性検査等の６つのステップでデータの前処理及び品質管理を実行する。

ステップ１．１：地域検査
実測海洋環境データを潜水機の周辺の一定の範囲内にするために、地域検査を行う必要があり、本発明に用いられる検査方法は地域検索方法である。つまり、実測海洋環境データがユーザにより選択される緯度・経度の範囲内にあると、それを保存し、範囲外にあると、それを削除し、データ範囲はユーザのニーズに応じて設定することができる。

ステップ１．２：重複深さ検査
観測断面における重複深さデータを除去するために、重複深さ検査を行う必要がある。本発明に用いられる検査方法は線形検索方法である。つまり、各観測断面ごとに検査を行う。該断面に重複深さデータがあると検索したと、それをマーキングして、マンマシンインタラクション検査を行い、同一深さ層における海洋環境データと上下層との連続性及び整合性の原則に基づき、そのうちの１つの重複断面を除去する。

ステップ１．３：深さ逆転検査
観測データの中に偽の深さ反転（Ｉｎｖｅｒｓｅ）が存在するかどうかを検査するために、深さ逆転検査を行う必要がある。本発明に用いられる検査方法は線形検索方法である。深さ反転データを検索したと、それをマーキングし、更にマンマシンインタラクション検査を行い、上下層の観測資料と周囲の観測資料を組み合わせ、有効な深さ層の海洋環境データを保存する。

ステップ１．４：温度・塩度データ範囲検査
温度・塩度観測データに極限値が現れないことを確保するために、データ範囲検査を行う必要がある。本発明に用いられる検査方法は対比検査方法である。つまり、各海域の温度・塩度の最大と最小値の範囲を予め統計して与え、上記範囲を超える温度・塩度観測データがあると、それを削除する。

ステップ１．５：温度・塩度勾配検査
温度・塩度観測断面の勾配に極限値が現れないことを確保するために、勾配検査を行う必要がある。本発明に用いられる方法は対比検査方法である。つまり、各海域の温度・塩度の最大勾配値範囲を予め統計して与え、上記範囲を超える温度・塩度観測データがあると、それを削除する。

ステップ１．６：密度安定性検査
観測断面データの中に偽の密度反転（Ｉｎｖｅｒｓｅ）が現れないことを確保するために、密度安定性検査を行う必要がある。本発明に用いられる検査方法は線形検索方法である。つまり、上から下へ各層の密度値を比較し、下層の密度が上層の密度より小さいと、それをマーキングし、更にマンマシンインタラクション検査を行い、観測領域における海洋環境特性及び該海域の履歴観測資料を組み合わせ、密度反転の真実性を判断し、真実であると、それを保存し、真実ではないと、それを削除する。

ステップ１．７：上記品質管理後の海洋環境データの前処理
ステップ１．１〜ステップ１．６で品質管理された実測データを同化に必要なフォーマット、例えば常用のＡＳＣＩＩフォーマット又はバイナリのストレージフォーマットに変換し、次のデータ同化プロセスの入力インタフェースとする。

ステップ２：順序再帰的フィルタリング３次元変分の目的関数値の計算
目的関数の最小値を計算して、観測データ状態変数の最適な推定値を得るために、予め目的関数値を計算する必要がある。具体的な計算方法は下記の通りである。

ステップ２．１：目的関数の定義
計算量が非常に大きな背景誤差共分散行列の反転（Ｉｎｖｅｒｓｅ）が現れることを避けるために、本発明に用いる目的関数は下記の通りである。

ここで、Ｊは目的関数であり、ｗ＝Ｂ^−１（ｘ−ｘ_ｂ）であり、Ｂは背景場誤差共分散行列であり、ｘ_ｂは背景場（通常は予報場である）ベクトルであり、ｘは状態変数である。Ｏは観測場の誤差共分散行列であり、Ｈはモードグリッドから観測点への補間演算子を表し、ｄ＝ｙ−Ｈｘ_ｂであり、ｙは観測場ベクトルである。背景場は海洋数値モードによってデータ同化ステップで前進積分を行って得られる。上記目的関数の制御変数はｗである。

ステップ２．２：ｗの初期推測値を与える
上記目的関数の表現形態に基づき、最適化アルゴリズムにより最適なｗを得る前に、ｗの初期推測値を与える必要がある。本発明において、ｗの初期推測値を０に設定し、ステップ２．１におけるｗの計算式により、ｗの初期推測値を０に設定する場合、状態ベクトルｘの初期推測値をその背景場の値ｘ_ｂに間接的に仮定する。

ステップ２．３：並列再帰的フィルタリングによるＢｗの計算
本発明は、並列再帰的フィルターを用いてＢｗの計算過程をシミュレートし、並列再帰的フィルターの伝達関数Ｓ（ｚ）は、Ｓ（ｚ）＝Ｓ_１（ｚ）＋Ｓ_２（ｚ）＋Ｓ_３（ｚ）＋Ｓ_４（ｚ）に定義され、
当該フィルターが、実質的には四次ｖａｎ Ｖｌｉｅｔ再帰的ガウスフィルターが低次並列接続の形式に変換されるものである。ここで、関数Ｓ_ｉ（ｚ）（ｉ＝１、２、３、４）を下記の式により計算する。

ここで、ｅ_ｉ、ｆ_ｉ（ｉ＝１、２、３、４）は未定係数であり、計算によって得られるように、
ｅ_１＝２Ｒｅ（ｃ_１） 、ｆ_１＝−２Ｒｅ（ｃ_１ｄ_２）、
ｅ_２＝２Ｒｅ（ｃ_２） 、ｆ_２＝−２Ｒｅ（ｃ_２ｄ_４）、
ｅ_３＝−２Ｒｅ（ｃ_３ｄ_２） 、ｆ_３＝２Ｒｅ（ｃ_３）、
ｅ_４＝−２Ｒｅ（ｃ_４ｄ_４） 、ｆ_４＝２Ｒｅ（ｃ_４）、
ここで、Ｒｅ（・）は複素数の実部を取ることを示し、ｃ_ｉ、ｄ_ｉ（ｉ＝１、２、３、４）は未定係数であり、計算によって得られるように、

ここで、σ_０＝２であり、σ＞０であり、上付き文字「＊」は共役を表し、σはフィルタリングパラメータである。

ここで、Ｓ_−１（ｚ）（ｚ＝ｄ_１、ｄ_３）は四次ｖａｎ Ｖｌｉｅｔ再帰的ガウスフィルター伝達関数であり、その表現式は下記の通りである。

上式中のｂ_ｉ（ｉ＝１、２、３、４）及びαの表現式は下記の通りである。

ステップ２．４：観測増分の計算 ｄ＝ｙ−Ｈｘ_ｂ
Ｈは双線形補間演算子であり、Ｈによってモード空間を観測空間にマッピングする。本発明では、Ｈは双線形距離反比例法で計算され、即ちモードの１つのセルにおける４つの格子点の状態値を該セル内の観測位置に補間する時、まず４つの格子点と観測点との物理的距離を計算し、次に距離の反比例を重みとして４つのモード格子点における状態値に対して加重平均を行うことにより、観測位置の状態値を得る。

ステップ２．５：目的関数値Ｊ（ｗ）の計算
ステップ２．２、２．３、２．４で算出されたｗ、Ｂｗ、ｄの値をステップ２．１で決定された目的関数に代入し、目的関数値を算出する。

ステップ３：目的関数の勾配値の計算
目的関数の最小値を計算し、更に制御変数の最適な推定値を得るために、ｗに対する目的関数の勾配値を予め計算する必要がある。具体的な計算方法は下記の通りである。

ステップ３．１：目的関数の勾配の定義
本発明に採用した目的関数勾配の計算式は、下記の通りである。

ここで、上付き文字「Ｔ」は転置を表し、ほかの符号の意味はステップ２．１で定義された目的関数計算式と同様である。

ステップ３．２：正接線形・随伴モデルコンパイラ（Ｔａｎｇｅｎｔ ａｎｄ ａｄｊｏｉｎｔ ｍｏｄｅｌ ｃｏｍｐｉｌｅｒ、ＴＡＭＣ）による目的関数勾配のソースコードの生成
目的関数値を計算する際に、再帰的フィルタリングアルゴリズムを用いる必要があるため、定義を直接用いて目的関数の勾配を計算することは困難である。計算の便宜上、目的関数値を計算するＦｏｒｔｒａｎプログラム計算ソースコードをＴＡＭＣコンパイラに入力し、当該コンパイラは目的関数勾配のＦｏｒｔｒａｎソースコードを自動的に返す。

ステップ３．３：目的関数の勾配値∇Ｊの計算
ステップ３．２で得られた、目的関数勾配を計算するＦｏｒｔｒａｎソースコードを実行し、ｗに対する目的関数の勾配値を得る。

ステップ４：目的関数の最小値の計算
制御変数の最適な推定値を得るために、目的関数の最小値を計算する必要がある。具体的な計算方法は下記の通りである。

ステップ４．１：目的関数値及び目的関数の勾配値をＬＢＦＧＳ（ｌｉｍｉｔｅｄ ｍｅｍｏｒｙ Ｂｒｏｙｄｅｎ−Ｆｌｅｔｃｈｅｒ−Ｇｌｏｄｆａｒｂ−Ｓｈａｎｎｏ）最小化アルゴリズムに入力する。

ステップ４．２：勾配値と目的関数値に基づいて最小化アルゴリズムにより制御変数での実現可能な降下方向と検索ステップ幅を自動的に決定する。実現可能な降下方向と検索ステップ幅は、最適化アルゴリズムにおける次回の線形検索の制御変数の初期値を決定するためのものである。

ステップ４．３：上記の最適化アルゴリズムで得られた解析場を用い、制御変数ｗを更新し、目的関数の勾配値∇Ｊを改めて計算し、この時の勾配が収束基準を満たすかどうかを判断する。満たさないと、収束基準を満たすまで、ステップ４．１〜４．２を繰り返して計算し続ける。満たすと、計算を終了し、この時に得られた最小値に対応する状態変数ｘの値は背景場の最適な推定値、即ち解析場である。当該最適な推定値は最小二乗法の意味の下で最適である。

本発明に用いられる収束基準は、

であり、
ここで、∇Ｊ_０は１回目のイテレーションにおける目的関数勾配値である。

ステップ５：実測データにおける空間マルチスケール情報の抽出
観測におけるマルチスケール海洋環境情報を抽出するために、順序再帰的フィルタリングプロセスを実行する必要がある。具体的な計算方法は下記の通りである。

ステップ５．１：ステップ２．３における再帰的ガウスフィルタリングにより、ｘ^ｎｅｗ＝Ｂｗ^ｎｅｗを計算する。

ステップ５．２：観測資料からステップ５．１で得られた長波推定値を減算して新たな観測とし、即ち、ｄ＝ｄ−Ｈｘ^ｎｅｗである。ここでの観測は、ステップ１で得られた前処理後の観測値からステップ５．１で得られた長波推定値を減算した残量である。

ステップ５．３：σの値を小さく調整し、即ちσ＝σ^（ｍ）＜σ^{（ｍ−１）}とする。

σは指数減衰の形式で小さくなり、例えばσ＝σ／２^Ｍ−１であり、ここでＭはサイクル回数であり、一般的には８回である。

ステップ５．４：終了条件を満たすかどうかを判断する。満たさないと、現在の状態変数ｘを記録し、ステップ２に移行する。終了条件を満たすと、ステップ６に移行する。

サイクル回数が設定された回数Ｍを満たすかどうかを終了条件とする。サイクル回数Ｍとして、一般的には８回であり、この時のα値は初期σ値の０．７％まで小さくなる。

ステップ６：背景場データの更新
ステップ５でＭ回のサイクルのうちの各サイクルにおいて得られた状態変数ｘを加算し、最終的な解析場を得て、一回のデータ同化プロセスを完了する。該解析場を現在の海洋環境状態を反映するデータ場とすることができる。

Claims (4)

1. ステップ１：先ず実測データにおける無効データと異常データを取り除き、次に同化に必要なフォーマットに変換する、海洋環境実測データの品質管理及び前処理ステップと、
   ステップ２：前のステップで前処理した実測データと背景場データを用い、並列再帰的フィルターで、且つフィルタリングパラメータを与えることによって、目的関数値を計算する、目的関数値の計算ステップと、
   ステップ３：目的関数値と正接線形・随伴モデルコンパイラにより生成される目的関数勾配のソースコードに基づき、目的関数の勾配値を計算する、目的関数の勾配値の計算ステップと、
   ステップ４：目的関数値と目的関数の勾配値を準ニュートン最小化アルゴリズムに入力し、目的関数の最小値を計算し、最小値に対応する状態変数が観測データの長波推定値である、目的関数の最小値の計算ステップと、
   ステップ５：実測海洋環境データからステップ４で得られた長波推定値を減算して新たな実測データとし、フィルタリングパラメータを小さく調整し、同時にステップ２における実測海洋環境データの代わりに、新たな実測データを使用し、ステップ２〜ステップ４を繰り返し、実測データにおける中波及び短波情報を順に抽出する、実測データにおける空間マルチスケール情報の抽出ステップと、
   ステップ６：背景場データの更新を行い、最終的なデータ同化解析場を形成する、背景場データの更新ステップと、
   を含むことを特徴とする順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法。
2. 前記目的関数値の計算は、具体的には、下記のステップ２．１〜ステップ２．５を含み、
   ステップ２．１では、目的関数を定義し、
   

   

   
   ここで、Ｊは目的関数であり、ｗ＝Ｂ^−１（ｘ−ｘ_ｂ）であり、Ｂは背景場誤差共分散行列であり、ｘ_ｂは背景場ベクトルであり、ｘは状態変数であり、Ｏは観測場誤差共分散行列であり、Ｈはモードグリッドから観測点への補間演算子を表し、ｄ＝ｙ−Ｈｘ_ｂであり、ｙは観測場ベクトルであり、背景場は海洋数値モードによってデータ同化ステップで前進積分を行って得られ、
   ステップ２．２では、ｗの初期推測値を与え、
   ｗの初期推測値を０に設定し、状態ベクトルｘの初期推測値をその背景場の値ｘ_ｂに間接的に仮定し、
   ステップ２．３では、並列再帰的フィルタリングによってＢｗを計算し、
   並列再帰的フィルターの伝達関数Ｓ（ｚ）は、Ｓ（ｚ）＝Ｓ_１（ｚ）＋Ｓ_２（ｚ）＋Ｓ_３（ｚ）＋Ｓ_４（ｚ）に定義され、
   当該フィルターは、実質的には四次ｖａｎ Ｖｌｉｅｔ再帰的ガウスフィルターから低次並列接続の形式に変換されるものであり、ここで、関数Ｓｉ（ｚ）、ｉ＝１、２、３、４を、下記の式により計算し、
   

   

   
   ここで、ｅ_ｉ、ｆ_ｉ、ｉ＝１、２、３、４は未定係数であり、計算によって下記の値を得て、
   ｅ_１＝２Ｒｅ（ｃ_１）、ｆ_１＝−２Ｒｅ（ｃ_１ｄ_２）、
   ｅ_２＝２Ｒｅ（ｃ_２）、ｆ_２＝−２Ｒｅ（ｃ_２ｄ_４）、
   ｅ_３＝−２Ｒｅ（ｃ_３ｄ_２）、ｆ_３＝２Ｒｅ（ｃ_３）、
   ｅ_４＝−２Ｒｅ（ｃ_４ｄ_４）、ｆ_４＝２Ｒｅ（ｃ_４）、
   ここで、Ｒｅ（・）は複素数の実部を取ることを示し、ｃ_ｉ、ｄ_ｉ（ｉ＝１、２、３、４）は未定係数であり、計算によって下記の値を得て、
   

   

   
   ここで、σ０＝２であり、σ＞０であり、上付き文字「＊」は共役を表し、σはフィルタリングパラメータであり、ここで、Ｓ_−１（ｚ）ｚ＝ｄ_１、ｄ_３は四次ｖａｎ Ｖｌｉｅｔ再帰的ガウスフィルター伝達関数であり、その表現式は、下記の通りであり
   

   

   
   上式中のｂ_ｉ、ｉ＝１、２、３、４及びαの表現式は、下記の通りであり
   

   

   
   ステップ２．４では、観測増分をｄ＝ｙ−Ｈｘ_ｂにより計算し、
   ステップ２．５では、目的関数値Ｊ（ｗ）を計算し、
   ステップ２．２、２．３、２．４で算出されたｗ、Ｂｗ、ｄの値をステップ２．１で決定された目的関数に代入し、目的関数値を算出することを特徴とする請求項１に記載の順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法。
3. 前記目的関数の勾配値の計算は、具体的には、下記のステップ３．１〜ステップ３．３を含み、
   ステップ３．１では、目的関数の勾配を定義し、
   

   

   
   ここで、上付き文字「Ｔ」は転置を表し、
   ステップ３．２では、正接線形・随伴モデルコンパイラにより目的関数勾配のソースコードを生成し、
   ステップ３．３では、目的関数の勾配値∇Ｊを計算し、
   ステップ３．２で得られた、目的関数勾配を計算するＦｏｒｔｒａｎソースコードを実行し、ｗに対する目的関数の勾配値を得ることを特徴とする請求項２に記載の順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法。
4. 前記実測データにおける空間マルチスケール情報の抽出は、具体的には、下記のステップ５．１〜ステップ５．４を含み、
   ステップ５．１では、ステップ２．３における再帰的ガウスフィルタリングによりｘ^ｎｅｗ＝Ｂｗ^ｎｅｗを計算し、
   ステップ５．２では、観測資料からステップ５．１で得られた長波推定値を減算して新たな観測とし、即ち、ｄ＝ｄ−Ｈｘ^ｎｅｗであり、
   ステップ５．３では、σの値を小さく調整し、即ちσ＝σ^（ｍ）＜σ^{（ｍ−１）}とし、
   σは指数減衰の形式で小さくなり、σ＝σ／２^Ｍ−１であり、ここでＭはサイクル回数で、８回であり、
   ステップ５．４では、終了条件を満たすかどうかを判断し、満たさないと、現在の状態変数ｘを記録し、目的関数値の計算に移行し、終了条件を満たすと、背景場データの更新に移行することを特徴とする請求項３に記載の順序再帰的フィルタリング３次元変分に基づく実測海洋環境データ同化方法。

Images