(57)【要約】【課題】本発明は経験的直交関数分解法による静止海色衛星データの再構築方法を提供する。【解決手段】この経験的直交関数分解法による静止海色衛星データの再構築方法の手順は次のとおりである。１）元の海色リモートセンシング領域に対して集中度検出及びエッジ検出により異常ピクセルを除去する；２）経験的直交関数分解法においてラプラス平滑化フィルタリングによって時間共分散行列に対して平滑化フィルタリングを行う；３）経験的直交関数分解法により取得された時間特徴モード値をサブ時間特徴モードに分解し、次に経験的モード分解を行い、最後に補間して、その元の精度を維持し時間と空間が完全にカバーされた水色リモートセンシング再分析データセットを取得する。本方法の利点は、静止海色衛星データの高時間分解能の特徴及びデータ欠損に連続した高欠損率と完全欠損があるという状況を十分に考慮し、静止海色衛星データの再構築に適され、精度が高く、解釈率が高いことにある。本発明は、静止海色衛星データの再構築にとって非常に重要な実際の応用価値がある。  【選択図】図１
2021533427.tif 000002
【特許請求の範囲】
【請求項１】
経験的直交関数分解法による静止海色衛星データの再構築方法であって、
再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像を前処理して、異常ピクセルをデータ欠損点とするとともに、データセットにおける時間的又は空間的な欠損率が上限より高い極高欠損率のデータを除去するステップＳ１と、
前処理後の海色リモートセンシングデータセットの時空間フィールドで二次元データ行列を構築し、経験的直交関数の分解と再構築を行って、各反復分解前にラプラス平滑化フィルタリングにより時間共分散行列をフィルタリングし、再構築後に、最適な保留モード数と、対応する時間特徴モード及び空間特徴モードとを取得するステップＳ２と、
前記時間特徴モードを列によって特徴モードベクトルに分解し、各列の特徴モードベクトルをデータの属する時点によって複数のサブ時間特徴モードに分解し、さらに各サブ時間特徴モードに対して経験的モード分解を行い、三次スプライン補間法により分解後の固有モード関数成分及び残量における極高欠損率のデータに対応する箇所のデータを補完し、改めて重ね合わせてマージして時系列的に連続した完全時間特徴モードを取得するステップＳ３と、
Ｓ２で取得した空間特徴モード及びＳ３で取得した完全時間特徴モードを利用して、時間と空間が完全にカバーされた静止海色衛星データセットを再構築して取得するステップＳ４とを含むことを特徴とする再構築方法。
【請求項２】
前記ステップＳ１において、それぞれの静止軌道海色衛星リモートセンシング画像に対して、
再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像に対して、それぞれ画像における各有効ピクセルを中心として辺長Ｌの空間窓を設定するステップＳ１１、
各有効ピクセルのそれぞれに対して、該ピクセルに対応する空間窓におけるすべての有効ピクセルをサンプルデータセットＯとして、該ピクセルの集中度Ｏｃｏｎｃを算出し、集中度Ｏｃｏｎｃは、下式の通りであり、
【数１】
2021533427.tif 000080
式中、Ｏｉが集中度を算出すべき現在ピクセル値であり、ＯｍがサンプルデータセットＯの中央値であり、標準偏差δの算出式がδ=ｋ*ＭＡＤ(Ｏ)であり、ｋが比例因子定数であり、ＭＡＤ（Ｏ）がサンプルデータセットＯの中央値絶対偏差であるステップＳ１２、
各有効ピクセルのそれぞれに対して、該ピクセルのエッジ検出値Ｏｐｒｏｘを算出し、該ピクセルの周囲に無効ピクセルがあれば、Ｏｐｒｏｘが正数値を取り、該ピクセルの周囲に無効ピクセルがなければ、Ｏｐｒｏｘがゼロを取るステップＳ１３、
各有効ピクセルのそれぞれに対して、該ピクセルの異常検出結果Ｏfｉｎａｌを算出し、異常検出結果Ｏfｉｎａｌは、下式の通りであり、
【数２】
2021533427.tif 000081
式中、wｐｒｏｘがエッジ検出重みであり、wｃｏｎｃが集中度検出重みであり、且つwｐｒｏｘ + wｃｏｎｃ =１を満足するステップＳ１４、
予め設定された異常ピクセル検出閾値Ｏｔｈｒｅｓｈｏｌｄに基づき、Ｏfｉｎａｌが閾値Ｏｔｈｒｅｓｈｏｌｄより大きいピクセルを異常ピクセルとして、異常ピクセルをデータ欠損点とするＳ１５に従って異常ピクセル識別処理を行うことを特徴とする請求項１に記載の再構築方法。
【請求項３】
前記ステップＳ１において、前記極高欠損率のデータの除去方法は、
時間的に、いずれか１つの時点の静止軌道海色衛星リモートセンシング画像の欠損率が欠損率の上限より大きい場合、該時点に対応する画像を除去し、空間的に、すべての静止軌道海色衛星リモートセンシング画像におけるいずれか１つの位置点の時系列データの欠損率が欠損率の上限より大きい場合、該位置点のすべての画像でのピクセル値を除去することを特徴とする請求項１に記載の再構築方法。
【請求項４】
前記欠損率の上限が９５％であることを特徴とする請求項３に記載の再構築方法。
【請求項５】
前記ステップＳ２の具体的な方法は、
前処理後の海色リモートセンシングデータセットの時空間フィールドをｍ×ｎの二次元データ行列Ｘ０とし、ｍが静止軌道海色衛星リモートセンシング画像における空間位置点の個数であり、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、且つｍ>ｎであり、行列Ｘ０のｍ行がそれぞれｍ個の異なる空間位置点の時系列値を代表し、ｎ列がそれぞれｎ個の異なる時点のすべての空間位置点の値を代表するステップＳ２１と、
行列Ｘ０におけるすべての有効データの平均値
2021533427.tif 000082
を算出し、
2021533427.tif 000083
にし、行列Ｘから一部の有効データを交差検証セットＸｃとしてランダムに抽出し、行列Ｘにおけるすべてのデータ欠損点及び交差検証セットＸｃに抽出された要素にいずれも０の値を割り当てるステップＳ２２と、
保留モード数Ｐを反復して、異なるＰ値における空間特徴モードＵＰ、時間特徴モードＶＰ及び交差検証セットの再構築値と元の値との二乗平均平方根誤差Ｒを算出し、反復するたびにいずれもステップＳ２３１〜Ｓ２３４を実行するステップＳ２３と、
保留モード数Ｐ=１, ２,…, ｋＭＡＸにおける二乗平均平方根誤差ＲＰを比較し、ＲＰに対応する最小Ｐ値を最適な保留モード数
2021533427.tif 000084
として、最適な時間特徴モード
2021533427.tif 000085
及び最適な空間特徴モード
2021533427.tif 000086
を取得するステップＳ２４とを含み、
ステップＳ２３１〜Ｓ２３４は以下のとおりであり、
ステップＳ２３１：行列Ｘに対して特異値分解を行い、特徴ベクトル分解によってＰ個の最大特異値及び特異ベクトルを算出し、更にｍ×Ｐの空間特徴モードＵＰ、Ｐ×Ｐの特異値行列SＰ及びｎ×Ｐの時間特徴モードＶＰを取得し、計算式は、下式の通りであり、
【数３】
2021533427.tif 000087
【数４】
2021533427.tif 000088
式中、uｉが空間特徴モードＵｐの第ｉ列であり、ｖｉが時間特徴モードＶＰの第ｉ列であり、ρｉがuｉ及びｖｉに対応する特異値であり、ｉ = １,...,Ｐであり、ＸＴＸがｎ×ｎの時間共分散行列であり、
ステップＳ２３２：ラプラス平滑化フィルタリングによって前記時間共分散行列の列をフィルタリングし、次に前記時間共分散行列の行をフィルタリングし、時間共分散行列の行又は列がいずれもｎ個の時点データを含む時系列データセットＣ={ｃ１,ｃ２,…，ｃｎ}であると見なされ、Ｃに対応する時間セットがＴ={ｔ１,ｔ２,…, ｔｎ}であり、ｃｌがｌ番目の時点ｔｌの時系列データであり、ｌ=１,２,…, ｎであり、且つ
2021533427.tif 000089
を満足し、ラプラス平滑化フィルタリングの式は、下式の通りであり、
【数５】
2021533427.tif 000090
【数６】
2021533427.tif 000091
【数７】
2021533427.tif 000092
式中、ｃｌ(ｓ)がラプラス平滑化フィルタリング前のデータｃｌを示し、ｃｌ(ｓ+１)がラプラス平滑化フィルタリング後のデータｃｌを示し、αが平滑化強度であり、０＜α＜１であり、
ステップＳ２３３：１回のラプラス平滑化フィルタリングを完了した後、行列Ｘにおけるすべてのデータ欠損点の再構築値を算出し、行列Ｘの第ｉ行第ｊ列の欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数８】
2021533427.tif 000093
式中、(uｐ)ｉが空間特徴モードＵＰの第ｐ列のｉ番目の要素であり、(ｖｐＴ)ｊが時間特徴モードＶＰの第ｐ列のｊ番目の要素であり、ρｐが対応する特異値であり、
ステップＳ２３４：前記交差検証セットにおけるすべての要素の再構築値と真理値との二乗平均平方根誤差Ｒを算出し、二乗平均平方根誤差Ｒは、下式の通りであり、
【数９】
2021533427.tif 000094
式中、Ｎが前記交差検証セットにおけるデータの総個数であり、Ｘｒｅｐ、Ｘｃｐがそれぞれ交差検証セットにおけるｐ番目の要素の再構築値、真理値であり、
二乗平均平方根誤差Ｒが誤差閾値より大きい場合、Ｓ２３３において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、次に更新後の行列ＸでステップＳ２３１〜Ｓ２３３を繰り返し、二乗平均平方根誤差Ｒが閾値以下に収束するまで、繰り返しを停止して、収束された二乗平均平方根誤差Ｒを現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記すことを特徴とする請求項１に記載の再構築方法。
【請求項６】
前記ステップＳ２３において、保留モード数Ｐの反復範囲が［１,ｋＭＡＸ］であり、最大反復回数がｋＭＡＸ≦ｍｉｎ（ｍ,ｎ）であることを特徴とする請求項５に記載の再構築方法。
【請求項７】
前記ステップＳ２３において、ステップＳ２３１〜Ｓ２３３を繰り返して特異値分解及びラプラス平滑化フィルタリングを行う回数が最大回数に達すると、繰り返しを停止して、最終的な二乗平均平方根誤差Ｒを現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記すことを特徴とする請求項５に記載の再構築方法。
【請求項８】
前記ステップＳ３の具体的な方法は、
最適な時間特徴モード
2021533427.tif 000095
における各列の特徴モードベクトルを分解し、分解時に日を周期として異なる日における同じ時点に属するデータを１つのサブ時間特徴モードに分解し、第ｉ列の特徴モードベクトルｖｉ（ｔ）は、次のように分解され、
【数１０】
2021533427.tif 000096
式中、Ｋが静止軌道海色衛星リモートセンシング画像の日内データの頻度総数であり、ｖｉｋ（ｄｋ）がｖｉ（ｔ）におけるｋ番目の時点に対応するサブ時間特徴モードであり、ｖｉ（ｔ）におけるすべての日付のｋ番目の時点のデータを含み、
2021533427.tif 000097
であり、ｋ=１,２,…, Kであり、ｄｋがｖｉｋ（ｄｋ）におけるデータに含まれる日付を示し、ｄｋに前記極高欠損率のデータに対応する欠損日付を含まないステップＳ３１と、
各列の特徴モードベクトルにおける各サブ時間特徴モードｖｉｋ（ｄｋ）に対して経験的分解を行い、
【数１１】
2021533427.tif 000098
式中、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、ｃｊ（ｄｋ）がｊ番目の固有モード関数成分であり、ｒｎ（ｄｋ）が分解した残量であるステップＳ３２と、
それぞれ各サブ時間特徴モードｖｉｋ（ｄｋ）に対して、三次スプライン補間法によりすべての固有モード関数成分ｃｊ（ｄｋ）及び残量信号ｒｎ（ｄｋ）におけるデータ欠損日付位置に対して補間補充を行い、次にすべての補間補充後の固有モード関数成分及び残量信号を改めて重ね合わせて、時系列が完全なサブ時間特徴モードｖｉｋ（Ｄｋ）を取得し、Ｄｋにはｋ番目の時点のすべての連続日付が含まれるステップＳ３３と、
同じ列の特徴モードベクトルのK個のサブ時間特徴モードｖｉｋ（Ｄｋ）を改めて重ね合わせてマージして、すべての時点を含む完全特徴モードベクトルｖｉ（Ｔ）を取得するステップＳ３４と、
すべての
2021533427.tif 000099
列の完全特徴モードベクトルｖｉ（Ｔ）を改めて重ね合わせてマージして、時系列的に連続した完全時間特徴モードＶＴを取得するステップＳ３５とを含むことを特徴とする請求項１に記載の再構築方法。
【請求項９】
前記静止軌道海色衛星リモートセンシング画像の日内データの頻度総数Kが８であることを特徴とする請求項８に記載の再構築方法。
【請求項１０】
前記ステップＳ４の具体的な方法は、
Ｓ２において取得された空間特徴モード
2021533427.tif 000100
及びＳ３において取得された完全時間特徴モードＶＴを利用して、行列Ｘにおけるすべてのデータ欠損点の再構築値を改めて算出し、行列Ｘの第ｉ行第ｊ列の欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数１２】
2021533427.tif 000101
式中、
2021533427.tif 000102
が空間特徴モード
2021533427.tif 000103
の第ｐ列のｉ番目の要素であり、
2021533427.tif 000104
が時間特徴モードＶＴの第ｐ列のｊ番目の要素であり、
2021533427.tif 000105
が対応する特異値であるステップＳ４１と、
Ｓ４１において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、行列Ｘの残りの要素値が変化せず、欠損値を補完した海色リモートセンシングデータセットの時空間フィールドを取得し、更に時間と空間が完全にカバーされた静止海色衛星データセットを取得するステップＳ４２とを含むことを特徴とする請求項１に記載の再構築方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は海洋リモートセンシング分野に属し、具体的に経験的直交関数分解法による静止海色衛星データの再構築方法に関する。
【背景技術】
【０００２】
海洋リモートセンシング衛星は全地球の海洋を広範囲且つ長期間に観測することができ、海洋を理解・認識するために他の観測方式が代替できないデータソースを人間に提供する。静止軌道海色衛星の出現は、海洋リモートセンシングについての研究に、海色の変化を時々モニタリングするという真新しい視角を提供する。世界の１番目の静止海色センサーであるＧＯＣＩ（ＧＯＣＩ、Ｇｅｏｓｔａｔｉｏｎａｒｙ  Ｏｃｅａｎ  Ｃｏｌｏｒ  Ｉｍａｇｅｒ）は初めて時間分解能が時間レベルに達する海色データを提供し、これにより海洋環境の時間ごとの変化をモニタリングすることが可能になる。ところが、海洋上空の雲、霧及び煙霧の影響のため、データに連続した高欠損率ひいては完全欠損がある場合が多く、データの時間的・空間的連続性に悪影響を与え、逆解析した海色製品の利用可能性を大幅に低下させ、その価値を十分に発揮することができない。
【０００３】
データの欠損は海洋リモートセンシングデータのよくある問題であり、欠損したリモートセンシングデータを補完するために、国内外ではデータの再構築を行う多くのデータ処理方法、例えば最適補間法、クリギング補間法、経験的直交関数分解法等を提供したが、経験的直交関数分解法（ＤＩＮＥＯＦ、Ｄａｔａ  ＩＮｔｅｒｐｏｌａｔｉｎｇ  Ｅｍｐｉｒｉｃａｌ  Ｏｒｔｈｏｇｏｎａｌ  Ｆｕｎｃｔｉｏｎｓ）は現在より効果的な海洋リモートセンシングデータの再構築方法であり、様々な海洋リモートセンシングデータ、例えば海表面温度、葉緑素ａ濃度、総浮遊物濃度、海表面塩度等の再構築に応用されている。経験的直交関数分解法は、その自身の制限のため、静止海色衛星データの再構築に多くの適合性問題があり、Ｌｉｕ  ＸｉａｏｍｉｎｇとＷａｎｇ  Ｍｅｎｇｈｕａが実施した再構築調査では、その時間幅が四半期より短く、且つデータ欠損率が比較的低い夏と秋の期間に集中している。
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本発明の目的は、従来技術の欠点を解決し、経験的直交関数分解法による静止海色衛星データの再構築方法を提供することである。
【課題を解決するための手段】
【０００５】
本発明の目的を実現するために、以下の技術案を提供する。   
【０００６】
経験的直交関数分解法による静止海色衛星データの再構築方法であって、
再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像を前処理して、異常ピクセルをデータ欠損点とするとともに、データセットにおける時間的又は空間的な欠損率が上限より高い極高欠損率のデータを除去するステップＳ１と、
前処理後の海色リモートセンシングデータセットの時空間フィールドで二次元データ行列を構築し、経験的直交関数の分解と再構築を行って、各反復分解前にラプラス平滑化フィルタリングにより時間共分散行列をフィルタリングし、再構築後に、最適な保留モード数と、対応する時間特徴モード及び空間特徴モードとを取得するステップＳ２と、
前記時間特徴モードを列によって特徴モードベクトルに分解し、各列の特徴モードベクトルをデータの属する時点によって複数のサブ時間特徴モードに分解し、さらに各サブ時間特徴モードに対して経験的モード分解を行い、三次スプライン補間法により分解後の固有モード関数成分及び残量における極高欠損率のデータに対応する箇所のデータを補完し、改めて重ね合わせてマージして時系列的に連続した完全時間特徴モードを取得するステップＳ３と、
Ｓ２で取得した空間特徴モード及びＳ３で取得した完全時間特徴モードを利用して、時間と空間が完全にカバーされた静止海色衛星データセットを再構築して取得するステップＳ４とを含む、ことを特徴とする再構築方法。
【０００７】
上記技術案に基づき、各ステップは以下の具体的な態様で実現されることが好ましい。ただし、各ステップの好ましい実現態様は、衝突しない限り、いずれも対応して組み合わせられてもよく、制限にならない。
【０００８】
好ましくは、前記ステップＳ１において、それぞれの静止軌道海色衛星リモートセンシング画像に対して、
再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像に対して、それぞれ画像における各有効ピクセルを中心として辺長Ｌの空間窓を設定するステップＳ１１、
各有効ピクセルのそれぞれに対して、該ピクセルに対応する空間窓におけるすべての有効ピクセルをサンプルデータセットＯとして、該ピクセルの集中度Ｏｃｏｎｃを算出し、集中度Ｏｃｏｎｃは、下式の通りであり、
【数１】
2021533427.tif 000003
式中、Ｏｉが集中度を算出べき現在ピクセル値であり、ＯｍがサンプルデータセットＯの中央値であり、標準偏差δの計算式がδ=ｋ*ＭＡＤ(Ｏ)であり、ｋが比例因子定数であり、ＭＡＤ（Ｏ）がサンプルデータセットＯの中央値絶対偏差であるステップＳ１２、
各有効ピクセルのそれぞれに対して、該ピクセルのエッジ検出値Ｏｐｒｏｘを算出し、該ピクセルの周囲に無効ピクセルがあれば、Ｏｐｒｏｘが正数値を取り、該ピクセルの周囲に無効ピクセルがなければ、Ｏｐｒｏｘがゼロを取るステップＳ１３、
各有効ピクセルのそれぞれに対して、該ピクセルの異常検出結果Ｏfｉｎａｌを算出し、異常検出結果Ｏfｉｎａｌは、下式の通りであり、
【数２】
2021533427.tif 000004
式中、wｐｒｏｘがエッジ検出重みであり、wｃｏｎｃが集中度検出重みであり、且つwｐｒｏｘ + wｃｏｎｃ =１を満足するステップＳ１４、
予め設定された異常ピクセル検出閾値Ｏｔｈｒｅｓｈｏｌｄに基づき、Ｏfｉｎａｌが閾値Ｏｔｈｒｅｓｈｏｌｄより大きいピクセルを異常ピクセルとして、異常ピクセルをデータ欠損点とするステップＳ１５に従って異常ピクセル識別処理を行う。
【０００９】
好ましくは、前記ステップＳ１において、前記極高欠損率のデータの除去方法は、
時間的に、いずれか１つの時点の静止軌道海色衛星リモートセンシング画像の欠損率が欠損率の上限より大きい場合、該時点に対応する画像を除去し、空間的に、すべての静止軌道海色衛星リモートセンシング画像におけるいずれか１つの位置点の時系列データの欠損率が欠損率の上限より大きい場合、該位置点のすべての画像でのピクセル値を除去する。
【００１０】
より好ましくは、前記欠損率の上限が９５％である。
【００１１】
好ましくは、前記ステップＳ２の具体的な方法は、
前処理後の海色リモートセンシングデータセットの時空間フィールドをｍ×ｎの二次元データ行列Ｘ０とし、ｍが静止軌道海色衛星リモートセンシング画像における空間位置点の個数であり、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、且つｍ>ｎであり、行列Ｘ０のｍ行がそれぞれｍ個の異なる空間位置点の時系列値を代表し、ｎ列がそれぞれｎ個の異なる時点のすべての空間位置点の値を代表するステップＳ２１と、
行列Ｘ０におけるすべての有効データの平均値
2021533427.tif 000005
を算出し、
2021533427.tif 000006
にし、行列Ｘから一部の有効データを交差検証セットＸｃとしてランダムに抽出し、行列Ｘにおけるすべてのデータ欠損点及び交差検証セットＸｃに抽出された要素にいずれも０の値を割り当てるステップＳ２２と、
保留モード数Ｐを反復して、異なるＰ値における空間特徴モードＵＰ、時間特徴モードＶＰ及び交差検証セットの再構築値と元の値との二乗平均平方根誤差Ｒを算出し、反復するたびにいずれものステップＳ２３１〜Ｓ２３４を実行するステップＳ２３と、
保留モード数Ｐ=１, ２,…, ｋＭＡＸにおける二乗平均平方根誤差ＲＰを比較し、ＲＰに対応する最小Ｐ値を最適な保留モード数`
2021533427.tif 000007
として、最適な時間特徴モード
2021533427.tif 000008
及び最適な空間特徴モード
2021533427.tif 000009
を取得するステップＳ２４とを含み、
ステップＳ２３１〜Ｓ２３４は以下のとおりであり、
ステップＳ２３１：行列Ｘに対して特異値分解を行い、特徴ベクトル分解によってＰ個の最大特異値及び特異ベクトルを算出し、更にｍ×Ｐの空間特徴モードＵＰ、Ｐ×Ｐの特異値行列SＰ及びｎ×Ｐの時間特徴モードＶＰを取得し、計算式は、下式の通りであり、
【数３】
2021533427.tif 000010
【数４】
2021533427.tif 000011
式中、uｉが空間特徴モードＵＰの第ｉ列であり、ｖｉが時間特徴モードＶＰの第ｉ列であり、ρｉがuｉ及びｖｉに対応する特異値であり、ｉ = １,..., Ｐであり、ＸＴ Ｘがｎ×ｎの時間共分散行列であり、
ステップＳ２３２：ラプラス平滑化フィルタリングによって前記時間共分散行列の列をフィルタリングし、次に前記時間共分散行列の行をフィルタリングし、時間共分散行列の行又は列がいずれもｎ個の時点データを含む時系列データセットＣ={ｃ１,ｃ２,…，ｃｎ}であると見なされ、Ｃに対応する時間セットがＴ={ｔ１,ｔ２,…,ｔｎ}であり、ｃｌがｌ番目の時点ｔｌの時系列データであり、ｌ=１,２,…,ｎであり、且つ
2021533427.tif 000012
を満足し、ラプラス平滑化フィルタリングの式は、下式の通りであり、
【数５】
2021533427.tif 000013
【数６】
2021533427.tif 000014
【数７】
2021533427.tif 000015
式中、ｃｌ(ｓ)がラプラス平滑化フィルタリング前のデータｃｌを示し、ｃｌ(ｓ+１)がラプラス平滑化フィルタリング後のデータｃｌを示し、αが平滑化強度であり、０＜α＜１であり、
ステップＳ２３３：１回のラプラス平滑化フィルタリングを完了した後、行列Ｘにおけるすべてのデータ欠損点の再構築値を算出し、行列Ｘの第ｉ行第ｊ列の欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数８】
2021533427.tif 000016
式中、(uｐ)ｉが空間特徴モードＵＰの第ｐ列のｉ番目の要素であり、(ｖｐＴ)ｊが時間特徴モードＶＰの第ｐ列のｊ番目の要素であり、ρｐが対応する特異値であり、
ステップＳ２３４：前記交差検証セットにおけるすべての要素の再構築値と真理値との二乗平均平方根誤差Ｒを算出し、二乗平均平方根誤差Ｒは、下式の通りであり、
【数９】
2021533427.tif 000017
式中、Ｎが前記交差検証セットにおけるデータの総個数であり、Ｘｒｅｐ、Ｘｃｐがそれぞれ交差検証セットにおけるｐ番目の要素の再構築値、真理値であり、
二乗平均平方根誤差Ｒが誤差閾値より大きい場合、Ｓ２３３において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、次に更新後の行列ＸでステップＳ２３１〜Ｓ２３３を繰り返し、二乗平均平方根誤差Ｒが閾値以下に収束するまで、繰り返しを停止して、収束された二乗平均平方根誤差Ｒを現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記す。
【００１２】
好ましくは、前記ステップＳ２３において、保留モード数Ｐの反復範囲が［１,ｋＭＡＸ］であり、最大反復回数がｋＭＡＸ≦ｍｉｎ（ｍ,ｎ）である。
【００１３】
より好ましくは、前記ステップＳ２３において、ステップＳ２３１〜Ｓ２３３を繰り返して特異値分解及びラプラス平滑化フィルタリングを行う回数が最大回数に達すると、繰り返しを停止して、最終的な二乗平均平方根誤差Ｒを現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記す。
【００１４】
好ましくは、前記ステップＳ３の具体的な方法は、
最適な時間特徴モード
2021533427.tif 000018
における各列の特徴モードベクトルを分解し、分解時に日を周期として異なる日における同じ時点に属するデータを１つのサブ時間特徴モードに分解し、第ｉ列の特徴モードベクトルｖｉ（ｔ）は、次のように分解され、
【数１０】
2021533427.tif 000019
式中、Ｋが静止軌道海色衛星リモートセンシング画像の日内データの頻度総数であり、ｖｉｋ（ｄｋ）がｖｉ（ｔ）におけるｋ番目の時点に対応するサブ時間特徴モードであり、ｖｉ（ｔ）におけるすべての日付のｋ番目の時点のデータを含み、
2021533427.tif 000020
であり、ｋ=１,２,…, Kであり、ｄｋがｖｉｋ（ｄｋ）におけるデータに含まれる日付を示し、ｄｋに前記極高欠損率のデータに対応する欠損日付を含まないステップＳ３１と、
各列の特徴モードベクトルにおける各サブ時間特徴モードｖｉｋ（ｄｋ）に対して経験的分解を行い、
【数１１】
2021533427.tif 000021
式中、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、ｃｊ（ｄｋ）がｊ番目の固有モード関数成分であり、ｒｎ（ｄｋ）が分解した残量であるステップＳ３２と、
それぞれ各サブ時間特徴モードｖｉｋ（ｄｋ）に対して、三次スプライン補間法によりすべての固有モード関数成分ｃｊ（ｄｋ）及び残量信号ｒｎ（ｄｋ）におけるデータ欠損日付位置に対して補間補充を行い、次にすべての補間補充後の固有モード関数成分及び残量信号を改めて重ね合わせて、時系列が完全なサブ時間特徴モードｖｉｋ（Ｄｋ）を取得し、Ｄｋにはｋ番目の時点のすべての連続日付が含まれるステップＳ３３と、
同じ列の特徴モードベクトルのK個のサブ時間特徴モードｖｉｋ（Ｄｋ）を改めて重ね合わせてマージして、すべての時点を含む完全特徴モードベクトルｖｉ（Ｔ）を取得し、Ｔが極高欠損率のデータを含むデータに対応する時点を示すステップＳ３４と、
すべての`
2021533427.tif 000022
列の完全特徴モードベクトルｖｉ（Ｔ）を改めて重ね合わせてマージして、時系列的に連続した完全時間特徴モードＶＴを取得するステップＳ３５とを含む。
【００１５】
より好ましくは、前記静止軌道海色衛星リモートセンシング画像の日内データの頻度総数Ｋが８である。
【００１６】
好ましくは、前記ステップＳ４の具体的な方法は、
Ｓ２において取得された空間特徴モード
2021533427.tif 000023
及びＳ３において取得された完全時間特徴モードＶＴを利用して、行列Ｘにおけるすべてのデータ欠損点の再構築値を改めて算出し、行列Ｘの第ｉ行第ｊ列の欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数１２】
2021533427.tif 000024
式中、
2021533427.tif 000025
が空間特徴モード
2021533427.tif 000026
の第ｐ列のｉ番目の要素であり、
2021533427.tif 000027
が時間特徴モードＶＴの第ｐ列のｊ番目の要素であり、
2021533427.tif 000028
が対応する特異値であるステップＳ４１と、
Ｓ４１において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、行列Ｘの残りの要素値が変化せず、欠損値を補完した海色リモートセンシングデータセットの時空間フィールドを取得し、更に時間と空間が完全にカバーされた静止海色衛星データセットを取得するステップＳ４２とを含む。
【発明の効果】
【００１７】
従来技術に比べて、本発明は以下の有益な効果を有する。
本発明は、ＤＩＮＥＯＦ方法の制限に対して、高時間分解能を有する特徴及びデータ欠損に連続した高欠損率と完全欠損があるという状況を十分に考慮し、まず集中度検出及びエッジ検出によってリモートセンシングデータにおける異常ピクセルを検査することにより、異常値がＤＩＮＥＯＦ方法における交差検証に干渉することを回避し、次にラプラス平滑化フィルタリングによってＤＩＮＥＯＦ分解における時間共分散行列をフィルタリングすることにより、夜間間隔及びデータ連続欠損によるデータギャップを減少させ、次に再構築により取得された時間特徴モードを二次分解してから補間することにより、高欠損率のデータ及び連続欠損データの再構築を実現し、最終的にその元の精度を維持し時間と空間が完全にカバーされた水色リモートセンシング再分析データセットを取得する。本発明は、保留モード数（Ｍｏｄｅ）、再構築精度（ＭＳＥ、ＲＭＳＥ、ＭＡＥ）、信号解釈率及び再構築率の面でいずれも改良され、操作しやすく、再構築結果が直感的である。本発明は、他の、高い時間分解能特徴を有するリモートセンシングデータ又は時間間隔の変化がより大きい水色衛星データの再構築に模範及び参考を提供する。
【図面の簡単な説明】
【００１８】
【図１】経験的直交関数分解法による静止海色衛星データの再構築方法のフローチャートである。
【図２】２０１７年６月２日のＧＯＣＩによる総浮遊物濃度の元の画像と再構築された画像との比較図である。
【図３】２０１７年６月２日のＧＯＣＩによる総浮遊物濃度の再構築精度結果を示す図である。
【発明を実施するための形態】
【００１９】
以下、図面を参照しながら具体的な実施形態によって本発明を更に説明する。
【００２０】
図１は本発明の好適な実施例の経験的直交関数分解法による静止海色衛星データの再構築方法のフローチャートである。該方法は、主に、
再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像（ＧＯＣＩ画像）を前処理して、異常ピクセルを識別して異常ピクセルをデータ欠損点とするとともに、データセットにおける時間的又は空間的な欠損率が上限より高い極高欠損率のデータを除去するＳ１と、
前処理後の海色リモートセンシングデータセットの時空間フィールドで二次元データ行列を構築し、経験的直交関数の分解と再構築を行って、各反復分解前にラプラス平滑化フィルタリングにより時間共分散行列をフィルタリングし、再構築後に、最適な保留モード数と、対応する時間特徴モード及び空間特徴モードとを取得するＳ２と、
Ｓ２で取得された時間特徴モードを列によって特徴モードベクトルに分解し、各列の特徴モードベクトルをデータの属する時点によって複数のサブ時間特徴モードに分解し、さらに各サブ時間特徴モードに対して経験的モード分解を行い、三次スプライン補間法により分解後の固有モード関数成分及び残量における極高欠損率のデータに対応する箇所のデータを補完し、改めて重ね合わせてマージして時系列的に連続した完全時間特徴モードを取得するＳ３と、
Ｓ２で取得された空間特徴モード及びＳ３で取得された完全時間特徴モードを利用して、再構築によって時間と空間が完全にカバーされた静止海色衛星データセットを取得するＳ４と、を含む。
【００２１】
以下、本実施例におけるＳ１〜Ｓ４の具体的な実現方式及びその効果を詳しく説明する。
【００２２】
まず、ＧＯＣＩデータはデータ品質制御を受けていますが、依然として合理的な値域範囲内に異常値がある可能性がある。このため、ステップＳ１の目的の１つはこれらの異常ピクセルを処理することにより、再構築過程にノイズが生じることを防止することである。一般的な集中度検出方式は、すべてのデータの平均値及び標準偏差に基づいてデータが「異常値」であるかどうかを判断する。該方法は、算出しやすいが、検出対象とするデータセットに多くの異常値があれば、平均値及び標準偏差によるデータの全体特徴の説明に大きく影響し、それにより異常値の検出に不確定性がより高い。上記状況を回避するために、中央値及び中央値絶対偏差（ＭＡＤ、Ｍｅｄｉａｎ  Ａｂｓｏｌｕｔｅ  Ｄｅｖｉａｔｉｏｎ）で平均値及び標準偏差を置き換えることはより確実な検出方式である。
【００２３】
本実施例では、それぞれの静止軌道海色衛星リモートセンシング画像に対して、前処理過程において集中度検出及びエッジ検出の重み付け結果を基準として、Ｓ１１〜Ｓ１５によって異常ピクセルの識別及び処理を行う。
Ｓ１１：再構築対象とする、元の海色リモートセンシングデータセットにおける各静止軌道海色衛星リモートセンシング画像に対して、それぞれ画像における各有効ピクセルを中心として辺長Lの空間窓を設定する。Lは、実際の状況に応じて調整してもよい。
Ｓ１２：各有効ピクセルのそれぞれに対して、該ピクセルに対応する空間窓におけるすべての有効ピクセルをサンプルデータセットＯとして、該ピクセルの集中度Ｏｃｏｎｃを算出し、集中度Ｏｃｏｎｃは、下式の通りであり、
【数１３】
2021533427.tif 000029
式中、Ｏｉが集中度を算出べき現在のピクセル値であり、ＯｍがサンプルデータセットＯの中央値であり、標準偏差δの計算式がδ=ｋ*ＭＡＤ(Ｏ)であり、ｋが比例因子定数であり、ＭＡＤ（Ｏ）がサンプルデータセットＯの中央値絶対偏差である。
正規分布データについては、正常値が中間の５０％の領域にあると仮定すれば、異常値が両側の５０％の領域にあり、換算によってｋ＝１．４８２６を取得できる。３σ基準に従い、Ｏｃｏｎｃが３より大きいピクセルは異常ピクセルであると見なされてもよい。
Ｓ１３：本発明は、異常ピクセルモニタリングの安定性を強化するために、異常ピクセルの発生要因、つまり異常ピクセルの雲層エッジでの発生確率がより高いことに基づき、簡単なエッジ検出方法で雲層エッジ位置のピクセルを記す。つまり、各有効ピクセルのそれぞれに対して該ピクセルのエッジ検出値Ｏｐｒｏｘを算出し、該ピクセルの周囲の８つのピクセルに無効ピクセルがあれば、Ｏｐｒｏｘが正数値Ｍを取り、該ピクセルの周囲の８つのピクセルに無効ピクセルがなければ、Ｏｐｒｏｘがゼロを取る。本実施例では、下記の通り、集中度検出結果の数値と一致するように、ピクセルの周囲に無効ピクセルがあれば、Ｍ＝３と記し、そうでなければ、０と記す。
【数１４】
2021533427.tif 000030
Ｓ１４：各有効ピクセルのそれぞれに対して、該ピクセルの異常検出結果Ｏfｉｎａｌを算出し、異常検出結果Ｏfｉｎａｌは、下式の通りであり、
【数１５】
2021533427.tif 000031
式中、wｐｒｏｘがエッジ検出重みであり、wｃｏｎｃが集中度検出重みであり、且つwｐｒｏｘ + wｃｏｎｃ =１を満足する。wｐｒｏｘ、wｃｏｎｃは、実際の状況に応じて調整してもよい。
重み付け結果を使用して異常ピクセルを検出することにより、検出安定性を向上させ、単一の方法の過度懲罰によるデータ欠損がより多く、データの元の変化特徴を失うことを回避する一方、更に所定の原データの異常ピクセルのソース及び特徴に対して重みを調整して異常値の検出効果を向上させることができる。
Ｓ１５：予め設定された異常ピクセル検出閾値Ｏｔｈｒｅｓｈｏｌｄに基づき、Ｏfｉｎａｌ が閾値Ｏｔｈｒｅｓｈｏｌｄより大きいピクセルを異常ピクセルとして、異常ピクセルをデータ欠損点とする。閾値Ｏｔｈｒｅｓｈｏｌｄの大きさは除去された異常ピクセルの数を決定し、本実施例では検出閾値Ｏｔｈｒｅｓｈｏｌｄを３として設定する。
【００２４】
また、原データの時間的次元又は空間的次元での欠損率が高すぎる（一般的に９５％を超える）と、これらのデータは有効情報を提供できないだけでなく、最終結果に影響する恐れもある。従って、経験的直交関数分解法（以下、ＤＩＮＥＯＦと略称される）は、これらの極高欠損率のデータを再構築できず、再構築前に除去する必要があり、そうでなければ、高欠損率画像の再構築結果は元のリモートセンシングデータセットの平均値で代替される。
【００２５】
本実施例では、極高欠損率のデータの除去方法は、時間的に、いずれか１つの時点の静止軌道海色衛星リモートセンシング画像の欠損率が欠損率の上限より大きい場合、該時点に対応する画像を除去し、空間的に、すべての静止軌道海色衛星リモートセンシング画像におけるいずれか１つの位置点の時系列データ欠損率が欠損率の上限より大きい場合、該位置点のすべての画像でのピクセル値を除去し、除去されたピクセルをデータ欠損点としない。一般的に欠損率の上限を９５％として設定してもよい。
【００２６】
ＤＩＮＥＯＦ方法で再構築された原データは時間的に連続し、且つデータの間隔時間が同じであることを保持すべきである。静止海色衛星データについては、そのデータ自体は当日で連続したものであって時間間隔が一致するが、高頻度の昼間データの後で、長期間にデータがない夜間隙間があり、夜間隙間が一般的に１２時間を超え、且つ冬でより長くなる。一方、雲層の被覆時間は一般的に静止海色衛星のモニタリング周期より長く、データには領域データの連続した高欠損率ひいては完全欠損が生じる場合がより多く、このため、データの連続性を更に低下させる。連続性の悪いデータには、いくつかのデータの非正常な変化特徴があるが、これらの特徴はデータ全体の時空間変動性を代表できず、このため、ＤＩＮＥＯＦ方法は少量の原データと非正常な変化とを過度にフィッティングし、結果として再構築画像が必ずしも該所定の画像の有すべき時空間変化を反映するとは限らず、再構築データに明らかなギャップがあり、非常に不合理な連続時点画像が生成してしまう。本発明はフィルタリングにより元のリモートセンシングデータを平滑化することを考慮し、一般的な平滑化方式（例えば、中央値フィルタリング、平均値フィルタリング等）と異なり、該フィルタリングは突然変化及びノイズをなくすとともに、データの元の変化特徴を維持する必要がある。これに対して、本発明はラプラス平滑化方法を用いて時系列データの平滑化に適するフィルタリングを構築して、ＤＩＮＥＯＦの再構築過程に応用される。以下、まず該フィルタリング過程の理論を説明する。
【００２７】
データ行列Ｘ（ｍ×ｎ）のｍ個の時系列にラプラスフィルタリングを用いることは、一次元時系列に作用するラプラス平滑化フィルタリング（ｎ×１）を二次元に拡張すると理解されてもよく、構築されたラプラスフィルタリングＬがｎ×ｎの行列であり、行列Ｘを時間域においてフィルタリングし、
【数１６】
2021533427.tif 000032
2021533427.tif 000033
がラプラス平滑化フィルタリング後の元の行列である。理論的に該フィルタリングによりＸを直接フィルタリングしてもよいが、元の行列に対する処理計算量が大きすぎて、ｍ回の計算が必要である。
2021533427.tif 000034
は、特異値分解（ＳＶＤ）の計算において、最も重要な最初の数個の時間特徴モードと空間特徴モードのみを必要とし、特徴ベクトル分解によってＰ個の最大特異値及び特異ベクトルを算出することを考慮して、方程式は、下式の通りであり、
【数１７】
2021533427.tif 000035
ｖｉが時間特徴モードＶの第ｉ列であり、ρｉが対応する特異値であり、ｉ＝１，...，Ｐである。
2021533427.tif 000036
がｎ×ｎの行列であり、時間共分散行列Ａとして定義する。該フィルタリングＬを
2021533427.tif 000037
のＳＶＤ分解と組み合わせて、元の行列に対するフィルタリングを時間共分散行列に対するフィルタリングで代替する。ラプラス平滑化フィルタリングによりフィルタリングした行列
2021533427.tif 000038
に対して特徴ベクトル分解表現を行うことにより、ラプラス平滑化フィルタリング後の時間共分散行列を把握でき、
【数１８】
2021533427.tif 000039
上記式によれば、まずにラプラスフィルタリングＬにより時間共分散行列Ａの列をフィルタリングし、次に時間共分散行列Ａの行をフィルタリングする。これにより、算出量を大幅に減少させる（ｍ回から２ｎ回に減少させる）とともに、欠損データへの感度も低下させる（時間共分散行列に欠損値がない）。ＤＩＮＥＯＦ分解された二乗平均平方根誤差によって最適なαとＰの値を決定することができる。一般的にα値が一定を維持し、Ｐ値を最適化すればよい。
【００２８】
上記理論に基づき、本実施例の経験的直交関数分解再構築過程は、主にステップＳ２により実現され、以下にその具体的な方法を詳しく説明する。
【００２９】
Ｓ２１：前処理後の海色リモートセンシングデータセットの時空間フィールドをｍ×ｎの二次元データ行列Ｘ０とし、ｍが静止軌道海色衛星リモートセンシング画像における空間位置点の個数であり、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、且つｍ＞ｎであり、行列Ｘ０の各行がそれぞれ１つの空間位置点の異なる時点での画素値からなる時系列を代表し、各列がそれぞれ１つの時点のすべての空間位置点の値、すなわちある時点の１幅のＧＯＣＩ画像における画素値を代表する。従って、行列Ｘ０のｍ行がそれぞれｍ個の異なる空間位置点の時系列値を代表し、ｎ列がそれぞれｎ個の異なる時点のすべての空間位置点の値を代表する。Iをデータ欠損点セット（元に欠損したもの及びＳ１で識別された異常ピクセルを含む）と記し、行列における欠損値がＮａＮで示され、Ｐを保留モード数と記し、Ｐ≦ｍｉｎ（ｍ,ｎ）である。
【００３０】
Ｓ２２：行列Ｘ０におけるすべての有効データの平均値
2021533427.tif 000040
を算出し、
2021533427.tif 000041
にし、Ｘ０におけるすべての有効要素の値から平均値
2021533427.tif 000042
を引くことを示し、行列Ｘから一部の有効データを交差検証セットＸｃとしてランダムに抽出し、一般的に有効データ総量の１％とし、最適な保留モード数の判断基準とする。Ｘにおける交差検証セットの同じ位置にある点も欠損点と見なし、値がＮａＮである。行列Ｘにおけるすべてのデータ欠損点及び交差検証セットＸｃに抽出された要素にいずれも０の値を割り当て、つまり０でＸにおけるすべての値がＮａＮである点を代替する。
【００３１】
Ｓ２３：保留モード数Ｐを反復して、異なるＰ値における空間特徴モードＵＰ、時間特徴モードＶＰ及び交差検証セットの再構築値と元の値との二乗平均平方根誤差Ｒを算出する。Ｐが正の自然数であり、［１，ｋＭＡＸ］の値区間内において反復を始める前にＰ＝１に初期化し、次に１回反復するたびにＰ＝Ｐ＋１にし、最大反復回数ｋＭＡＸ≦ｍｉｎ（ｍ,ｎ）である。
【００３２】
反復するたびにいずれもステップＳ２３１〜Ｓ２３４を実行する。
【００３３】
Ｓ２３１：行列Ｘに対して特異値分解を行い、特徴ベクトル分解によってＰ個の最大特異値及び特異ベクトルを算出し、更にｍ×Ｐの空間特徴モードＵＰ、Ｐ×Ｐの特異値行列SＰ及びｎ×Ｐの時間特徴モードＶＰを取得し、計算式は、下式の通りであり、
【数１９】
2021533427.tif 000043
【数２０】
2021533427.tif 000044
式中、uｉが空間特徴モードＵＰの第ｉ列であり、ｖｉが時間特徴モードＶＰの第ｉ列であり、ρｉがuｉ及びｖｉに対応する特異値であり、ｉ= １,...,Ｐであり、ＸＴＸがｎ×ｎの時間共分散行列であり、上付き文字Ｔが転置を代表する（以下同様）。
これにより、現在の保留モード数Ｐにおける空間特徴モードＵＰ、時間特徴モードＶＰ及び特異値行列SＰを取得することができ、後続の呼び出しに備える。
【００３４】
Ｓ２３２：ラプラス平滑化フィルタリングによって前記時間共分散行列の列をフィルタリングし、次に前記時間共分散行列の行をフィルタリングする。時間共分散行列の行又は列をフィルタリングするとき、いずれもｎ個の時点データを含む時系列データセットＣ={ｃ１,ｃ２,…，ｃｎ}と見なし、Ｃに対応する時間セットがＴ={ｔ１,ｔ２,…, ｔｎ}であり、ｃｌがｌ番目の時点ｔｌの時系列データであり、ｌ=１,２,…, ｎであり、且つ
2021533427.tif 000045
を満足する。ラプラス平滑化フィルタリングの式は、下式の通りであり、
【数２１】
2021533427.tif 000046
【数２２】
2021533427.tif 000047
【数２３】
2021533427.tif 000048
式中、ｃｌ(ｓ)がラプラス平滑化フィルタリング前のデータｃｌを示し、ｃｌ(ｓ+１)がラプラス平滑化フィルタリング後のデータｃｌを示し、αが平滑化強度であり、０＜α＜１である。上付き文字ｓがｓ回目のラプラス平滑化フィルタリングを示し、ラプラス平滑化フィルタリング回数ｓの値が１、２・・・ｓｍａｘである。ｓの初期値が１であり、フィルタリングするたびにｓ＝ｓ＋１にする。
フィルタリングによる平滑化効果はα及び最大反復回数ｓｍａｘの値によって決定される。該フィルタリングは、時系列変化特徴を維持するとともに時系列データの平滑化を実現し、時間相対性の影響を考慮する。
【００３５】
Ｓ２３３：ラプラス平滑化フィルタリングを完了したたびに、行列Ｘにおけるすべてのデータ欠損点の再構築値を算出する。行列Ｘにおける行列Ｘの第ｉ行第ｊ列の要素が欠損点である場合、該欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数２４】
2021533427.tif 000049
式中、ｉ、ｊがＸの空間及び時間の下付き文字であり、(uｐ)ｉが空間特徴モードＵＰの第ｐ列のｉ番目の要素であり、（ｖｐＴ)ｊが時間特徴モードＶＰの第ｐ列のｊ番目の要素であり、ρｐが対応する特異値である。
これにより、行列Ｘにおけるすべてのデータ欠損点の再構築値を取得することができ、上記ステップにおいて行列Ｘから一部のデータを交差検証セットとして取り出したため、フィルタリング後に交差検証セットの真理値及び予測した再構築値に基づいて誤差を算出してもよい。
【００３６】
Ｓ２３４：前記交差検証セットにおけるすべての要素の再構築値と真理値との二乗平均平方根誤差Ｒを算出し、二乗平均平方根誤差Ｒは、下式の通りであり、
【数２５】
2021533427.tif 000050
式中、Ｎが前記交差検証セットにおけるデータの総個数であり、Ｘｒｅｐ、Ｘｃｐがそれぞれ交差検証セットにおけるｐ番目の要素の再構築値、真理値であり、
二乗平均平方根誤差Ｒが誤差閾値より大きい場合、Ｓ２３３において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、該更新過程において、データ欠損点のみを更新するが、他の位置が変化しないように維持する。従って、更新後の行列Ｘｒｅは、次のように表すことができる。
【数２６】
2021533427.tif 000051
式中、∂Ｘが欠損点の再構築値からなる行列である。
次に更新後のＸｒｅを改めて行列Ｘとし、ステップＳ２３１〜Ｓ２３３を繰り返す。繰り返すたびにいずれも二乗平均平方根誤差Ｒを算出する必要があり、二乗平均平方根誤差Ｒが閾値以下に収束するまで、繰り返しを停止して、収束された二乗平均平方根誤差Ｒを現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記す。しかしながら、算出効率を考慮し、繰り返し回数上限すなわち上記ｓｍａｘを設定する必要がある。ステップＳ２３１〜Ｓ２３３を繰り返して特異値分解及びラプラス平滑化フィルタリングを行う回数が最大回数ｓｍａｘに達すると、繰り返しを停止し、最終的な二乗平均平方根誤差Ｒを直接現在の保留モード数Ｐにおける二乗平均平方根誤差ＲＰと記す。
【００３７】
Ｓ２４：上記ステップＳ２３の実行中において、保留モード数Ｐの反復範囲が［１，ｋＭＡＸ］であり、異なる保留モード数Ｐ＝１, ２,…,ｋＭＡＸにおける二乗平均平方根誤差ＲＰを取得することができ、ＲＰを比較し、ＲＰに対応する最小Ｐ値を最適な保留モード数
2021533427.tif 000052
とし、最適な時間特徴モード
2021533427.tif 000053
及び最適な空間特徴モード
2021533427.tif 000054
を取得する。Ｐ＝ａの場合、二乗平均平方根誤差Ｒが最も小さく、従って最適な保留モード数`
2021533427.tif 000055
である。
【００３８】
上記Ｓ２ステップを経て、最終的に最適な保留モード数
2021533427.tif 000056
を決定し、最適な特異値行列
2021533427.tif 000057
、時間特徴モード
2021533427.tif 000058
及び最適な空間特徴モード
2021533427.tif 000059
を取得することができる。以上のように、ＤＩＮＥＯＦが原データの統計情報に依存して再構築し、欠損率が９５％より大きい極高欠損率のデータを正確に再構築できず、従って、再構築前に除去する必要があり、この部分の除去されたデータはＳ２ステップにおけるラプラスフィルタリング過程に参加せず、ステップＳ３に依存して補間再構築を行う必要がある。
【００３９】
ＤＩＮＥＯＦ方法により再構築して取得された時間特徴モードそのものは、データセットの時系列変化特徴を反映しており、且つ隣接する時点のデータは、同じモードに対応する時間特徴モード値において関連性と連続性を有しているため、時間特徴モードを補間して極高欠損率の時点に対応する時間特徴モードベクトルを取得することは実現可能な方法である。しかしながら、静止海色衛星データについては、該方法は一定の欠陥があり、極高欠損率の時点のデータの再構築に直接応用できない。その理由は、主に、静止海色衛星データの再構築により取得された時間特徴モードについては、その時間分解能標準での連続性が１日の中に限られるためである。つまり、異なる日の隣接するデータを利用して補間するために、時間特徴モードを時間分解能に対応して、不等間隔補間を行う必要があり、周期性を有する特徴モードデータについては、従来の方法による補間効果がより低い。一方、当日に雲や霧に覆われていた場合、１日の中のすべての時点データがいずれも連続した完全欠損データである。このような場合、隣接する数値を使用して時間間隔で補間すると、その補間結果が実際の状況に一致することが困難である。これらの問題に対して、本発明は、時間分解能を低下させるやり方を用い、高時間分解能リモートセンシングデータが低時間分解能において依然として良好な周期性を有するという特徴を利用して、再構築により取得された時間特徴モードをより低い時間分解能によって分解し、連続性がより高く時間間隔が一致する複数のサブ時間特徴モードを構築し、次に経験的モード分解方法を用いてサブ時間特徴モードを再び分解し、分解結果における欠損信号を三次スプライン法で補間補充し、最後に時間的に連続し間隔が一致するサブ時間特徴モードにマージして、更にすべての時点を含む新しい時間特徴モードを取得することを提案する。該方法は時間特徴モードを２回分解してから補間するため、本発明ではこれを「時間モード二次分解補間」と呼ぶ。該方法の趣旨は、静止海色衛星データのデータ取得時間が低時間分解能において等時間間隔であり、つまり連続して複数日の同じ時点に取得したデータが同じ時間間隔（２４時間）を有し、データが周期性を有する点である。一方、サブ時間特徴モードベクトルを直接補間するのではなく、経験的モード分解により取得された各固有モード関数成分を補間することを選択する点であり、その理由は、分解された信号の波動性が原信号より低く、補間補充に役立つとともに、固有モード関数成分にはより多くの元のベクトルの有しない詳細が含まれており、一層高い精度を実現することができるためである。
【００４０】
本発明では、説明しやすくするために、異なる日付の同じ時点（例えば、異なる日付の朝８時）に取得したデータが同じ時点のデータと称される。隣接する２日間の同じ時点のデータ間隔が一定、すなわち２４ｈである。
【００４１】
本実施例では、ステップＳ３において「時間モード二次分解補間」を行う具体的な方法は、次のとおりである。
Ｓ３１：最適な時間特徴モード
2021533427.tif 000060
`における各列の特徴モードベクトルを分解し、分解時に日を周期として異なる日における同じ時点に属するデータを１つのサブ時間特徴モードに分解し、第ｉ列の特徴モードベクトルｖｉ（ｔ）は、次のように分解され、
【数２７】
2021533427.tif 000061
式中、
2021533427.tif 000062
であり、ｋ=１,２,…, Kである。Ｋが静止軌道海色衛星リモートセンシング画像の日内データの頻度総数であり、本発明に用いられるＧＯＣＩデータは、毎日北京時間８：３０〜１５：３０に該領域を８回走査して、時間間隔が１時間である８幅の画像を取得し、従って、Ｋの値が８である。ｖｉｋ（ｄｋ）がｖｉ（ｔ）におけるｋ番目の時点に対応するサブ時間特徴モードであり、ｖｉ（ｔ）におけるすべての日付のｋ番目の時点のデータを含み、ｄｋがｖｉｋ（ｄｋ）におけるデータに含まれる日付を示し、且つデータにおける極高欠損率のデータが除去されたため、分解後のサブ時間特徴モードにおけるこれらのデータに対応する日付が欠損したものであり、従ってｄｋで示され、つまりｄｋには前記極高欠損率のデータに対応する欠損日付が含まれていないが、ｖｉ（ｔ）におけるｔも同様にデータ欠損を示すことに用いられる。後続のステップによりこれらの欠損日付位置のデータを補間補完する。
【００４２】
Ｓ３２：各列の特徴モードベクトルにおける各サブ時間特徴モードｖｉｋ（ｄｋ）に対して経験的分解を行い、
【数２８】
2021533427.tif 000063
式中、ｎが海色リモートセンシングデータセットにおける画像で代表される時点総数であり、ｃｊ（ｄｋ）がｊ番目の固有モード関数成分であり、ｒｎ（ｄｋ）が分解して残量である。
【００４３】
Ｓ３３：それぞれ各サブ時間特徴モードｖｉｋ（ｄｋ）に対して、三次スプライン補間法によりすべての固有モード関数成分ｃｊ（ｄｋ）及び残量信号ｒｎ（ｄｋ）におけるデータ欠損日付位置を補間補充し、次に逆順序で上記の経験的分解公式を用い、すべての補間補充後の固有モード関数成分及び残量信号を改めて重ね合わせて、時系列が完全なサブ時間特徴モードｖｉｋ（Ｄｋ）を取得し、この時のサブ時間特徴モードにおけるデータ欠損日付位置には、既に対応するデータを補間補充したため、Dkで示され、つまりDkにはｋ番目の時点のすべての連続日付が含まれる。
【００４４】
Ｓ３４：同じ列の特徴モードベクトルのＫ個のサブ時間特徴モードｖｉｋ（Ｄｋ）を改めて重ね合わせてマージして、すべての時点を含む完全特徴モードベクトルｖｉ（Ｔ）を取得し、Ｔが極高欠損率のデータを含むデータに対応する時点を示す。
【００４５】
Ｓ３５：すべての
2021533427.tif 000064
列の完全特徴モードベクトルｖｉ（Ｔ）を改めて順に重ね合わせて行列形式にマージして、時系列的に連続した完全時間特徴モードＶＴを取得する。
【００４６】
これにより、上記Ｓ１〜Ｓ３によって、完全空間特徴モード及び完全時間特徴モードを取得した。つまり、ＤＩＮＥＯＦ方法で時間と空間が完全にカバーされた静止海色衛星データセットを再構築して取得することができる。以下、ステップＳ４の具体的な再構築方法を簡単に説明する。
【００４７】
Ｓ４１：Ｓ２において取得された空間特徴モード
2021533427.tif 000065
及びＳ３において取得された完全時間特徴モードＶＴを利用して、行列Ｘにおけるすべてのデータ欠損点の再構築値を改めて算出し、行列Ｘの第ｉ行第ｊ列の欠損点の再構築値Ｘｒｅｉ,ｊは、下式の通りであり、
【数２９】
2021533427.tif 000066
式中、上付き文字Ｔがいずれも転置を示し、「極高欠損率のデータを含むデータに対応する時点Ｔ」と異なり、
2021533427.tif 000067
が空間特徴モード
2021533427.tif 000068
の第ｐ列のｉ番目の要素であり、
2021533427.tif 000069
が時間特徴モードＶＴの第ｐ列のｊ番目の要素であり、
2021533427.tif 000070
が
2021533427.tif 000071
及び
2021533427.tif 000072
に対応する特異値である。
【００４８】
Ｓ４２：Ｓ４１において取得された行列Ｘにおけるすべてのデータ欠損点の再構築値を行列Ｘに更新し、該更新過程において、データ欠損点のみを更新するが、行列Ｘの他の位置が変化しないように維持し、これにより取得された更新後の行列が欠損値を補完した海色リモートセンシングデータセットの時空間フィールドである。該更新後の行列によって、改めて分解して復元して、時間と空間が完全にカバーされた静止海色衛星データセットを取得することができ、該データセットにおけるすべての静止軌道海色衛星リモートセンシング画像が時間的且つ空間的にいずれも連続し完全なものである。
【００４９】
以下、上記実施例の方法に基づき、具体的な実例に応用してその効果を示す。具体的な過程は以上のとおりであり、ここで詳細な説明は省略し、以下に主にその具体的なパラメータの設定及び実現効果を示す。
応用例
【００５０】
以下に杭州湾で２０１７年のＧＯＣＩ総浮遊物濃度データを例として、本発明を具体的に説明し、その具体的なステップは以下のとおりである。
【００５１】
１）、韓国海洋衛星センターが提供するＧＯＣＩ  Ｌｅｖｅｌ−１Ｂデータを用い、紫外線大気校正方法（ＵＶ−ＡＣ）により大気校正を行い、経験的モデルを使用して総浮遊物濃度データを逆解析する。ＧＯＣＩデータをＴＳＭ逆解析モデルによって逆解析して、杭州湾２０１７年の総浮遊物（Ｔｏｔａｌ  Ｓｕｓｐｅｎｄｅｄ  Ｍａｔｔｅｒ、ＴＳＭ）の逆解析結果製品（ＧＯＣＩ  ＴＳＭ）を２９２０幅取得する。各ＧＯＣＩ  ＴＳＭデータ製品のカバレッジ範囲が１２０．５５°Ｅ〜１２２．３５°Ｅ、２９．８５°Ｎ〜３０．９５°Ｎであり、ピクセルサイズが０．００５°×０．００５°である。２９２０幅のＧＯＣＩ  ＴＳＭデータの総欠損率が６３．５％であり、空間的に、研究領域における有効位置点の時系列データの欠損率（欠損時点）が高くとも９９．９３％（２９１８幅）であり、低くとも５３．９４％（１５７５幅）である。時間的に、２９２０幅のＧＯＣＩ  ＴＳＭデータにおける空間的な欠損率が高くとも１００％であり、低くとも０．８９％であり、欠損率が９５％より大きいデータが１１２４幅あり、占有率が３８．４９％である。
【００５２】
２）上記ステップＳ１に基づき、２０１７年の２９２０幅のＧＯＣＩ総浮遊物濃度データについては、Ｌ＝９にし、つまり９×９の空間窓を用いて集中度検出を行い、９×９の空間窓を用いて集中度検出を行い、集中度検出重みwｃｏｎｃとエッジ検出重みwｐｒｏｘがいずれも１/ ２であり、検出閾値Ｏｔｈｒｅｓｈｏｌｄを３として設定し、異常ピクセルをデータ欠損点としてその中の原データを去除する。ＤＩＮＥＯＦの原データの欠損率への制限のため、更に画像欠損率が９５％より大きい対応する時点の原データ及び時系列欠損率が９５％より大きい対応する位置点のすべての画素を除去する。データ除去操作によって使用可能なＧＯＣＩ  ＴＳＭデータの数を２９２０から１７９４まで減少させ、即ち３８．５６％減少させる。海洋画素が３９２７４から３７１５２まで減少し、即ち５．４％減少し、主に除去した画素が主に沿岸水域に集中している。除去後に再構築データの総欠損率が３８．１８％である。ＤＩＮＥＯＦ内の算出基礎は原データの平均値及び共分散に基づくものであり、ＴＳＭデータが正規分布を呈さず、ＤＩＮＥＯＦの前に原データの変換を実行する必要がある。従って、本発明は原データに対して１０を底とする対数化処理を行う。
【００５３】
３）上記Ｓ２ステップに基づいて前処理後の海色リモートセンシングデータセットの時空間フィールドで二次元データ行列を構築し、経験的直交関数の分解と再構築を行って、各反復分解前にラプラス平滑化フィルタリングにより時間共分散行列をフィルタリングする。再構築後に、最適な保留モード数と、対応する時間特徴モード及び空間特徴モードとを取得する。本実施例では、αを０．０１とし、ｓｍａｘを３０とする。最終的に、Ｐ＝１６９の場合、二乗平均平方根誤差Ｒが最も小さく、従って最適な保留モード数
2021533427.tif 000073
`であり、該モード数における特異値行列
2021533427.tif 000074
、時間特徴モード
2021533427.tif 000075
及び最適な空間特徴モード
2021533427.tif 000076
を取得する。
【００５４】
４）上記Ｓ３ステップに基づき、時間特徴モードＶ`Ｐを列によって特徴モードベクトルに分解し、各列の特徴モードベクトルをデータの属する時点によって８つのサブ時間特徴モードに分解し、次に各サブ時間特徴モードに対して経験的モード分解を行い、三次スプライン補間法により分解後の固有モード関数成分及び残量における極高欠損率のデータに対応する箇所のデータを補完し、改めて重ね合わせてマージして時系列的に連続した完全時間特徴モードを取得する。
【００５５】
５）上記Ｓ４ステップに基づき、Ｓ２において取得された空間特徴モード及びＳ３において取得された完全時間特徴モードを利用して、再構築によって時間と空間が完全にカバーされた静止海色衛星データセットを取得する。
【００５６】
本発明の再構築方法をＤＩＮＥＯＦ−Ｇと命名し、従来の代表的なＤＩＮＥＯＦ方法（
2021533427.tif 000077
、ＢＡＲＴＨ  Ａ、ＲＩＸＥＮ  Ｍ、ｅｔ  ａｌ．Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  ｏｆ  ｉｎｃｏｍｐｌｅｔｅ  ｏｃｅａｎｏｇｒａｐｈｉｃ  ｄａｔａ  ｓｅｔｓ  ｕｓｉｎｇ  ｅｍｐｉｒｉｃａｌ  ｏｒｔｈｏｇｏｎａｌ  ｆｕｎｃｔｉｏｎｓ：  ａｐｐｌｉｃａｔｉｏｎ  ｔｏ  ｔｈｅ  Ａｄｒｉａｔｉｃ  Ｓｅａ  ｓｕｒｆａｃｅ  ｔｅｍｐｅｒａｔｕｒｅ［Ｊ］．Ｏｃｅａｎ  Ｍｏｄｅｌｌｉｎｇ．２００５、９（４）：３２５−３４６．）に比べて、保留モード数（Ｍｏｄｅ）、再構築精度（ＭＳＥ、ＲＭＳＥ、ＭＡＥ）、信号解釈率及び再構築率の面でいずれも改良され、下記表１に示される。再構築精度の面で８％増加し（ＭＳＥで算出する）、データ再構築率が５８％から９４％に増加し、３６％増加した。
【００５７】
表１
2021533427.tif 000078
【００５８】
本発明が提供するＤＩＮＥＯＦ−Ｇ方法における各ステップの効果及び再構築結果への影響を更に分析するために、ステップ及び実験パラメータを調整する方式を用い、異なるやり方の比較実験を設計し、具体的な計画及び試験結果を表２に示し、表中の試験パラメータが実行するステップを示す（欠損率パラメータを含む）。
【００５９】
表２
2021533427.tif 000079
【００６０】
それらのうち、試験２のパラメータが代表的なＤＩＮＥＯＦに一致し、試験５のパラメータが本発明が提供するＤＩＮＥＯＦ−Ｇに一致する。試験パラメータの順序が試験中の操作プロセスの順序に一致する。具体的な試験パラメータ及び試験結果についての分析によれば、試験２と試験３とを比較すると、異常ピクセルの検出によって再構築精度を向上させることができ、一部のデータを失ったが、再構築率への影響がより小さく、つまり異常ピクセルの検出による原データの時空間欠損率への影響が限られていることが分かる。試験３と試験４とを比較すると、ラプラス平滑化フィルタリングによって更に再構築精度を向上させるとともに、保留モード数を大幅に増加させることに対応して、信号解釈率も向上させる。ラプラス平滑化フィルタリングによって夜間隙間及びデータ連続欠損によるデータギャップをよく平滑化することができ、原データの一部の特徴を回復することが分かる。試験４と試験５とを比較すると、試験５は試験４の再構築結果に基づき、時間モード二次分解補間を利用して空間的な欠損率が９５％時点より大きいデータを再構築することが分かる。試験５は試験４と比べて、再構築精度が少々低下し、その理由は主に追加のサンプル観測値（欠損率が９５％より高く１００％より低いデータ）を導入したためであるが、依然として同様のサンプルを使用する試験１より優れ、これは時間モード二次分解補間結果の精度がより高いことも反映する。
【００６１】
再構築効果を直感的に示すために、本発明のＤＩＮＥＯＦ−Ｇ再構築方法によって杭州湾で２０１７年の３６５日、毎日８時点のＧＯＣＩ  ＴＳＭデータ再構築結果を取得した。図２は２０１７年６月２日のＴＳＭの原データ及び再構築データ（上位８幅が元の画像であり、下位８幅が再構築画像である）であり、この日の８つの原データの総欠損率が７６．６８％であり、８：３０及び１０：３０のデータ欠損率が９５％を超える。以上によれば、時間モード二次分解補間によって非連続欠損データを再構築するとき、９：３０のデータに対応する時間特徴モードを使用して補間を行うことがないが、データが依然として良好な連続性を有し、動的特徴が明らかであり、これは、平均値を使用して高欠損率のデータを再構築する代表的なＤＩＮＥＯＦ方法では実現できない効果である。１５：３０のデータにおける雲跡尾流（中国語：云迹尾流）に対して、再構築データもよくフィルタリング処理される。原データの有効ピクセル値と再構築データとが同じである領域の再構築ピクセル値に対して精度評価を行い、評価結果を図３に示し、元の値及び再構築予測値がｙ＝ｘ直線付近に集中しており、ＲＳＭが０．００９８であり、ＲＭＳＥが０．０９８９であり、ＭＡＥが０．０６４５であり、これは、本発明の方法の各精度指標がより理想的なものであり、静止海色衛星データの再構築にとって非常に重要な実際の応用価値を有することを示している。
【００６２】
以上の実施例は本発明の好適な技術案であって、本発明を制限するためのものではない。当業者は、本発明の趣旨及び範囲から逸脱することなく、更に種々の変更や変形を行うことができる。従って、同等置換または等価変換によって得られた技術案は、いずれも本発明の保護範囲内に含まれる。
【図１】
2021533427.tif 000106
【図２】
2021533427.tif 000107
【図３】
2021533427.tif 000108
【国際調査報告】
2021533427.tif 000109
2021533427.tif 000110
2021533427.tif 000111
2021533427.tif 000112