JP2008058109A - 観測データ推定方法及び観測データ推定プログラム - Google Patents

Abstract

【課題】気象観測や大気質測定を観測点網の観測地点以外の任意の地点において精度良く推定する。
【解決手段】観測点網における観測地点以外の任意の対象地点において観測データの推定を行うに際し、対象地点から一定の範囲内にある観測地点で観測された時系列の観測データに対し、主成分分析をおこなって独立変動成分の分離を行い、変動成分毎に空間分布（主成分分析における固有ベクトル）の補間をし、補間誤差が最小となる補間値を求め、変動成分毎の補間値の合成をすることにより、補間誤差を最小にする。
【選択図】図２

Description

本発明は、観測データ推定方法及び観測データ推定プログラムに関する。さらに詳述すると、観測点網による観測地点以外の任意の地点における気象要素や大気環境濃度等の観測データの推定を可能とする観測データ推定方法及び観測データ推定プログラムに関する。

現在、気象要素の観測は、全国約１５０カ所の気象官署と全国約１３００カ所のＡＭｅＤＡＳ（Automated Meteorological Data Acquisition System：アメダス）で行われている。また、大気環境の測定は、全国約２０００カ所の国や自治体の大気質測定局により観測が行われており、一定時間毎の時系列の観測データをインターネット等を通じて入手することが可能となっている。以下、ＡＭｅＤＡＳ、大気質測定局等の気象観測、大気観測を行う観測システムを観測点網といい、観測点網により観測及び測定される各種のデータ（降水量、風向、風速、大気中の二酸化硫黄SO₂・二酸化窒素NO₂濃度等の大気環境濃度等）を総称して観測データという。

観測点網における観測地点数は、例えばＡＭｅＤＡＳの場合、降水量を観測するのは全国に約１３００カ所、風向、風速、気温等を観測するのは全国に約８００カ所に限られる。したがって、観測地点以外の任意の地点（以下、対象地点という）で観測データを必要とする場合は、最も近い観測地点でのデータをそのまま用いたり、対象地点で必要に応じて観測データの実測を行っているのが現状である。また、対象地点において観測データの推定を行う場合もあるが、当該推定は、近接する複数の観測地点での観測データを単に線形補間を行うだけのものである。

また、任意の対象地点での日射量を推定する日射量推定システムが提案されている（特許文献１）。特許文献１に記載の技術は、日射量が実測される地点において、最初にその実測値と、その地点に関するＧＭＳ画像データを使用して所定の計算式で求めたその地点の日射量の推定値とを比較し、当該推定値と当該実測値との誤差が所定の許容量以下になるまで、計算式内のパラメータの値を自動的に調整し、実測地点に近接する複数の他の任意の地点のそれぞれに関して、実測地点において調整されたパラメータの値を使用して、計算式により日射量の推定値を計算し、実測地点の誤差を減少させるパラメータの値を、実測地点に近接する他の任意の地点の日射量の計算に使用することによって、他の任意の地点の日射量を推定するものである。
特開平１１−２１１５６０

しかしながら、対象地点における観測データ（実測した場合）と近接の観測地点での観測データには、誤差が大きいため、そのまま用いることはできない。

また、対象地点において実測を行うには膨大な人件費等のコストがかかる。更に、一定期間（数ヶ月間や一年間）の観測データの累計を要することが通常であるので、当該対象地点での観測データを得るのに膨大な時間及びコストを要することとなる。例えば、発電所のリプレース時等においては、１年間の気象環境の観測が行われている。このような場合に、当該地点で実測を行っていたのでは、気象環境の観測のために期間を１年間余計に要することとなっていた。また、実測には多額のコストを要するという問題があった。

また、対象地点で観測データの推定を行う場合であっても、従来の推定手法は近接する複数の観測点での観測データを単に線形補間を行う（平均値をとる）という単純なものに過ぎず、精度良く推定を行うことはできなかった。

また、特許文献１に記載の技術は、その推定対象は日射量に限られるものであり、他の観測データの推定を行うことができなかった。

そこで、本発明は、観測点網における観測地点以外の任意の対象地点において、観測データの高精度の推定を可能とする観測データ推定方法及び観測データ推定プログラムを提供することを目的とする。

かかる目的を達成するため、請求項１記載の観測データの推定方法は、観測点網における観測地点以外の任意の対象地点において、対象地点から一定の範囲内にある観測地点での時系列の観測データについて主成分分析を行うことにより固有値番号毎の空間分布を求める主成分分析処理と、対象地点において固有値番号毎に、空間分布を、補間誤差が最小となるように補間値を求める最適空間補間処理を行い、更に、求めた固有値番号毎の補間値を合成することにより、対象地点において、補間誤差が最小となる観測データの推定を行うようにしている。

また、請求項２に記載の観測データの推定プログラムは、任意に選択した対象地点の位置情報、観測点網における該対象地点の近傍の観測地点の位置情報、該観測地点での時系列の観測データを記憶装置に記憶させ、対象地点から一定の範囲内にある観測地点での時系列の観測データを主成分分析することにより固有値番号毎の空間分布を抽出し記憶する主成分分析処理と、対象地点において固有値番号毎に、空間分布を、補間誤差が最小となる補間値を求め記憶する最適空間補間処理と、固有値番号毎の補間値を合成する処理とをコンピュータに実行させることにより対象地点での観測データの推定を行うものである。

したがって、観測点網における観測地点以外の任意の対象地点において観測データの推定を行うに際し、対象地点から一定の範囲内にある観測地点で観測された時系列の観測データに対し、主成分分析をおこなって独立変動成分の分離を行い、変動成分毎に空間分布（主成分分析における固有ベクトル）の補間をし、補間誤差が最小となる補間値を求め、変動成分毎の補間値の合成をすることにより、補間誤差が最小となる観測データの推定を行っている。

以上説明したように、本発明にかかる観測データ推定方法及び観測データ推定プログラムによれば、任意の対象地点において、大気環境濃度、風速、風向などの各種の観測データを実測することなく、補間誤差を最小とする高精度の推定を可能とすることができる。

また、実測を省くことで、大幅なコストの削減を図ることが可能となる。また、実測にかかる期間を削減することが可能となる。例えば、発電所の建設予定地であれば、迅速な着工を図ることが可能となる。

以下、本発明の構成を図面に示す実施の形態に基づいて詳細に説明する。

本発明の観測データ推定方法は、観測点網における観測地点以外の任意の対象地点において、対象地点から一定の範囲内にある観測地点での時系列の観測データについて主成分分析を行うことにより固有値番号毎の空間分布を求める主成分分析処理と、対象地点において固有値番号毎に、空間分布を、補間誤差が最小となるように補間値を求める最適空間補間処理を行い、更に、求めた固有値番号毎の補間値を合成することにより、対象地点において、補間誤差が最小となる観測データの推定を行うようにしている。尚、「一定の範囲内にある観測地点」とは、任意に設定することが可能であり、例えば、１枚の地図データ上に表せる範囲にある観測地点、対象地点からの一定の距離内に含まれる観測地点、対象地点の近傍のＸ箇所（Ｘは任意に設定可能）の観測地点等とすることができる。

従来の観測データの推定では、近接する観測地点での観測データを単純に線形補間することにより推定するだけで、精度良く観測点間での推定を行うことはできなかった。したがって、実際にはある一定以上の推定精度を要求される場合には、当該地点で独自に観測データを実測する必要があり迂遠であった。

そこで、本願発明者は、主成分分析を用いて各観測地点での時系列の観測データから主成分を導出し（以下、主成分分析処理という）、さらに導出した主成分毎に空間補間を行う（以下、最適空間補間処理という）ことで、従来、不可能であった任意の対象地点での観測データを、従来に比べ精度よく推定することが可能となることを知見した。

現在、ＡＭｅＤＡＳでは、降水量、気温、日照時間、風向、風速、積雪等の気象要素の観測を行っている。また、自治体大気質測定局では大気中の二酸化硫黄SO₂、二酸化窒素NO₂濃度等の大気環境濃度の測定を行っている。本発明の観測データ推定方法によれば、上記すべての推定が可能となるが、次の３つの条件を満たすデータの推定に特に秀でているという特徴を有する。
（１）空間的に離れた観測地点間でも観測データに有る程度の相関があること。
（２）距離が遠くなるにつれ、その相関が弱くなること。
（３）複数の要因により変動するものであること。
具体的には、大気環境濃度、風速、風向のデータの推定に最も適している。尚、上記観測データは、例示でありこれに限られるものではない。

本発明の観測データ推定方法の概要は、次のようなステップとなる。
（ステップ１）各観測地点の時系列の変動分を相関のない独立な成分の重ね合わせで表す。
（ステップ２）重ね合わせの係数よりその空間的な分布スケールを求める。
（ステップ３）成分毎に分布スケールに対応した補間を行う。
（ステップ４）対象地点の値を成分の合計として求める。
このように、主成分（変動要因）毎に分離して補間を行うことにより、補間誤差を最小とし、精度の高い推定を行うことが可能となるものである。

以下に主成分分析および最適空間補間の方法について詳細に述べる。先ず、主成分分析について説明する。主成分分析とは、相関関係にあるいくつかの要因を合成して、いくつかの成分にし（独立変動成分分離）、その総合力や特性を求める方法であり、本実施形態においては、各観測点での時系列データを成分に分けること、換言すれば、各観測点のデータを時系列に重ね合わせることで、互いに相関のない主成分（Z_i(t)）を合成するものである。

本実施形態における観測データである気象データ、大気質データ等は、主成分分析の例題として用いられる「身長」と「体重」との関係のように、概念的に複数の主成分（この場合、「体の大きさ」、「肥満度」等）を導出することができるものではない。しかしながら、各観測点の測定値は様々な要因の重ね合わせであると考えられる。例えば、大気質のデータであれば、複数の汚染源の影響を受けており、しかもその影響は測定局（観測地点のこと）及び測定時刻により異なるといえる。また、主成分分析処理によれば、その要因（即ち、主成分）毎にその影響の分布状況を把握できる。

図１に主成分分析を概念的に説明するための図を示す。例えば、図１（ａ）に示すように、対象地点１０、観測地点（X₁〜X₃）での観測の時系列データがある場合に３つの主成分、例えば観測初期に比較的緩慢な変動のある第１主成分Z₁、中期に比較的速い変動のある第２主成分Z₂、後期に比較的緩慢な変動のある第３主成分Z₃とが存在するとする（図１（ｂ）、（ｃ））。

測定時系列は、数式１に示すように主成分の重ね合わせで表される。ここで、重ね合わせで用いられる係数からなるベクトル（数式２）を主成分ベクトルという。

数式２で示される主成分ベクトルは、各観測点での時系列間の共分散 <x_ix_j>_ｔ を要素とする分散共分散行列の固有ベクトルとなる。固有ベクトルであるので、主成分ベクトルは互いに直交し、また通常は単位ベクトルとする（数式３）。尚、Z_i(t)²の時間平均λ_iが固定ベクトルu_iに対応する固有値であり、全地点での振幅の時間平均に対応する。

上記主成分ベクトルは、対応する独立変動成分の空間分布パターンを示すものである。したがって、これを空間的に補間することにより、その変動成分の対象地点での値を推定することができる。

次いで、最適空間補間について説明する。図２（ａ）〜（ｃ）に示すように、主成分によってその空間スケールが異なる。ここで、実線は、等値線（主成分ベクトル、空間分布）を、波線１１は、線形補間可能域（波線の半径が重み半径）を示す。例えば、図２（ｂ）に示す第２主成分のように対象地点１０近傍で急峻な変化をする空間分布をしている場合、即ち、対象地点１０の近傍に変動の中心がある場合は、変動の中心１２を含めるように広い範囲の観測値を用いて補間を行うと急峻な変化の影響を受け、補間誤差が大きくなる。また、図２（ｃ）に示す第３主成分のように対象地点１０付近でなだらかな変化をするような広い範囲で空間的次関数で近似可能な場合に、狭い範囲で補間近似すると、各データの誤差の影響を受けやすく、最適な補間とはならない。即ち、補間可能域（波線１１）は、変動の中心１２を含まないように半径を設定することが必要となる。

最適空間補間は、対象地点近傍での変動の空間的スケールに応じて、最適な半径を持った重み関数を用いてデータに重み付けをすることで補間精度を向上させるものである。

具体的には以下のような処理を行う。重み半径に対して重み関数 w(r;q) を数式４により定義する。尚、qは位置ベクトルである。

更に、各測定局i毎の値{y_i}_i=1,Nに対して、原点（対象地点）での w(r) の重み付空間平均 <y>_w(r) を数式５により定義する。

ここで、<y>_w(r) の値はrが最遠測定局までの距離に比べて十分に大きい場合は全測定点の平均値<y> 、rが最寄測定局までの距離に比べて十分に小さい場合は、最寄測定点での値となる。

また、rを次第に増大させたときの<y>_w(r)の変化を指数関数 a exp(-r/ρ)+<y> で近似した場合のρの最適値ρ^＊が、yの空間分布のスケールとなる。尚、後述するように最適値ρ^＊の設定方法は、当該手法に限られるものではない。

得られた最適値ρ^*に基づく重み関数（数式６）を用いた場合、空間分布y=(y1,y2...)^tを位置ベクトルqの１次関数で補間する数式７により表される。

更に、係数ベクトルaは、位置の共分散行列などを用いて数式８により表される。

求めるのは対象地点（位置ベクトルo）での推定値であるので、数式９により表すことができる。

また、この補間による推定誤差ε^２は、数式１０で与えられる。尚、後述のように重み半径を検索し、推定誤差をε^２最小にするような最適値ρ^＊を求めるようにしても良い。

本実施形態では、主成分ベクトルu_jをyに当てはめることで、空間分布スケールρ_j ^*を求め、対象地点での補間値u^_j(o)および補間誤差ε_j ^２を得ることができる。全ての変動成分を合成した値が、最終的に対象地点での観測値の推定値となるから、その値x^（o）及び補間誤差σ^２は数式１１により表される。

次いで、本発明の観測データ推定プログラムが行う処理を図４〜図６のフローチャートを用いて説明する。

本発明の観測データ推定プログラムは、例えば、次のような観測データ推定装置１により実行される。図３に観測データ推定装置１の構成の一例を示す。観測データ推定装置１は、ディスプレイ等の出力装置２と、キーボード、マウス等の入力装置３と、演算処理を行う中央処理演算装置（ＣＰＵ）４と、計算中のデータ、パラメータ等が記憶される主記憶装置（ＲＡＭ）５と、計算結果等が記録される補助記憶装置としてのハードディスク６、外部との通信を行う通信インタフェース７等を備えている。尚、主記憶装置５及び補助記憶装置６を総称して、単に記憶装置ともいう。上記のハードウェア資源は例えばバス８を通じて電気的に接続されている。

また、本発明の観測データ推定プログラムは、補助記憶装置６に記録されており、当該プログラムがＣＰＵ４に読み込まれ実行されることによって、コンピュータが観測データ推定装置１として機能する。その実行の際に必要なデータは、ＲＡＭ５にロードされる。また、図示はしないが、観測データ推定装置１は、主成分分析処理を行う主成分分析手段、最適空間補間処理を行う最適空間補間手段を備えている。尚、上述のハードウェア構成は一例であってこれに限られるものではない。

プログラムの実行にあたっては、ＡＭｅＤＡＳなどの観測点網から観測地点の位置情報（緯度、経度、高度等）及び各観測地点での時系列の観測データ等を取得し、初期データとして入力する。入手可能なデータは、提供元により異なるが、少なくとも位置情報（緯度、経度）及び時系列の観測データが入手可能であればよい。例えば、測定局の位置は、一部は住所のみ入手可能なものも存在するが、そのような場合でも住所から緯度、経度といった位置情報を求めておくことで実行することが可能となる。本実施形態では取得する時系列データは例えば１時間毎とするが、これには限られず、例えば、ＡＭｅＤＡＳによる風向、風速の発表間隔である１０分間隔で取得するようにしてもよい。

本発明の観測データ推定プログラムでは、図４に示すように、先ず初期データの入力を行う（Ｓ１）。本実施形態においては、以下に挙げる初期データを補助記憶装置６等に記憶させる。初期データとしては、x₀,y₀，{x_i,y_i,p(i,t）｝_i,tを入力する。ここで、iは地点番号、tは時刻を示す。また、x₀,y₀は、補間対象位置の座標を、x_i,y_iはｉ番目の観測地点を、p(i,t）は測定値を示す。即ち、x₀,y₀は、対象位置の画像データ上でのx,y座標（xは緯度、yは経度）を、{x_i,y_i,p(i,t）｝_i,tは、各測定局における時系列の測定値を意味する。本実施形態では、測定値pは、大気質データ（大気中のSO₂濃度）とするが、他のデータの推測を行う場合は、pを当該データに変更すればよい。

次に、主成分分析処理（Ｓ２）を行う。主成分分析処理の詳細フローを図５に示す。

主成分分析処理（Ｓ２）では、先ず、測定値p(i,t)を入力引数とし（Ｓ２０１）、p(i,t)の時間平均値である<p>(i)と、p(i,t)とp(j,t)の共分散であるP(i,j)を数式１２により求める時間的統計処理を行う（Ｓ２０２〜Ｓ２０３）。
＜数１２＞
<p>(i)＝Σ_t p(i,t)/T
P(i,j)＝Σ_t (p(i,t)-<p>(i))×(p(j,t)-<p>(j))/T
尚、i,jは地点番号(=1...N)、tは、時刻(=1...T)を示す。

得られた平均<p>(i)と共分散P(i,j)に基づいて、Ｐの固有値解析を行い、固有値λ_k、固有ベクトルu_k(i)を求める（Ｓ２０４〜Ｓ２０５）。本実施形態における固有値解析は、対象行列を入力として固有値及びそれに対応する固有ベクトルを出力とするものであり、解析手法は既知（ハウスホルダー法、ヤコビ法等）または新規の手法を用いることとすれば良く、特に限られるものではない。尚、Ｐは分散共分散行列(i,j 成分がP(i,j)である行列)を示し、固有ベクトルu_k(i)は、空間部分を示す。尚、kは固有値番号(=1...N;但し、固有値の大きい順とする)を示す。

得られた固有値λ_kと固有ベクトルu_k(i)に基づいて数式１３により主成分の計算を行い、k番目の主成分時系列Z_k(t)を求める（Ｓ２０６〜Ｓ２０７）。
＜数１３＞
Z_k(t)＝Σ_i p(i,t) u_k(i)
{λ_k,u_k(i),Z_k(t)}_i,k,tを出力引数として（Ｓ２０８）、主成分分析処理（Ｓ２）は終了する。

次に、求めた主成分毎に最適空間補間処理（Ｓ３〜Ｓ６）を行う。先ず、固有値番号をk=1として初期化を行なう（Ｓ３）。

次に、最適空間補間処理（Ｓ４）を行う。最適空間補間処理の詳細フローを図６に示す。最適空間補間処理（Ｓ４）は、固有値番号kをインクリメントしながらループ処理を行うものである。

最適空間補間処理では、先ず、対象地点の座標x_o,y₀、及び{x_i,y_i,u(i)}_iを入力引数とする（Ｓ４０１）。尚、x_i,y_i はi番目の観測地点を、u(i)はi番目の観測地点での値(補間値)を示す。

次に、重み半径r及び補間誤差E_Mの初期化を行う（Ｓ４０２）。ここで、rは最適値を検索するため走査する重み半径であり、r_iは重み半径初期値を示す。また、E_Mは最適重み半径の候補値r_Mに対応する、最小二乗補間の誤差である（Ｓ４０６参照）。重み半径初期値及び補間誤差初期値は数式１４で定義される。尚、βは、十分な大きさの重み半径初期値を得ることを目的として任意に設定可能なパラメータであり、３程度が好ましい。また、重み半径rに依存するi番目の観測地点での重みw_iを数式１５により定義し、全ての観測地点のうち最大の重み値w_max及び全ての観測地点での重みの合計値Wを数式１６により求める（Ｓ４０３〜Ｓ４０４）。また、R_iは補間対象位置までの距離であり、数式１７により定義される。尚、上記「十分な大きさの重み半径を設定する」とは、各観測地点の重みがほぼ同じになることを意味しており、数式１５で定義される重みw_iは、最大値が１であるので、値の一番小さくなる最遠方の観測地点でも値が１に近くなるようにr_iを設定すれば良いこととなる。即ち、r_xを補間対象地点から最遠方の観測地点までの距離、w_xをその場合の重みとすると、数式１８の近似が成立する。よって、(R_x/r_i)²<0.1の程度であれば、w_xは１に近い値となり、r_iはr_xの３倍程度でよいこととなる。また、本実施形態では、重みw_iを数式１５により定義しているが、重みの設定方法の一例に過ぎず、数式１９が、収束するような単調減少、正値の関数であればよいこととなる。
＜数１４＞
r_i = β×最遠方の観測地点までの距離
E_M=Σ_i u_i ²（または、E_M=処理システムにおける実数の最大値）
＜数１５＞
w_i = exp(-(R_i/r)²)
＜数１６＞
w_max = Max_i w_i
W = Σ_i w_i
＜数１７＞
Ri = ((x₀-x_i)²+(y₀-y_i)²)^1/2
＜数１８＞
w_x ~ 1-(R_x/r_i)²

次に、数式２０により収束判定を行う（Ｓ４０５）。尚、数式２０により収束判定を行うのは以下の理由による。最小二乗法による近似式を作る際に、パラメータを決定するための測定値の自由度は、パラメータ数よりも多いことが必要となる。本実施形態では、補間の近似式として２次元の位置（緯度、経度）の線形関数としているため（数式２２参照）、パラメータはu,a,bの３つである。
＜数２０＞
W > 3w_max

更に、最寄地点の重みw_maxに対するw_iの重みの地点の補間への貢献はw_i /w_maxであり、重みの小さい地点は、データの自由度への貢献も、小さいこととなる。遠くに重みが限りなく０である地点があったとしても、近似パラメータの安定性には何の貢献もないからである。また、重みは相対値でしかないため、一番大きな重みを持つ最寄地点の測定データが自由度１相当の貢献があるとして、w_i/w_maxが補間の自由度への貢献となると考えられる。よって、全地点の重みつきデータとしての自由度は数式２１で示されることになり、重みつきデータとしての自由度が、補間パラメータ数３を越える必要があるため、w/w_max > 3が、安定補間の条件となり、これを下回る場合は補間が不安定となり、不適な重み半径となるからである（即ち、数式２０）。
＜数２１＞
Σ_i (w_i/w_max) = (Σ_i w_i)/w_max = w/w_max

Wが3w_maxを超える場合（Ｓ４０５；Ｙｅｓ）は、数式２１によりiに関する重み付き平均値xの計算を行う（Ｓ４０６）。同様に平均値y,uについても求める。
＜数２１＞
x = Σ_i w_i x_i

更に、最小二乗法による係数及び誤差の計算を行う（Ｓ４０６）。尚、a,bは補間係数u'iを数式２２により定義したときに補間誤差Eを最小にする係数である。補間誤差Eは数式２３により表される。係数a,bは数式２４で定義される。
＜数２２＞
u'_i＝u＋a(x_i−x)＋b(y_i−y)
＜数２３＞
E＝Σ_i w_i (u'_i−u)²
＜数２４＞
a＝(S_YY S_UX − S_XY S_UY)/D, b＝(S_XX S_UY − S_XY S_UX)/D,
但し、S_XX＝Σ_i w_i(x_i−x) ², S_YY＝Σ_i w_i(y_i−y) ², S_XY＝Σ_i w_i(x_i−x)(y_i−y),
S_UX＝Σ_i w_i(x_i−x)(u_i−u), S_UY＝Σ_i w_i(y_i−y)(u_i−u), D＝S_XX S_YY−S_XY ²

次に、E_M>Eの場合（Ｓ４０７；Ｙｅｓ）は、最適重み半径の候補値r_Mを当該ループ内での計算結果に更新する。併せて、最適重み半径の候補値r_Mに対応する、最小二乗補間の係数a_M,b_M,誤差E_Mの更新を行う（Ｓ４０８）。一方、E_MがE以下の場合（Ｓ４０７；Ｎｏ）は、最適重み半径の候補値r_Mの更新は行わない。

次いで、数式２５により重み半径rの更新を行う（Ｓ４０９）。尚、αは重み半径走査のために重み半径を減少させる係数パラメータであり、０＞α＞１とする。尚、α＝０．７〜０．８程度が好ましい。
＜数２５＞
r = α×r

更に、終了判定を行う。本ステップでは、数式２６により終了判定を行うようにしている。rがr_min / 2以上の場合（Ｓ４１０；Ｎｏ）は、Ｓ４０３へ戻り処理が継続する。ここでr_minは、最寄り地点までの距離を表す。これは、最寄り地点までの半分以下となるような重み半径を設定し、補間を行っても意味がないからである。
＜数２６＞
r > r_min / 2

一方、上記収束条件を満たした場合（Ｓ４０５；Ｎｏ）及び、上記終了判定を満たした場合（Ｓ４１０；Ｙｅｓ）は、数式２７により補間推定を行う（Ｓ４１１）。ここでu₀は、補間対象位置（推定位置）でのuの推定値を示している。
＜数２７＞
u₀＝u＋a(x₀−x)＋b(y₀−y)
最終的に、u₀,r_M,a_M,b_M,E_Mを出力引数として最適空間補間処理は終了する（Ｓ４１２）。

予め設定された使用する主成分数Kを超えない場合（Ｓ５；Ｎｏ）は、ｋに１を加えて（Ｓ６）、最適空間補間処理（Ｓ４）を行う。一方、使用する主成分数Kを超える場合は、数式２８により時系列合成（Ｓ７）を行い、最終的な対象地点での推定値p₀(t)を得るものである。
＜数２８＞
p₀(t)＝Σ_k=1,K Z(k,t) u_k(0)
以上で、本発明の観測データ推定プログラムが実行する処理を終了する。

尚、上述の実施形態は本発明の好適な実施の例ではあるがこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。また、上述の演算式は一例であり、本発明の要旨を逸脱しない範囲において種々変形実施可能である。例えば、最小２乗法による補間（Ｓ４０６）において、数式２９，３０に示すような非線形関数や、数式３１，３２に示すような緯度、経度以外のパラメータを含めたものを用いるようにしても良い。
＜数２９＞
ax + by + c + px² + qy² + rxy
＜数３０＞
ax + by + c + p sin(μx + θ)+ q sin(νx + φ)
＜数３１＞
ax + by + c + pz
但し、zは標高示す。
＜数３２＞
ax + by + c + pn
但し、nは人口密度を示す。

以上、本発明の観測データ推定方法、推定プログラムによる大気質データの推定について説明したが、風データ（風向、風速）の推定方法について以下に説明する。

ＡＭｅＤＡＳの風の観測点は一般にそれぞれ地上およそ１０mの高度（地上風）において１時間毎に風向・風速の値(１０分平均値)を取得している。値の分解能は、風向が１６方位、風速が１m/sである。観測点の設置密度は平均２５km程度の間隔となっている。

風データの推定では、先ず、AMeDAS の風向(北を０とする方位D[度])、 風速(S[m/s])より各地点の風速の東西成分U_iおよび南北成分V_iを数式３３、３４により求める。尚、、i は測定地点番号(1〜N)を示す。
＜数式３３＞
Ui = −S_i sin(D_i)
＜数式３４＞
Vi = −S_i cos(D_i)

次に、上述の観測データ推定方法における測定値p_iの代わりにU_iとし、これについて緯度、経度について空間補間し、対象地点での値Uを推定する。同様に、測定値p_iの代わりにV_iとして空間補間し、対象地点での値Vを推定する。

最後に、風速Sび風向Dは、東西、南北両成分より合成して求める（数式３５）。
＜数式３５＞
S = (U²+V²)^1/2
D=cos^-1(-V/S)= sin^-1(-U/S)

AMeDASは連続測定を行っており、１月程度の遅れでデータが一般に提供されているため、本発明の観測データ推定方法によれば１年分のデータをすぐに得ることが可能となる。

また、本発明の観測データ推定方法により推定された風速、風向（地上風）に、公知または新規の地上気象観測からの上層気象観測の推定手法を組み合わせることにより、地上風だけでなく上層風の推定を行うことが可能となる。

更に、本発明の観測データ推定方法により大気安定度の推定が可能となる。一般に大気安定度は、全天日射量、放射収支量及び風速より求めることが可能である。よって、本発明に公知又は新規の全天日射量及び放射収支量の推定手法を組み合わせることにより、大気安定度の推定が可能となる。

例えば、大気環境アセスメントにおいては、発電所のリプレース時に環境負荷が低減される場合や、新鋭の類似施設と同等の実行可能なよりよい排煙処理等の技術を取り入れ、十分な環境保全措置を計画している発電所の建設コスト削減、期間短縮が要求されている。このような場合に、本発明の観測データ推定方法によれば、発電所排煙の拡散シミュレーションを実施する際に上層拡散場の設定に用いる風向・風速および大気安定度を、発電所計画地点で気象観測を実施せずに推定することが可能となり、発電所の建設コスト削減、期間短縮が可能となる。

本発明の観測データ推定方法による観測データの推定値と、実測値との比較実験を行った。A、Bの両火力発電所を対象に半径５０km以内の大気質測定局のSO₂データについて、本発明の観測データ推定方法により推定を行った。

A発電所周辺は、工業地帯でSO₂の排出源が多く存在するため、測定局も密に分布している。本実験では、A発電所に最寄りの測定局を対象地点と想定し(以下、ケースA、データセットAとする)、それ近傍の測定局での観測データを補間して推定した値の比較を行った。併せて、必要な測定局の密度を調べるため、推定には、以下に上げるように複数の測定局のセットパターンを用いて行った。
Bセット:５０km以内の全測定局（２km程度の間隔）のセット
Cセット:７km程度の間隔となる測定局のセット
Dセット:１５km程度の間隔となる測定局のセット
各セットにつき、本手法による観測データの推定と最寄りの測定局の値そのもので代替する方法とを、評価地点(ケースA)と１年間比較した結果を表１に示す。

すべてのケースを通じて平均値は５．６〜７．６の範囲に収まっており年平均値としてはいずれのセットを用いても環境濃度としては十分な推定精度があると考えられる。推定誤差、相関係数を見ると、どのセットを用いても単なる最寄り地点の値に比べ本発明の推定方法を用いる方が良好である。また、誤差の違いはBセットとCセットとでは０．１ppb程度とそれほど大きくないがDセットになると やや大きい。相関係数の違いでは、Cセットでは０．８を越えているもののDセットでは０．７を下回ってしまうことから測定局の密度としてはCセット、つまり間隔は７km 程度以下であることが必要と考えられる。

以上より、本発明の観測データ推定方法により、大気質測定局のSO₂データについて精度良く推定を行うことが可能であることを確認した。また、約２km間隔の全地点（測定局）を用いた補間と、７km間隔に間引いた地点での補間は、ほぼ同様の精度で、１５kmまで間引いたデータではやや精度が落ちることを確認した。

また、B発電所周辺は、SO₂の排出源が少ないため、測定局の密度も平均的な間隔が２０km程度と低い。大気濃度の場合５０kmを超えて遠方の地点が主成分として大きな影響を与える可能性は小さい。そのため、５０km以遠の測定局を用いる場合は局数を２０程度とした。B発電所に最寄りの測定局を擬似評価地点としてそれ以外の測定局の値を補間して推定した値および擬似評価地点の最寄り測定局の値の比較を行った(表２)。

平均値は評価地点の５．１６ppbに対して両代替値とも違いは小さい。また、本手法での推定値と評価地点の値の相関係数は０．４と小さい値にとどまるものの、推定誤差も、２．３ppbとA発電所に比べても小さい。これらの差は、いずれも環境基準に比べ十分に小さく、十分な推定精度があると考えられる。

以上より、本発明の観測データ推定方法により空間的に補間した値を評価地点での環境濃度観測に対する代替として用いることが可能であることを確認できた。

主成分分析処理を説明するための図及びグラフであり、（ａ）は、対象地点と観測地点との位置関係を示す図、（ｂ）は、各観測地点での時系列毎の観測地の変動を示すグラフ、（ｃ）は、主成分分析により抽出した各主成分を示すグラフである。 最適空間補間処理を説明するための図であり、（ａ）は、対象地点の左下方に変動のピークがある図、（ｂ）は、対象地点の左直上に変動のピークがある図、（ｃ）は、対象地点の右下遠方に変動のピークがある図を示している。 観測データ推定装置のハードウェア構成図の一例である。 本発明の観測データ推定プログラムが行う処理の一例を示すフローチャートである。 主成分分析処理の処理の一例を示すフローチャートである。 最適空間補間処理の処理の一例を示すフローチャートである。

符号の説明

１０ 対象地点
１１ 空間分布
１２ 変動の中心

Claims (2)

1. 観測点網における観測地点以外の任意の対象地点において、前記対象地点から一定の範囲内にある前記観測地点での時系列の観測データについて主成分分析を行うことにより固有値番号毎の空間分布を求める主成分分析処理と、前記対象地点において前記固有値番号毎に、前記空間分布を、補間誤差が最小となるように補間値を求める最適空間補間処理を行い、更に、求めた前記固有値番号毎の前記補間値を合成することにより、前記対象地点において、補間誤差が最小となる観測データの推定を行うことを特徴とする観測データ推定方法。
2. 任意に選択した対象地点の位置情報、観測点網における該対象地点の近傍の観測地点の位置情報、該観測地点での時系列の観測データを記憶装置に記憶させ、前記対象地点から一定の範囲内にある前記観測地点での前記時系列の観測データを主成分分析することにより固有値番号毎の空間分布を抽出し記憶する主成分分析処理と、前記対象地点において前記固有値番号毎に、前記空間分布を、補間誤差が最小となる補間値を求め記憶する最適空間補間処理と、前記固有値番号毎の前記補間値を合成する処理とをコンピュータに実行させることにより前記対象地点での観測データの推定を行うことを特徴とする観測データ推定プログラム。
   

Images