JP2003296702A - Method for denoising earth observation satellite data, denoising processing program and recording medium recorded with denoising processing program - Google Patents
Method for denoising earth observation satellite data, denoising processing program and recording medium recorded with denoising processing programInfo
- Publication number
- JP2003296702A JP2003296702A JP2002103966A JP2002103966A JP2003296702A JP 2003296702 A JP2003296702 A JP 2003296702A JP 2002103966 A JP2002103966 A JP 2002103966A JP 2002103966 A JP2002103966 A JP 2002103966A JP 2003296702 A JP2003296702 A JP 2003296702A
- Authority
- JP
- Japan
- Prior art keywords
- data
- time
- observation
- observation data
- coefficient
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims description 73
- 238000012545 processing Methods 0.000 claims abstract description 112
- 230000001932 seasonal effect Effects 0.000 claims description 32
- 230000000737 periodic effect Effects 0.000 claims description 19
- 230000008569 process Effects 0.000 claims description 18
- 230000007774 longterm Effects 0.000 claims description 6
- 230000035945 sensitivity Effects 0.000 claims description 4
- 230000006870 function Effects 0.000 description 47
- 238000004364 calculation method Methods 0.000 description 17
- 230000008859 change Effects 0.000 description 13
- 238000001514 detection method Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 13
- 238000004422 calculation algorithm Methods 0.000 description 12
- 230000002265 prevention Effects 0.000 description 9
- 238000003672 processing method Methods 0.000 description 8
- 238000001914 filtration Methods 0.000 description 7
- 238000004088 simulation Methods 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 6
- 238000013480 data collection Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 6
- 239000000284 extract Substances 0.000 description 5
- DFPOZTRSOAQFIK-UHFFFAOYSA-N S,S-dimethyl-beta-propiothetin Chemical compound C[S+](C)CCC([O-])=O DFPOZTRSOAQFIK-UHFFFAOYSA-N 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 4
- 238000003786 synthesis reaction Methods 0.000 description 4
- 238000012731 temporal analysis Methods 0.000 description 4
- 238000000700 time series analysis Methods 0.000 description 4
- 239000002131 composite material Substances 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 241000196324 Embryophyta Species 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000005295 random walk Methods 0.000 description 2
- 238000012502 risk assessment Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- LQEOTHINZCAFAC-UHFFFAOYSA-N (4-methylpiperazin-1-yl)-(3,4,5-trimethoxyphenyl)methanone Chemical compound COC1=C(OC)C(OC)=CC(C(=O)N2CCN(C)CC2)=C1 LQEOTHINZCAFAC-UHFFFAOYSA-N 0.000 description 1
- 239000002028 Biomass Substances 0.000 description 1
- 102100033040 Carbonic anhydrase 12 Human genes 0.000 description 1
- 241000905957 Channa melasoma Species 0.000 description 1
- 238000006424 Flood reaction Methods 0.000 description 1
- 101000867855 Homo sapiens Carbonic anhydrase 12 Proteins 0.000 description 1
- 101000680262 Homo sapiens Transmembrane protein 60 Proteins 0.000 description 1
- 102100022076 Transmembrane protein 60 Human genes 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000012010 growth Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000037211 monthly cycles Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 208000036460 primary closed-angle glaucoma Diseases 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000000779 smoke Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000009628 steelmaking Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Landscapes
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
Description
【0001】[0001]
【発明の属する技術分野】本発明は、衛星データのノイ
ズ除去処理方法、ノイズ除去処理プログラム、ノイズ除
去処理プログラムを記録した記録媒体に係り、特に、高
頻度観測衛星データのノイズ除去処理による植生の季節
変動観測法に関する。BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a satellite data noise removal processing method, a noise removal processing program, and a recording medium on which the noise removal processing program is recorded. Regarding seasonal variation observation method.
【0002】1997年から1998年にかけてインド
ネシアで大規模な森林火災が発生し、煙害等で周辺諸国
も含めて甚大な影響を及ぼした(文献1,2参照。な
お、文献は後述する。以下同様。)。その後も東南アジ
ア地域では森林火災が頻発し、大規模火災の再発生が懸
念されている。現在、国際協力事業団(JICA)の
「インドネシア森林火災予防計画プロジェクト」の第二
フェーズが進行中であり、その中で、「森林火災早期発
見システム」が稼働している。これは、火災発見を衛星
データを用いて行い、森林管理を担当する機関に通報す
るものである(文献2参照)。一方、火災予防の観点か
らは火災危険度の評価が待たれている。A large-scale forest fire broke out in Indonesia from 1997 to 1998, and caused a great impact due to smoke damage including surrounding countries (see References 1 and 2. Note that the references will be described later. The same applies below. .). After that, forest fires frequently occurred in Southeast Asia, and there is concern that a large-scale fire will reoccur. Currently, the second phase of the Japan International Cooperation Agency (JICA) "Indonesia Forest Fire Prevention Plan Project" is in progress, and the "Forest Fire Early Detection System" is in operation. This is to detect fires using satellite data and notify the agency in charge of forest management (see Reference 2). On the other hand, from the viewpoint of fire prevention, evaluation of fire risk is awaited.
【0003】火災危険度評価の基礎となるのは森林植生
図であるが、乾燥に対する森林の感受性を評価する必要
があることを考えると、単なる森林区分では不十分で、
時間を含んだ情報に基づいて危険度を算出することが期
待されている。Forest vegetation maps are the basis of fire risk assessment. However, considering that forest susceptibility to drought needs to be assessed, mere forest classification is not enough,
It is expected to calculate the risk level based on information including time.
【0004】従来の植生の季節変動観測法のための時系
列処理手法の開発に対して、その課題について以下に説
明する。With respect to the development of the conventional time series processing method for the seasonal variation observation method of vegetation, its problems will be described below.
【0005】従来のNOAA/AVHRR あるいはS
POT/vegetation、MODISなどの高頻
度観測衛星からのデータによって広域での植生観測が行
われている(文献3−13参照)。特に、大規模森林火災
や乾燥害などの農業・林業災害は、植生の乾燥の程度に
影響されるため、植生の状態を迅速かつ的確に把握する
手法の開発は防災面で大きな意味を持っている。しか
し、光学センサーでは、雲やヘイズの影響が避けられ
ず、的確な植生分布およびフェノロジーの把握が困難で
ある。特に、熱帯を対象とする場合、数ヶ月間にもわた
って雲の影響があり、雨期にはほとんど有効なデータが
観測できないこともある。Conventional NOAA / AVHRR or S
Vegetation observation is carried out over a wide area using data from high-frequency observation satellites such as POT / vegetation and MODIS (see References 3-13). In particular, agricultural and forestry disasters such as large-scale forest fires and drought damage are affected by the degree of vegetation dryness, so the development of a method for quickly and accurately grasping the status of vegetation has great significance in terms of disaster prevention. There is. However, optical sensors cannot avoid the effects of clouds and haze, and it is difficult to accurately grasp the vegetation distribution and phenology. Especially in the tropical region, cloud influences may last for several months, and in the rainy season almost no useful data may be observed.
【0006】また、従来より広域での植生分布について
も数多くの研究がある(文献3−7参照)。例えば、GL
CC(Global Land Cover Char
acteristics data base)がある
が、日本にサバンナが出現するなど誤分類が見受けられ
る(文献7参照)。これも雲の影響を十分に取り除くこと
ができなかったのが原因である。Further, there have been many studies on the distribution of vegetation in a wider area than before (see References 3-7). For example, GL
CC (Global Land Cover Char
Although there is an artificial data base, misclassification can be seen such as the appearance of savanna in Japan (see Reference 7). This is also due to the fact that the influence of clouds could not be sufficiently removed.
【0007】さらに、これまでにも、高頻度観測衛星の
データから雲の影響を取り除いてNDVI(植生指数)
の季節変化を抽出する試みがなされている(文献8参
照)。例えば、解像度を落として、複数の画素の最大値
を取るものや、10日〜1ヶ月のデータの最大値を取る
ものなどがある。しかし、これらは、解像度が低かった
り、データの取得間隔が一ヶ月であったりと即時性の要
求される防災データとしては使いにくい。データの合成
期間が長くなれば、時間的な分解能が低下する。また、
これらは、植生の変動を一週間から10日の間隔で把握
する必要のある農林業分野では応用できないという課題
がある。Furthermore, until now, NDVI (vegetation index) has been removed by removing the influence of clouds from the data of high-frequency observation satellites.
Attempts have been made to extract the seasonal changes in (see Reference 8). For example, there is a method in which the resolution is reduced and the maximum value of a plurality of pixels is taken, and a method in which the maximum value of the data for 10 days to 1 month is taken. However, these are difficult to use as disaster prevention data that requires immediacy such as low resolution and data acquisition interval of one month. The longer the data synthesis period, the lower the temporal resolution. Also,
There is a problem that these cannot be applied in the agricultural and forestry field where it is necessary to grasp the fluctuation of vegetation at intervals of 1 week to 10 days.
【0008】また、Fourier展開類似の方法を用
いて、植生指数の季節変化を三角関数で表すこともなど
も行われている(文献9−14参照)。この方法にはどの
関数を使うかについて恣意性があり、様々な土地被覆を
含む広域データの処理には適していない。また、この方
法では、関数の当てはめにおいて、処理に用いる三角関
数群はあらかじめ決められており、必ずしも処理に適し
た関数群が選ばれていないという課題がある。[0008] Further, a method similar to the Fourier expansion is used to represent the seasonal change of the vegetation index by a trigonometric function (see References 9-14). This method is arbitrary about which function to use and is not suitable for processing wide area data including various land covers. In addition, this method has a problem that the trigonometric function group used for the processing is predetermined in the function fitting, and the function group suitable for the processing is not necessarily selected.
【0009】また、16×16などのウインドウを設定
し、その中での最大値をウインドウ内のすべての画素値
に置き換える方法もある。しかし、このような方法で
は、空間分解能が低下するという課題がある。There is also a method of setting a window such as 16 × 16 and replacing the maximum value in the window with all pixel values in the window. However, such a method has a problem that the spatial resolution is lowered.
【0010】本発明は、以上の点に鑑み、空間分解能を
保持したまま、任意の時間分解能を実現させながら、雲
やヘイズの影響を取り除くために、時系列モデルを用い
た衛星データ処理手法LMF(Local Maxim
um Fitting:局所最大値フィッティング)を
提案することを目的とする。すなわち、本発明は、画素
毎の時系列的な特徴を保存し、原データの時間分解能を
完全に保持したまま、雲などの影響およびシステムノイ
ズの除去を目的とするものである。In view of the above points, the present invention provides a satellite data processing method LMF using a time series model in order to remove the influence of clouds and haze while realizing an arbitrary time resolution while maintaining the spatial resolution. (Local Maxim
um Fitting: Local maximum fitting) is proposed. That is, the present invention aims at removing the influence of clouds and the like and system noise while preserving the time-series characteristics of each pixel and completely maintaining the time resolution of the original data.
【0011】また、本発明は、時間分解能を保持したま
ま植生などの変動の様子が捉えられるようにすること
で、防災、環境、農業分野での広域現状把握、および予
測が容易とすることを目的とする。さらに、本発明は、
適切な関数を自動的に選択し、ノイズ除去処理を自動で
行えるようにすることにより、定常業務の一部分とする
ことを可能とすることを目的とする。さらに本発明は、
水分指数など他の地球高頻度観測衛星のデータにもその
まま転用および応用を可能とすることを目的とする。Further, the present invention makes it possible to grasp the state of changes in vegetation and the like while maintaining time resolution, thereby facilitating disaster prevention, the environment, wide area current situation in the agricultural field, and prediction. To aim. Further, the present invention provides
It is an object of the present invention to enable a part of a routine work by automatically selecting an appropriate function and automatically performing noise removal processing. Further, the present invention is
The purpose is to make it possible to transfer the data to other high frequency earth observation satellite data such as the water content index and apply it as it is.
【0012】[0012]
【課題を解決するための手段】本発明の特徴のひとつと
しては、例えば、処理対象画素の時系列データについ
て、前後n時期について、局所最大値を抽出するフィル
ターを作用させる。そこで得られた局所最大値群に関数
当てはめを行い、画素値を推定する。ある時期に、推定
値と実際の画素値のずれが、設定したしきい値よりも大
きい場合には雲やシステムノイズの影響とみなしてその
時期のデータは採用しない。採用されたデータだけを用
いて再度、局所最大値の算出と関数当てはめを行った値
で画素値を置き換える。必要であればこの処理を繰り返
し、結果が収束するまで行う。As one of the features of the present invention, for example, a filter for extracting a local maximum value in time series data of a pixel to be processed is extracted for n times before and after. A function is applied to the obtained local maximum value group to estimate the pixel value. If the difference between the estimated value and the actual pixel value is larger than the set threshold value at a certain time, it is regarded as the influence of clouds or system noise and the data at that time is not used. The pixel value is replaced with the value obtained by calculating the local maximum value and performing the function fitting again using only the adopted data. If necessary, repeat this process until the results converge.
【0013】本発明の他の特徴としては、例えば、画素
毎の時系列的な局所最大値を求め、それらに対して関数
の当てはめをくりかえし行うことで、有効なデータだけ
を用いた推定を行うものである。上記の処理を画像全体
に対して行うため、雲やシステムノイズの影響が除去さ
れるという作用を有する。各画素毎の処理は独立してい
るため、原データの空間分解能は完全に保持される。As another feature of the present invention, for example, a local maximum value in time series for each pixel is obtained, and a function is repeatedly applied to them, thereby performing estimation using only effective data. It is a thing. Since the above processing is performed on the entire image, it has an effect of removing the influence of clouds and system noise. Since the processing for each pixel is independent, the spatial resolution of the original data is completely retained.
【0014】また、本発明は、時期k近傍の局所的最大
値を、時間遅れを少なくするように抽出する作用を有す
る。これにより、雲などの影響を受けた時期の画素値が
除去されるという作用を有する。Further, the present invention has the function of extracting the local maximum value near the time k so as to reduce the time delay. This has the effect of removing the pixel value at the time of being affected by clouds or the like.
【0015】本発明の第1の解決手段によると、処理部
は、観測対象エリア内の各画素に対応して、衛星による
植生の季節変動についての時系列の第1の観測データを
記憶した第1の観測データファイルから、画素毎の第1
の観測データを読み込むステップと、According to the first solution of the present invention, the processing unit stores the first observation data in time series regarding the seasonal variation of vegetation by the satellite, corresponding to each pixel in the observation target area. 1 observation data file for each pixel
Reading the observation data of
【0016】処理部は、各画素毎に、ある時期又は時点
に対応して、該時期又は時点の前後所定期間の局所最大
値を抽出することにより時系列の第2の観測データを求
め、それを各画素に対応して第2の観測データファイル
に記憶するステップと、The processing unit obtains the second observation data in time series by extracting the local maximum value in a predetermined period before and after the time or time for each pixel, corresponding to the time or time. Storing in the second observation data file for each pixel,
【0017】処理部は、第2の観測データファイルに記
憶された第2の観測データを、長期的非周期項及び季節
変化項を含む次式の関数を用い、各第1の係数ciを最
小二乗法又は他の推定法により決定することでフィッテ
ィングし、各第1の係数ciを第1のパラメータファイ
ルに記憶するステップと、The processing section uses the second observation data stored in the second observation data file as a function of the following equation including a long-term aperiodic term and a seasonal variation term, and calculates each first coefficient c i . Fitting by determining by least squares or other estimation method and storing each first coefficient c i in a first parameter file;
【0018】[0018]
【数10】
[ ここで、
ci:求める第1の係数
t:時期又は時間
N:周期関数のペア数(この場合sinとcosのペア
数):
l:周期関数の仮の番号
kl:12ヶ月が回帰の整数倍になる数列(例:1,
2,3,4,6,12)
M:1年間のデータ数(シーン数) ][Equation 10] [Where, c i : first coefficient to be obtained t: time or time N: number of pairs of periodic function (in this case, number of pairs of sin and cos): l: temporary number of periodic function k l : 12 months regression A sequence that is an integer multiple of (example: 1,
2, 3, 4, 6, 12) M: Number of data (number of scenes) in one year]
【0019】処理部は、第1の係数ciで定められた前
記関数に従い、第2の観測データから雑音データを除去
した第1の推定観測データを求め、それを第1の推定観
測データファイルに記憶する、又は出力部に出力、若し
くは表示するステップとを含む地球観測衛星データのノ
イズ除去方法、これらのステップをコンピュータに実行
させるためのノイズ除去処理プログラム、及びノイズ除
去処理プログラムを記録した記録媒体が提供される。The processing unit obtains first estimated observation data obtained by removing noise data from the second observation data according to the function defined by the first coefficient c i , and obtains the first estimated observation data file. A method for removing noise from earth observation satellite data, including a step of storing in, or outputting or displaying in an output unit, a noise removal processing program for causing a computer to execute these steps, and a record recording the noise removal processing program Media is provided.
【0020】[0020]
【発明の実施の形態】本実施の形態では、植生の季節変
動観測のために、時系列モデルを用いた新しい衛星デー
タ処理法を提案する。また、雲の影響を取り除いた植生
の季節変動に関する観測データから、火災危険度を算出
した。以下、火災危険度評価の基礎となる植生の季節変
動観測のための衛星データ処理法について述べる。な
お、植生の季節変動観測データとしては、例えば、植生
指数データ、水分指数、中間赤外線バンド、熱赤外線輝
度バンド、温度データ等である。また、観測データに含
まれる可能性がある雑音データは、例えば、雲、霧、
雨、ヘイズ等の大気状態に起因するもの、センサ感度等
のシステムに起因するもの等である。BEST MODE FOR CARRYING OUT THE INVENTION In this embodiment, a new satellite data processing method using a time series model is proposed for observing seasonal variation of vegetation. In addition, the fire risk was calculated from the observational data on the seasonal variation of vegetation excluding the influence of clouds. The satellite data processing method for seasonal observation of vegetation, which is the basis of fire risk assessment, is described below. The vegetation seasonal variation observation data includes, for example, vegetation index data, moisture index, mid-infrared band, thermal infrared brightness band, temperature data and the like. Noise data that may be included in the observation data includes, for example, clouds, fog,
Those caused by atmospheric conditions such as rain and haze, those caused by the system such as sensor sensitivity, and the like.
【0021】1. LMFの理論
LMF処理は、時系列フィルタリングと関数によるフィ
ッティング手法を組み合わせたアルゴリズムで、対象画
素の時系列データから、雲やヘイズなどの影響を除去
し、季節変化パターンを抽出するものである。フィッテ
ィングに用いる関数を、(1)式に示す。1. Theoretical LMF processing of LMF is an algorithm that combines time-series filtering and a fitting method using a function, and removes the influence of clouds, haze, etc. from the time-series data of the target pixel, and extracts the seasonal change pattern. The function used for fitting is shown in equation (1).
【0022】[0022]
【数11】
ここで、最初の二項は長期的非周期項、その他の項は季
節変化項である。また、各パラメータは、以下の通りで
ある。
ci:求めるパラメータ:c0は平均、c1はトレンド
係数、それ以外は周期関数のパラメータとなる。
t:利用するデータの時間間隔に基づく時間単位数(1
0日合成のNOAAデータを用いる場合なら、10日が
1単位)
N:周期関数のペア数(この場合sinとcosのペア
数):このソフトウェアの構造は周期関数の種類には依
存しないため、任意の周期関数と置き換えて利用するこ
とが可能。
l:周期関数の仮の番号
kl:12ヶ月が回帰の整数倍になる数列:1,2,
3,4,6,12。これによって、1ヶ月周期、2ヶ月
周期、3ヶ月周期、4ヶ月周期、6ヶ月周期および12
ヶ月周期の周期関数(sin、cos)が当てはめられ
る。
M:1年間のデータ数:10日間最大値合成のNOAA
データを用いる場合、1年間36シーンになるため、M
=36である。
なお、この場合、シーン数Nと1年間のデータ数Mが一
致しないので、一般的な離散 Fourier級数は
計算できない。[Equation 11] Here, the first two terms are long-term aperiodic terms, and the other terms are seasonal variation terms. The parameters are as follows. c i : Parameters to be obtained: c 0 is the average, c 1 is the trend coefficient, and other parameters are parameters of the periodic function. t: the number of time units (1 based on the time interval of the data to be used
When using NOAA data of 0-day synthesis, 10 days are 1 unit. N: Number of pairs of periodic functions (in this case, number of pairs of sin and cos): Since the structure of this software does not depend on the type of periodic functions, It can be used by replacing it with any periodic function. l: Temporary number of periodic function k l : Sequence in which 12 months is an integer multiple of regression: 1, 2,
3, 4, 6, 12. As a result, 1-month cycle, 2-month cycle, 3-month cycle, 4-month cycle, 6-month cycle and 12-month cycle
A periodic function (sin, cos) of the monthly cycle is fitted. M: Number of data per year: 10 days NOAA with maximum value synthesis
When using the data, since there are 36 scenes per year, M
= 36. In this case, since the number of scenes N does not match the number of data M for one year, a general discrete Fourier series cannot be calculated.
【0023】1.1. シミュレーションによる処理特
性の評価(文献15参照)
植生指数や熱赤外線輝度温度データでは、観測データに
は雲などの影響がすでに含まれており、真値は不明であ
る。そこで、本実施の形態では、シミュレーション的手
法を用いてフィルタ特性の評価を行った。まず、図1
に、シミュレーションによるフィルタ特性の評価図を示
す。例えば、(1)式の係数の組cに、発生した乱数
(S11)を代入して、時系列データを10000組発
生させる(S12)。解析のパラメータは、一例とし
て、N=54、M=36とした。次に、雲やヘイズの影
響を受けた状況をシミュレートするため、ノイズを発生
し(S13)、そのノイズを元の値の10%〜70%の
ものに置き換えて、ノイズが加わった状況をシミュレー
トした(S14)。ノイズが加算される確率は各時期に
ついて0.5とした。そして、ノイズ加算後の時系列デ
ータに対して、ノイズ除去処理し(S15)、その結果
を調べた(S16)。結果の解析は発生させた係数の組
cと、処理後に再生された係数の組c’とについて、単
回帰分析を行った。1.1. Evaluation of processing characteristics by simulation (Ref. 15) In the vegetation index and thermal infrared brightness temperature data, the observation data already include the effects of clouds and the true value is unknown. Therefore, in the present embodiment, the filter characteristics are evaluated using a simulation method. First, Fig. 1
Figure 3 shows an evaluation diagram of the filter characteristics by simulation. For example, the generated random number (S11) is substituted into the coefficient c of the equation (1) to generate 10,000 time-series data (S12). As an example, the analysis parameters are N = 54 and M = 36. Next, in order to simulate a situation affected by clouds and haze, noise is generated (S13), and the noise is replaced by 10% to 70% of the original value to determine the situation in which noise is added. It was simulated (S14). The probability that noise was added was 0.5 for each period. Then, noise removal processing was performed on the time-series data after noise addition (S15), and the result was examined (S16). In the analysis of the results, single regression analysis was performed on the generated coefficient set c and the coefficient set c ′ regenerated after the treatment.
【0024】1.2. 最小二乗法によるフィッティン
グ
一般に、ノイズ加算後のデータをそのまま(1)式でフ
ィッティングした場合、どの項も係数の約7割が保存さ
れる。これに、(2)式で表されるような局所最大値を
抽出するような時系列フィルタを作用させる。なお、こ
の局所最大値の計算式は一例であり、適宜の式を用いる
こともできる。1.2. Fitting by the method of least squares In general, when the data after noise addition is fitted as it is by the equation (1), about 70% of the coefficients are stored for each term. A time-series filter that extracts a local maximum value represented by equation (2) is applied to this. The formula for calculating the local maximum value is an example, and an appropriate formula may be used.
【0025】[0025]
【数12】
ここで、
dt’ は時間単位tの時の局所最大値 (第2の観測
データ)
dt は時間単位tの時のデータ値 (第1の観測
データ)
w は局所最大値を求める片側時間単位数(ウインドウ
サイズ)
である。[Equation 12] Here, d t 'is the local maximum value at the time unit t (second observation data) d t is the data value at the time unit t (first observation data) w is the one-sided time for which the local maximum value is obtained The number of units (window size).
【0026】(2)式では、dtにおける左右w間のデ
ータの最大値をもとめ、その小さい方を局所最大値とし
ている。この考え方が本発明のひとつの特徴である。本
実施の形態では高頻度観測衛星からの観測データd
t(第1の観測データ)に対して、(2)式のようなフ
ィルタリング処理を実行して局所最大値dt’(第2の
観測データ)を求め、さらに、(1)式の関数で最小二
乗法によりフィッティングして各係数ciを求める。こ
れによって、データに影響を与える雲やシステムノイズ
を軽減して、季節変化する地表の情報として時間単位t
に近い時期のデータを自動抽出することができる。In the equation (2), the maximum value of the data between the left and right w at d t is obtained, and the smaller one is set as the local maximum value. This idea is one of the features of the present invention. In the present embodiment, the observation data d from the high frequency observation satellite
For t (first observation data), a filtering process such as equation (2) is executed to obtain a local maximum value d t ′ (second observation data), and further, with the function of equation (1) Each coefficient c i is obtained by fitting by the method of least squares. This reduces clouds and system noise that affect the data, and the time unit t can be used as seasonally changing surface information.
It is possible to automatically extract the data in the period close to.
【0027】1.3. フィッティング関数の自動選択
上述のような式(2)によるフィルタと最小二乗法によ
るフィッティングを併用したものでは、定数項、一年周
期、半年周期項はほぼ回復できるが、それら以外の短周
期関数項の係数については保存されない場合がある。こ
れを改善するため、本発明は、さらに以下のように、
(1)式の関数のなかで、必要なもの以外はフィッティ
ングに用いないようにしてもよい。なお、最小二乗近似
に対しても、適宜のフィッティング方法を用いることが
できる。1.3. Automatic selection of fitting function With the combination of the filter by the equation (2) and the fitting by the least squares method, the constant term, the one-year period, and the half-year period term can be almost recovered, but other short-term function terms. The coefficient of may not be saved. In order to improve this, the present invention further includes the following:
Of the functions of equation (1), only the necessary functions may be used for fitting. An appropriate fitting method can be used for the least-squares approximation.
【0028】いま、最小二乗近似を用いて時系列データ
を関数でフィッティングすることを考える。フィッティ
ングに用いる関数の個数(種類)が多ければ残差は減少
するが、フィッティングした結果は不安定(結果がデー
タに大きく依存する)となる。赤池らによると、赤池の
情報量基準(AIC)が最小となるようなモデルが最適
であるとされている(文献16参照)。残差が正規分布に
従うと仮定されているから、残差の分散をσとするとA
ICは、対数尤度とパラメータ数の差で表される(3)
式のようになる。Now, consider fitting the time series data with a function using the least squares approximation. If the number of functions (types) used for fitting is large, the residual decreases, but the result of fitting becomes unstable (the result largely depends on the data). According to Akaike et al., The model that minimizes Akaike's information criterion (AIC) is the most suitable (see Reference 16). Since the residuals are assumed to follow a normal distribution, let σ be the variance of the residuals.
IC is represented by the difference between the log likelihood and the number of parameters (3).
It becomes like a formula.
【0029】[0029]
【数13】
となる。ただし、
N は所定期間のデータ数
j は用いる関数の個数(周期関数のペア数)
σはオリジナルデータとフィッティングによるデータの
残差の分散
本実施の形態ではAICの値が最小になるように、周期関
数の組み合わせを自動的に選択するようにした。即ち、
この値が最小になるようにjを選ぶようにした。なお、
正の定数倍は結果を変えないため問題ない。この例にお
いては、シミュレーションの結果、三ヶ月周期より短い
周期の項については改善がみられたが、逆に四ヶ月周期
の項についての相関係数は低くなる。[Equation 13] Becomes However, N is the number of data in a predetermined period j is the number of functions to be used (the number of pairs of periodic functions) σ is the variance of the residuals of the original data and the data due to fitting. In the present embodiment, the value of AIC is minimized. The combination of periodic functions is automatically selected. That is,
I chose j so that this value is the minimum. In addition,
There is no problem with multiplying the positive constant because it does not change the result. In this example, as a result of the simulation, an improvement was found for the term having a cycle shorter than the three-month cycle, but conversely, the correlation coefficient for the term having a four-month cycle was low.
【0030】図2は、シミュレーション結果の例を示す
図である。
1.4. フィッティング処理の繰り返しによる高精度
化
本実施の形態ではさらに、処理の高精度化を実現するた
め、フィッティングを繰り返すアルゴリズムを採用し
た。最初に、原データに局所最大値抽出フィルタを作用
させた後、フィッティングする。その結果(計算値)
と、原データとの差をとり、しきい値よりも大きなもの
を取り除く。データレンジの15〜20%をしきい値と
すれば、1年および半年に相当する周期を持つ関数の係
数は95%以上、それ以外の短周期成分も80%が回復
できることが明らかになった(文献15参照)。FIG. 2 is a diagram showing an example of the simulation result. 1.4. Improvement of accuracy by repeating fitting processing In the present embodiment, an algorithm for repeating fitting is adopted in order to further improve accuracy of processing. First, the local maximum value extraction filter is applied to the original data, and then fitting is performed. Result (calculated value)
And the difference from the original data are taken, and those larger than the threshold value are removed. It was revealed that if the threshold is set to 15 to 20% of the data range, the coefficient of the function having a period corresponding to one year and half a year can be recovered to 95% or more, and the other short period components can be recovered to 80%. (See Reference 15).
【0031】1.5. 熱赤外バンドデータの処理
NOAA/AVHRRの熱赤外バンドでは、雲のあるピ
クセルは温度が低くなる性質を用いれば、NDVIと同
じ計算方法を用いればよいと考えられる。したがって、
本発明は輝度温度についても処理が可能である。1.5. Processing of Thermal Infrared Band Data In the thermal infrared band of NOAA / AVHRR, it is considered that the same calculation method as that of NDVI may be used if the pixel with a cloud has the property that the temperature becomes low. Therefore,
The present invention can also process brightness temperatures.
【0032】2. カルマンフィルタを用いたLMFア
ルゴリズムの高精度化
2.1. LMF処理の問題点
LMFは高頻度観測衛星からの植生指数や輝度温度デー
タから雲の影響を軽減し、季節変化パターンを抽出する
のに有効な手法である。しかし、処理期間すべてにわた
ってフィッティング処理を行うため、どうしても平均化
された結果が得られてしまう。また、データ収集期間が
短い(1年程度)場合には、フィッティングの結果が不
安定になりやすく、年々変動の抽出は困難な場合があっ
た。そこで、本実施の形態ではさらに、データ処理に状
態空間モデルを採用し、データ処理を行う手法を開発し
た(文献22参照)。2. Improvement of accuracy of LMF algorithm using Kalman filter 2.1. Problems of LMF processing LMF is an effective method to reduce the influence of clouds from the vegetation index and brightness temperature data from high-frequency observation satellites and to extract seasonal change patterns. However, since the fitting process is performed over the entire processing period, an averaged result is inevitably obtained. Further, when the data collection period is short (about one year), the fitting result is likely to be unstable, and it may be difficult to extract yearly variations. Therefore, in the present embodiment, a state space model is further adopted for data processing, and a method of performing data processing is developed (see Reference 22).
【0033】また、データを防災に役立てようとするな
らば、単に現況を把握するだけでなく、短期的にせよ状
況の予測をしながら対策を講じなければならない。した
がって、時系列データの解析についても、予測が可能な
手法を採用する必要がある。状態空間モデルであれば、
この点も問題なく適用することができる。本実施の形態
では、一例として、状態空間モデルのパラメータの推定
には、Kalmanの提唱した手法であるかルマンフィ
ルタ(Kalman Filtering)を応用する
が、それに限らず、適宜の推定法を用いることができ
る。If data is to be used for disaster prevention, it is necessary to take measures not only by grasping the current situation but also by predicting the situation in the short term. Therefore, it is necessary to adopt a method capable of predicting the analysis of time series data. If it is a state space model,
This point can be applied without any problem. In the present embodiment, as an example, for estimating the parameters of the state space model, the method proposed by Kalman or the Leman filter (Kalman Filtering) is applied, but not limited to this, an appropriate estimation method may be used. it can.
【0034】カルマンフィルタは、人工衛星の軌道・姿
勢の予測と制御のように物理的に定式化が可能なものだ
けでなく、降水量からの洪水予測や経済状況の時間的変
化の予想のように、直接に原因と結果の関係が定式化で
きなくても統計量さえあればモデル構築が可能な点が大
きな特徴である(文献22参照)。さらに、
○ 少数の観測値しか得られなくても複雑な数値モデル
を構築できる。
○ 最新の観測値によってパラメータ群を更新・最適化
するので、周期的でない要因も考慮できる。
などの優れた特徴を有しているので、エル・ニーニョの
ように非定常的な擾乱要素の影響も含めたモデルの構築
が可能である(文献22参照)。The Kalman filter is not limited to those that can be physically formulated, such as prediction and control of the orbit and attitude of artificial satellites, but it can also be used to predict floods from precipitation and prediction of temporal changes in economic conditions. The major feature is that even if the cause-effect relationship cannot be directly formulated, a model can be constructed with statistics only (see Reference 22). Furthermore, a complicated numerical model can be constructed even if only a small number of observation values are obtained. ○ Parameter groups are updated and optimized with the latest observation values, so factors that are not periodic can be considered. Since it has excellent features such as, it is possible to construct a model that includes the effects of non-stationary disturbance elements such as El Niño (see Reference 22).
【0035】2.2. 時系列モデルの設定
高頻度観測衛星データを用いて得られるNDVIの季節
変化プロファイルは、三角関数の組み合わせで表現でき
ることがLMF処理の結果から明らかとなった。本実施
の形態では、この知見を生かし、離散時変スペクトルの
推定技術を応用した(文献23参照)。すなわち、(1)
式の係数ckの組を状態空間モデルの状態ベクトルと考
え、その値の推定(状態推定)をカルマンフィルタによ
って行う。2.2. Time-series model setting It became clear from the results of LMF processing that the seasonal variation profile of NDVI obtained using high-frequency observation satellite data can be expressed by a combination of trigonometric functions. In the present embodiment, this knowledge is used to apply the technique for estimating the discrete time-varying spectrum (see Reference 23). That is, (1)
The set of coefficients c k of the equation is considered as the state vector of the state space model, and its value is estimated (state estimation) by the Kalman filter.
【0036】LMFでは、雲やヘイズの影響のない季節
変化プロファイルを、(1)式の形で記述できるとし
た。ここで、右辺の第二項が非周期的変動を考慮した項
である。この式での各係数ckは、(1)式では固定で
あるが、これを時変として、季節変化プロファイルの変
動に対応できるようにした。右辺第一項が時変になった
ことで、トレンドも表現できることから、第二項は削除
できる。したがって、(4)式となる。In the LMF, the seasonal change profile without the influence of clouds and haze can be described in the form of equation (1). Here, the second term on the right side is a term in which aperiodic fluctuation is considered. Although each coefficient c k in this equation is fixed in the equation (1), it is set to be time-varying so as to be able to cope with the variation of the seasonal change profile. Since the first term on the right side is time-varying, the trend can also be expressed, so the second term can be deleted. Therefore, the equation (4) is obtained.
【0037】[0037]
【数14】
ここで、各パラメータは(1)式と同じ、ただしc
1(トレンド項)は無い。この式は先のLMFに対して、
パラメータcが全て時変変数としたものである。これに
よって、衛星が代わり、観測センサの感度が変わるなど
のシステム変動(システム雑音の一部)があっても、季
節変化を追うことができる。この点が、特に、このソフ
トウェアの画期的なところである。(4)式のc
k(t)の時間的変化は自己回帰モデルで現される。[Equation 14] Here, each parameter is the same as equation (1), but c
There is no 1 (trend term). For this LMF,
The parameters c are all time-varying variables. By this, even if there is a system change (a part of system noise) such as a change of satellite and a change in sensitivity of an observation sensor, it is possible to follow a seasonal change. This is especially the breakthrough for this software. C in equation (4)
The time change of k (t) is represented by an autoregressive model.
【0038】[0038]
【数15】
ただし、akは自己回帰係数、wKは平均値0、分散が
σkの白色雑音(システム雑音)、mは自己回帰モデル
の次数。ここでは、K=0のときのモデル、[Equation 15] Here, a k is an autoregressive coefficient, w K is a mean value 0, white noise (system noise) having a variance of σ k , and m is an order of the autoregressive model. Here, the model when K = 0,
【0039】[0039]
【数16】 は、酔歩散乱モデル[Equation 16] Is a random walk scatter model
【0040】[0040]
【数17】
と、時刻tと時刻t−1の変化率が等しいと仮定したモ
デル、[Equation 17] And a model assuming that the change rates at time t and time t−1 are equal,
【0041】[0041]
【数18】
を折衷したものを採用した。酔歩散乱モデルだけだとデ
ータの追従性が低下し、年々変動を表すことができない
場合がある。また、(7−2)式のモデルだけを採用す
ると、フィルタリング結果の値が大きく変動しすぎる場
合がある。(6−1)式を採用することによって、セン
サ感度の違いなどのノイズによるデータの急激な変化を
押さえられる。K=0以外の時のモデルは(6−2)式
である。[Equation 18] The eclectic version was adopted. If only the random walk scatter model is used, the data followability deteriorates and it may not be possible to represent year-to-year variations. Further, if only the model of the expression (7-2) is adopted, the value of the filtering result may fluctuate too much. By adopting the equation (6-1), it is possible to suppress a rapid change in data due to noise such as a difference in sensor sensitivity. The model when K = 0 is other than the equation (6-2).
【0042】[0042]
【数19】
wk(t)をシステム雑音とみなすと、観測される値y
(t)は、(8)式の状態空間モデルで表すことができ
る。[Formula 19] Observing w k (t) as system noise, the observed value y
(T) can be represented by the state space model of equation (8).
【0043】[0043]
【数20】
ここで、各パラメータは以下の通りである。
cは求めるパラメータ
Fは予め定められた行列
Λはシステム雑音が、一度しか作用しないようにするた
めの定数行列。
Hは周期関数を表す行列
v(t)は平均0、分散がσvの白色雑音(雲、ヘイズ
等の大気の状態な等に起因する観測雑音)[Equation 20] Here, each parameter is as follows. c is a parameter F to be obtained, and a predetermined matrix Λ is a constant matrix for preventing system noise from acting only once. H is a matrix representing the periodic function v (t) is white noise with mean 0 and variance σ v (observation noise caused by atmospheric conditions such as clouds and haze)
【0044】[0044]
【数21】
この点も本発明の汎用性を高める特徴である。これによ
って、例えば、異なる衛星センサを利用せざるを得ない
20年間の衛星データから環境変化を一度に分析できるよ
うになった。また、各状態変数の雑音の分散Q[Equation 21] This point is also a feature that enhances the versatility of the present invention. This forces us to use different satellite sensors, for example.
It has become possible to analyze environmental changes at once from 20 years of satellite data. Also, the noise variance Q of each state variable
【0045】[0045]
【数22】
は真値が不明であるが、この値が計算上必要であるの
で、このプログラムでは処理パラメータとして与えるこ
とにする。三角関数による時変係数推定問題では、複素
形式での表記が一般的であり、計算量も半減できるが
(文献23参照)、この例では、使用する関数を限定しな
いアルゴリズムを開発するため、敢えて複素形式にはし
なかった。[Equation 22] The true value is unknown, but since this value is necessary for calculation, it is given as a processing parameter in this program. In the time-varying coefficient estimation problem using trigonometric functions, the notation in complex form is common, and the calculation amount can be halved.
(Refer to Reference 23), in this example, the algorithm is not limited to the function to be used.
【0046】実際にはNDVIデータや熱赤外領域の輝
度温度データでは、雲の影響があれば、値が必ず低くな
る性質がある。これは、観測雑音v(t)が有色となる
可能性を含んでいる。したがって、この状態空間モデル
で求められる状態変数の値は、必ずしも最適な量とは言
えずに、準最適なものである可能性があるが、データ処
理には時系列局所最大値フィルタを併用しているので、
観測雑音の有色性は軽減されているものと考える。Actually, the NDVI data and the brightness temperature data in the thermal infrared region have a property that the value is always lowered if there is an influence of clouds. This includes the possibility that the observation noise v (t) will be colored. Therefore, the value of the state variable obtained by this state space model may not be necessarily the optimum amount and may be a sub-optimal value, but a time series local maximum filter is also used for data processing. Because
It is considered that the chromaticity of the observation noise is reduced.
【0047】2.3. カルマンフィルタ
次に、時変係数c(t)を、カルマンフィルタを用いて
推定することについて説明する。ここで、c(t|t)
を時刻tまでの観測値がわかっている場合の値とする
と、一期先の予測はc(t+1|t)となる。カルマン
フィルタは(9)式の通りである。観測値をy(t)と
する。2.3. Kalman Filter Next, estimation of the time-varying coefficient c (t) using a Kalman filter will be described. Where c (t | t)
Is the value when the observed value up to time t is known, the prediction one period ahead is c (t + 1 | t). The Kalman filter is as shown in Expression (9). Let the observed value be y (t).
【0048】[0048]
【数23】
ここで、V(t)は時刻tにおけるc(t)の共分散行
列である。これらの式を用いれば、H、F、Qのすべて
の要素があらかじめ定められているとすれば、c(t|
t−1)やV(t|t−1)が既知のとき、観測値y
(t)から他の量を計算することができる。また、y
(t)が何らかの理由で欠損になった場合には、カルマ
ンゲインが計算できないので、(10)式を用いて状態
更新を行う。[Equation 23] Here, V (t) is the covariance matrix of c (t) at time t. Using these equations, if all the elements of H, F, and Q are predetermined, then c (t |
When t-1) and V (t | t-1) are known, the observed value y
Other quantities can be calculated from (t). Also, y
If (t) becomes deficient for some reason, the Kalman gain cannot be calculated, and therefore the state is updated using equation (10).
【0049】[0049]
【数24】
ここで、システム雑音の分散σkおよび、観測雑音の分
散σvは未知なので、処理パラメータとみなす。初期値
として、[Equation 24] Here, since the system noise variance σ k and the observation noise variance σ v are unknown, they are regarded as processing parameters. As an initial value,
【0050】[0050]
【数25】
(ただし、ε>0)
ε0は正の数(0〜1)、Iは単位行列である。を用い
て最初の計算を行うことがしばしば行われる。ところ
で、カルマンフィルタのアルゴリズムを用いて、時間を
t=Tからt=1に時間を遡って処理すればデータの平
滑化(smoothing)も可能である。それらの式
は、[Equation 25] (However, ε> 0) ε 0 is a positive number (0 to 1), and I is an identity matrix. Often, the first calculation is performed using. By the way, if the Kalman filter algorithm is used to process the time backward from t = T to t = 1, it is possible to smooth the data. These formulas are
【0051】[0051]
【数26】
である。時期1からTまでのデータが取得できており、
なおかつ、カルマンフィルタにより、c(t|t)、V
(t|t)およびV(t+1|t)などの要素が計算で
きているとする。(12)式を、t=Tからtを1ずつ
減じてゆけば、t=1まで求めることができる。c(t
|T)は、得られたデータすべてを用いて推定している
ので、一般に、c(t|t)よりも精確な値であること
が知られている。ここで述べた手法をLMF−KF(L
ocal Maximum Fitting with
Kalman Filter algorithm)
と以後略記する。この方法で、ベクトルQの要素をすべ
て10-4程度にすれば、事実上、各状態変数ckは変化
しなくなるので、従来のLMFの結果を再現することも
可能である。計算量は増大するが、LMF−KFはLM
Fを完全に包含した手法であることがわかる。また、こ
のLMF−KF法の特徴として、各係数についての自己
回帰モデルである(5)式の次数mを上げれば、長期間
データ欠損があってもそれらを推定することが可能であ
り、また、そのままのアルゴリズムでデータ収集期間の
終わりからの数年間先の予測までが可能である。[Equation 26] Is. Data from period 1 to T have been acquired,
Moreover, by the Kalman filter, c (t | t), V
It is assumed that elements such as (t | t) and V (t + 1 | t) have been calculated. By subtracting t by 1 from t = T in the equation (12), it is possible to obtain up to t = 1. c (t
Since | T) is estimated using all the obtained data, it is generally known that it is a more accurate value than c (t | t). The method described here is applied to LMF-KF (L
ocal Maximum Fitting with
Kalman Filter algorithm)
Will be abbreviated below. In this method, if all the elements of the vector Q are set to about 10 −4 , each state variable c k is virtually unchanged, so that the result of the conventional LMF can be reproduced. Although the calculation amount increases, LMF-KF is LM
It can be seen that the method completely includes F. Further, as a feature of the LMF-KF method, if the order m of the autoregressive model (5) for each coefficient is increased, it is possible to estimate them even if there is data loss for a long time. , With the same algorithm, it is possible to predict several years ahead from the end of the data collection period.
【0052】2.4. 実際のデータ処理手法
画素毎に行われる、実際のデータ処理の流れを以下に説
明する。先にも触れたように、カルマンフィルタルゴリ
ズムの適用にあたって、未知量であるシステム雑音の分
散σKおよび観測雑音の分散σvはパラメータとして与
えることにするが、c(1|0)、V(1|0)も必要
である。これらの要素を画素毎に推定するのは大変なの
で、最初に(1)式を適用してLMFにより係数を求め
て初期値とする。2.4. Actual Data Processing Method The actual data processing flow performed for each pixel will be described below. As described above, in applying the Kalman filter algorithm, the variance σ K of the system noise and the variance σ v of the observation noise, which are unknown quantities, are given as parameters, but c (1 | 0), V (1 | 0) is also required. Since it is difficult to estimate these elements for each pixel, the equation (1) is first applied to obtain the coefficients by the LMF and set as the initial values.
【0053】[0053]
【数27】
(ただし、ε>0)
ε0は正の数(0〜1)、Iは単位行列である。本実施の
形態では、これらの初期値を用いてフィルタリングと平
滑化を一回行って、次回の計算の初期値とする。[Equation 27] (However, ε> 0) ε 0 is a positive number (0 to 1), and I is an identity matrix. In the present embodiment, filtering and smoothing are performed once using these initial values, and used as the initial values for the next calculation.
【0054】[0054]
【数28】
その後、フィルタリングと平滑化を繰り返し、得られた
y(t|T)がほとんど変化しなくなったところで値を
確定する。(t = 1,2,3,…. ,M(最大シ
ーン数))
このアルゴリズム(LMF−KF)は、LMFと同様
に、画素毎に処理が行われるため、LMFの時系列解析
のサブルーチンを入れ替えるだけでよい。もちろん、プ
ログラム内でのデータ構造も同じなので、並列化が容易
である。[Equation 28] After that, filtering and smoothing are repeated, and the value is determined when the obtained y (t | T) hardly changes. (T = 1, 2, 3, ..., M (maximum number of scenes)) Since this algorithm (LMF-KF) performs processing for each pixel similarly to LMF, a subroutine for time series analysis of LMF is executed. All you have to do is replace them. Of course, since the data structure in the program is the same, parallelization is easy.
【0055】3.ハードウェア
図3に、命令列スケジュール処理を実行するための装置
の構成図を示す。この装置は、処理部(CPU)1、入
力部2、出力部3、記憶部4を備える。処理部1は、C
PUであり、命令列スケジュールプログラムによりスケ
ジュールを実行する。記憶部4は、計算パラメータファ
イル41、第1の観測データファイル42、第2の観測
データファイルファイル43、第1の推定観測データァ
イル44、第2の推定観測データファイル45、第1の
パラメータファイル46、第2のパラメータファイル4
7を有する。3. Hardware FIG. 3 shows a block diagram of an apparatus for executing the instruction sequence schedule processing. This device includes a processing unit (CPU) 1, an input unit 2, an output unit 3, and a storage unit 4. The processing unit 1 is C
It is a PU and executes a schedule by an instruction sequence schedule program. The storage unit 4 includes a calculation parameter file 41, a first observation data file 42, a second observation data file 43, a first estimated observation data file 44, a second estimated observation data file 45, and a first parameter file. 46, second parameter file 4
Have 7.
【0056】4.フローチャート
図4に、実際のデータ処理のフローチャートを示す。な
お、破線枠内の処理は、画素毎に行われる。処理部1
は、入力部2又は記憶部4の計算パラメータファイル4
1から計算パラメータを読み込む(S101)。ここで
は例えば、計算パラメータとして、画像数、像ライン
数、画像コラム数、関数セット種類数、欠損シーン数、
カルマンフィルタ使用時の関数セット数、局所最大値算
出時の片側ウインドウサイズ雲検出閾値、カルマンフィ
ルタ処理時の最大繰り返し回数、カルマンフィルタ収束
判定値、システムノイズの分散パラメータなどがある。4. Flowchart FIG. 4 shows a flowchart of actual data processing. The processing within the broken line frame is performed for each pixel. Processing unit 1
Is the calculation parameter file 4 of the input unit 2 or the storage unit 4.
The calculation parameters are read from 1 (S101). Here, for example, as the calculation parameters, the number of images, the number of image lines, the number of image columns, the number of function set types, the number of missing scenes,
The number of function sets when the Kalman filter is used, the one-sided window size cloud detection threshold when calculating the local maximum value, the maximum number of iterations when the Kalman filter processing is performed, the Kalman filter convergence determination value, the dispersion parameter of the system noise, and the like.
【0057】処理部1は、衛星による画素データ毎の植
生の季節変動についての時系列な観測データを記憶した
記憶部4の第1の観測データファイル42から、画素デ
ータ毎の第1の観測データを読み込む(S103)。こ
の観測データは、例えば、植生指数データ、水分指数、
中間赤外線バンド、熱赤外線輝度バンド、温度データ等
がある。処理部1は、記憶部4の第1の観測データファ
イル42から、1ライン又は1カラムのデータを読み込
み、読み込まれたライン又はカラムのデータを並列演算
計算機のパラメータに設定してもよい。The processing unit 1 uses the first observation data file 42 of the storage unit 4 in which the time-series observation data of seasonal variation of vegetation for each pixel data by the satellite is stored to obtain the first observation data for each pixel data. Is read (S103). This observation data includes, for example, vegetation index data, moisture index,
There are mid-infrared band, thermal infrared brightness band, temperature data, etc. The processing unit 1 may read data of one line or one column from the first observation data file 42 of the storage unit 4 and set the read data of the line or column as a parameter of the parallel computing computer.
【0058】処理部1は、式(2)を用いて各画素毎
に、ある時期の前後所定期間の局所最大値を抽出し、第
2の観測データを得て、第2の観測データファイル43
に画素毎に記憶する(S105)。The processing unit 1 extracts the local maximum value in a predetermined period before and after a certain time for each pixel using the equation (2), obtains the second observation data, and the second observation data file 43.
Is stored for each pixel (S105).
【0059】処理部1は、第2の観測データファイル4
3から第2の観測データを読み出し、これに基づき式
(1)の関数にフィッティング処理を実行して各係数c
kを最小二乗法により決定し、その係数を第1パラメー
タファイル46に記憶する(S107)。なお、処理部
1は、式(3)で示される赤池の情報量基準AICが最
小となる周期関数の組み合わせを選択するステップをさ
らに含むことができる。次のような各ステップが実行さ
れる。即ち、
・処理部1は、全ての関数の組み合わせを発生させる第
1ステップと、
・処理部1は、各組み合わせに対して最小二乗法を実行
してフィッティングに用いられる係数の関数を求める第
2ステップと、
・処理部1は、最小二乗法により求めた結果について、
それぞれAICを求める第3ステップと、
・処理部1は、AICの最小の組み合わせを使用するフ
ィッティングに用いられる関数として選択する第4ステ
ップとを含む。The processing unit 1 uses the second observation data file 4
The second observation data is read from No. 3, and the fitting process is executed on the function of the equation (1) based on this to obtain each coefficient c.
k is determined by the method of least squares, and the coefficient is stored in the first parameter file 46 (S107). Note that the processing unit 1 can further include a step of selecting a combination of periodic functions that minimizes the Akaike information amount reference AIC expressed by Expression (3). The following steps are executed. That is, the processing unit 1 generates a combination of all the functions, and the processing unit 1 executes a least squares method for each combination to obtain a function of a coefficient used for fitting. Steps: The processing unit 1 calculates the result obtained by the method of least squares,
Each includes a third step of obtaining an AIC, and the processing unit 1 includes a fourth step of selecting a minimum combination of AICs as a function used for fitting.
【0060】また、処理部1は例えば、ステップS10
5の前又は後若しくはステップS107の前又は後等
に、欠測値の補間をしてもよい。さらに欠測値の補間後
に、補間修正をしてもよい。さらに、処理部1は、原デ
ータとフィッティングされた関数による計算値の差を取
り、差がしきい値よりも大きなものを取り除いても良
い。また、処理部1は採用すべきデータのみを用いてフ
ィッティング、フィルタリング(カルマンフィルタ)等
の各処理を繰り返すことができる。Further, the processing section 1 is, for example, step S10.
The missing value may be interpolated before or after step 5, or before or after step S107. Further, the interpolation correction may be performed after the interpolation of the missing value. Further, the processing unit 1 may take the difference between the original data and the calculated value by the fitted function, and remove the one whose difference is larger than the threshold value. Further, the processing unit 1 can repeat each processing such as fitting and filtering (Kalman filter) using only the data to be adopted.
【0061】つぎに、処理部1は、定められた係数ck
に従い、関数によりノイズが除去された推定値を第1推
定観測データファイル44に記憶する、又は出力部3に
出力若しくは表示する。また、処理部1は、第1パラメ
ータファイル46から係数を読み出し、カルマンフィル
タ処理の初期値を式(13)のように設定する(S10
9)。Next, the processing section 1 determines the determined coefficient c k.
Accordingly, the estimated value from which the noise is removed by the function is stored in the first estimated observation data file 44, or is output or displayed on the output unit 3. Further, the processing unit 1 reads the coefficient from the first parameter file 46 and sets the initial value of the Kalman filter processing as shown in Expression (13) (S10).
9).
【0062】処理部1は、y(t|T)old とKの
値を、ぞれぞれ初期値として0(y(t|T)old=
0,K=0)とする(S111)。つぎに、処理部1
は、式(9)に基づいて、カルマンフィルタ処理を行う
(S113)。さらに、処理部1は、式(12)を用い
て、データの平滑化する(S115)。The processing unit 1 sets the values of y (t | T) old and K as 0 (y (t | T) old =
0, K = 0) (S111). Next, the processing unit 1
Performs Kalman filter processing based on equation (9) (S113). Further, the processing unit 1 smoothes the data by using the equation (12) (S115).
【0063】処理部1は、Kの値を1繰り上げる(K=
K+1) (S117)。さらに、処理部1は、時期K
に、推定値と実際の画素値のずれを求め、この値によっ
て次に行う処理を判断する。即ち、カルマンフィルタ処
理が収束してなく、このずれが予め設定したしきい値ε
よりも大きい場合には、次にステップS121の処理を
行い、一方収束して、しきい値εよりも小さい場合には
次にステップS125の処理を行う(S119)。The processing unit 1 increments the value of K by 1 (K =
K + 1) (S117). Furthermore, the processing unit 1 determines the time K
Then, the difference between the estimated value and the actual pixel value is obtained, and the next process is determined based on this value. That is, the Kalman filter processing has not converged, and this deviation is due to a preset threshold value ε.
If it is larger than the threshold, the process of step S121 is performed next, while it converges, and if it is smaller than the threshold ε, then the process of step S125 is performed (S119).
【0064】即ち、収束していない場合、処理部1は、
計算パラメータc(1|0)とV(1|0)の初期値を
更新し(c(1|0)=c(1|T)、V(1|0)=
V(1|T))(S121)、さらに、y(t|T)
oldの値をy(t|T)ol d=y(t|T)とする
(S123)。一方、収束した場合には、処理部1は、
計算結果である推定された係数ci(t)を第2のパラ
メータファイル47に記憶し、ノイズを除去した推定さ
れた観測値を第2の推定観測データファイル45に記憶
する、又は、出力部3に出力若しくは表示する(S12
5)。That is, when the processing has not converged, the processing unit 1
The initial values of the calculation parameters c (1 | 0) and V (1 | 0) are updated (c (1 | 0) = c (1 | T), V (1 | 0) =
V (1 | T)) (S121), and further y (t | T)
the value of the old y (t | T) ol d = y | a (t T) (S123). On the other hand, if converged, the processing unit 1
The estimated coefficient c i (t), which is the calculation result, is stored in the second parameter file 47, and the estimated observation value with noise removed is stored in the second estimated observation data file 45, or the output unit Output to 3 or display (S12
5).
【0065】5. 実際の衛星データ処理例
LMF処理を行うためには、地理補正、センサー補正済
みであることが必要である。すべてのシーンで、それぞ
れの画素の地理座標が揃っている必要がある。そのよう
なデータは、農林水産研究計算センター衛星データアー
カイブシステムであるSIDaB(文献17,18参
照)、Pathfinder プロジェクト(文献19参
照)、USGS(文献20参照)などのwebサイトから
ネットワークを通じて入手できる。また、SPOT/v
egetationであれば、切り出し範囲と投影法が
揃っていれば、ほとんど前処理なしでLMF処理するこ
とが可能である。5. Example of actual satellite data processing In order to perform LMF processing, it is necessary to have completed geographic correction and sensor correction. All scenes must have the same geographic coordinates for each pixel. Such data can be obtained through a network from websites such as SIDaB (see References 17 and 18) which is a satellite data archiving system for agriculture, forestry and fisheries research and computation, Pathfinder project (see Reference 19), and USGS (see Reference 20). Also, SPOT / v
In the case of the acquisition, if the cutout range and the projection method are the same, the LMF processing can be performed with almost no preprocessing.
【0066】計算上、LMF処理を行うには、最低1年
分(10日間合成で36シーン)が必要である。特に雲
の影響が強い熱帯では、1年半(54シーン)以下で
は、結果が期待できない。3年以上のデータの蓄積があ
れば、一ヶ月程度のシーン欠損があってもほぼ問題なく
処理することが可能である。計算はかなり膨大なものに
なるので、並列機などの使用が望ましい。(並列機を用
いての処理については付記を参照のこと)From the calculation, at least one year (36 scenes for 10 days of composition) is required to perform the LMF processing. Especially in the tropics where the influence of clouds is strong, the results cannot be expected after one and a half years (54 scenes). If data is accumulated for three years or more, even if there is a scene loss for about a month, it can be processed with almost no problem. Since the calculation becomes enormous, it is desirable to use parallel machines. (Refer to the appendix for processing using a parallel machine.)
【0067】5.1. Pathfinder データ
図5に、LMF処理前後のNDVIデータ図を示す。こ
の図は、Pathfinderデータ(8km NOA
A/AVHRR 10日間合成)を用いてNDVIのL
MF処理を行った例である。対象画素としては、北緯3
6.31度、東経136.92度(岐阜県大野郡白川
村)を選んだ。データ収集期間は95年1月上旬〜99
年12月下旬の180シーンを用いた。5.1. Pathfinder data FIG. 5 shows NDVI data diagrams before and after the LMF process. This figure shows Pathfinder data (8km NOA
A / AVHRR 10 days synthesis)
This is an example of performing MF processing. The target pixel is north latitude 3
We chose 6.31 degrees and 136.92 degrees east longitude (Shirakawa-mura, Ono-gun, Gifu prefecture). Data collection period is from early January 1995 to 99
180 scenes in late December were used.
【0068】Pathfinder データには、画素
毎に、データ取得時の雲の有無についてのフラグ(CL
AVR データ)が添付されている。CLAVRの値
が、1〜11 は“cloudy”、12〜22は”m
ixed”、23〜31は”clear”にカテゴリが
分かれている(文献21参照)。雪の影響の無い、95年
5月〜9月のNDVIデータについて、日本付近を例に
取り、LMF処理の精度について検証した。In the Pathfinder data, a flag (CL
(AVR data) is attached. The value of CLAVR is “cloudy” for 1 to 11 and “m” for 12 to 22.
"ixed" and 23-31 are divided into "clear" categories (see Reference 21.) For the NDVI data of May-September 1995, which is not affected by snow, taking the vicinity of Japan as an example, the LMF processing The accuracy was verified.
【0069】図6に、LMF処理前後のデータの相関図
を示す。この図は、処理データのNDVIを横軸にと
り、LMF処理済みデータを縦軸としたものである。処
理前のNDVIについて0.01毎に区間を区切り、そ
の区間に属する同じ時期、同一画素のLMF処理後のN
DVIの平均値と標準偏差をプロットした。雲の影響の
ない“Clear”カテゴリに属する画素では、LMF
処理前後の値に高い相関が見られるのに対し、部分的に
雲が混ざった“Mixed”や雲に覆われた”Clou
dy”カテゴリに属する画素ではLMF処理後の値が高
くなっている。FIG. 6 shows a correlation diagram of data before and after the LMF processing. In this figure, the NDVI of the processed data is plotted on the horizontal axis, and the LMF-processed data is plotted on the vertical axis. Sections are separated by 0.01 for NDVI before processing, and N after the LMF processing of the same pixel at the same time belonging to the section.
The average value and standard deviation of DVI were plotted. For pixels belonging to the “Clear” category, which is not affected by clouds, LMF is used.
There is a high correlation between the values before and after the treatment, while "Mixed" where the clouds are partially mixed and "Clou covered by the clouds".
The pixel after the LMF process has a high value in the pixels belonging to the “dy” category.
【0070】雲に覆われた画素のNDVIの真値を知る
ことは不可能であるが、少なくとも、雲の影響の無い画
素については、LMF処理の前後でほぼ値が保存されて
いる。It is impossible to know the true value of the NDVI of the pixel covered by the cloud, but at least for the pixel not affected by the cloud, the values are almost saved before and after the LMF processing.
【0071】熱赤外バンドデータの処理:図7に、AV
HRR/Ch4のデータをLMF処理したデータ図を示
す。この図によると、NDVIの時と同様、雲などの影
響が軽減されていることがわかる。
最小値画像:これまで、NDVI最大値画像は作られて
きたが、NDVI最小値画像は雲の影響を受けた画素ば
かりを集めてしまうために、意味づけをすることはでき
なかった。LMF処理により雲の影響を除去することで
はじめて、意味のある最小値画像が得られたことにな
る。ここで得られたNDVI最小値画像は、年間を通じ
て植物が多く存在するか否かを示している。また、ND
VIだけでなく、熱赤外バンドの輝度温度に関しても、
LMF処理済みデータであれば、最小値の意味づけが可
能である。Processing of thermal infrared band data: AV in FIG.
The data figure which carried out the LMF process of the data of HRR / Ch4 is shown. According to this figure, it can be seen that the influence of clouds and the like is reduced as in the case of NDVI. Minimum value image: Until now, an NDVI maximum value image was created, but it could not be meaningful because the NDVI minimum value image collects only the pixels affected by clouds. A meaningful minimum value image is obtained only by removing the influence of clouds by the LMF process. The NDVI minimum value image obtained here shows whether or not there are many plants throughout the year. Also, ND
Not only VI, but also the brightness temperature of the thermal infrared band,
With LMF-processed data, the minimum value can be given meaning.
【0072】5.2. USGS 1km 全球AVH
RRデータ
米国地質調査所(USGS)にアーカイブされている1
km解像度の AVHRR のNDVIについて、タイ
周辺のデータの処理前と処理後のものを比較したとこ
ろ、タイは10月〜3月までが乾期、4月〜9月までが
雨期であるが、雨期の雲の影響を抑えて植生が繁茂する
様子がよく捉えられている。また、データの欠損や原デ
ータにおけるパス間の質の違いも吸収できていることが
わかる。また、1画素毎に、NDVIの最小値と最大
値、季節変化パターンが得られるので、植生タイプの分
類だけでなく、植生の乾燥度や植生由来の地上可燃物量
についての知見を得ることができる。5.2. USGS 1km Global AVH
RR data archived by US Geological Survey (USGS) 1
Comparing the NDVI of AVHRR with km resolution between before and after processing the data around Thailand, Thailand shows that the dry season is from October to March and the rainy season is from April to September, but the rainy season is from Thailand to October. It is often seen that vegetation grows with the influence of clouds suppressed. In addition, it can be seen that the loss of data and the difference in quality between the paths in the original data can be absorbed. In addition, since the minimum and maximum values of NDVI and the seasonal variation pattern are obtained for each pixel, it is possible to obtain not only the classification of vegetation types but also the degree of vegetation dryness and the amount of terrestrial combustibles derived from vegetation. .
【0073】5.3. 実際のデータ処理例:Path
finderデータの処理
図8は、PathfinderデータのNDVIをLM
F−KF処理した結果を示す図である。LMFとLMF
−KFの結果を比較すると、LMF処理の結果が95年
〜99年の5年間の平均的な季節変化を表しているのに
対し、LMF−KFでは、雲の影響を取り除きながら、
年々変動にも追従している。注目すべきは、96年早春
は値の立ち上がりが遅く、98年では逆に値の立ち上が
りが早いのがはっきり現れていることである。過去の気
象データによると、96年は北陸地方の積雪量は例年と
比較して多かったが、98年早春はエル・ニーニョ現象
によって暖冬となり、積雪量が異常に少なかったことが
報告されている。5.3. Actual data processing example: Path
Processing of Finder Data FIG. 8 shows the NDVI of Pathfinder data as an LM.
It is a figure which shows the result of F-KF processing. LMF and LMF
-Comparing the results of -KF, the results of LMF treatment represent the average seasonal changes for 5 years from 1995 to 1999, while LMF-KF removes the influence of clouds,
It also follows year-to-year fluctuations. It should be noted that the rise of prices is slow in early spring of 1996, and on the contrary, the rise of prices is rapid in 1998. According to past meteorological data, the amount of snow in the Hokuriku region was higher in 1996 than in the average year, but it was reported that early spring 1998 was a mild winter due to the El Niño phenomenon, and the amount of snow was unusually low. .
【0074】また、図9は、AVHRR/Ch4輝度温
度データをLMF−KF処理した結果を示す図である。
NDVIデータと同じ傾向を示しており、LMF−KF
では96年、98年の特徴がはっきりと現れている。FIG. 9 is a diagram showing the result of LMF-KF processing of AVHRR / Ch4 brightness temperature data.
It shows the same tendency as NDVI data, and LMF-KF
Then, the characteristics of 1996 and 1998 are clearly shown.
【0075】本実施の形態では、ボルネオ周辺のNDV
Iデータについて、95年〜99年の5年間のLMF−
KF処理を行い、1年毎にカラー合成をしたものを求め
たところ、95年〜99年になるにしたがって、年間を
通じてNDVIの高い場所が減少し、代わって、植生指
数の低い赤色の部分が増加している。また、季節変動の
大きいことを示す白っぽい部分も増えている。98年以
降、ボルネオ島東部に植生の少ない地域が増えている。
これは98年にエル・ニーニョ現象のため大規模森林火
災が発生したが、その延焼地域を表している。In this embodiment, the NDV around Borneo is used.
For I data, LMF for 5 years from 1995 to 1999
When KF processing was performed and color composites were obtained every year, as 195-99 years, places with high NDVI decreased throughout the year, and instead, red parts with low vegetation index were found. It has increased. In addition, the number of whitish areas showing large seasonal fluctuations is increasing. Since 1998, the number of vegetated areas in the eastern part of Borneo is increasing.
This represents the spread of a large-scale forest fire that occurred in 1998 due to the El Niño phenomenon.
【0076】LMF処理を用いれば年々変動の解析も可
能であるが、データを期間毎に切り出して個別に処理す
る必要があった。さらには、処理の性質としてデータの
収集期間が短くなると結果に雲の影響が含まれやすくな
る(モデルが不安定化する)ことや、データ収集期間の
つなぎ目での計算値の整合性がとれない場合があった。It is possible to analyze yearly fluctuations by using the LMF processing, but it was necessary to cut out the data for each period and process them individually. Furthermore, as the nature of the processing, if the data collection period is shortened, cloud effects are likely to be included in the result (the model becomes unstable), and the calculated values at the joints of the data collection period are not consistent. There were cases.
【0077】しかし、LMF−KFでは、長期間のデー
タを一度に処理しても、それぞれの年の特徴が抽出でき
ることが示された。データの収集期間が長いので、モデ
ルの不安定性が軽減されるためであると考えられる。However, it was shown that the LMF-KF can extract the characteristics of each year even if the long-term data is processed at once. This is because the instability of the model is reduced because the data collection period is long.
【0078】6. 時系列データ処理法のまとめ
今回開発した、時系列フィルタと状態空間モデルを用い
た、高頻度観測衛星データ処理法では、データ処理期間
全体の平均的な季節変化を捉えるだけでなく、雲やヘイ
ズ、システムノイズなどの影響を抑えながら、年々変動
の抽出も可能であることがわかった。6. Summary of time-series data processing method The high-frequency observation satellite data processing method developed this time, which uses a time-series filter and a state space model, not only captures the average seasonal changes over the entire data processing period, but also clouds and haze. , It was found that it is possible to extract yearly variations while suppressing the influence of system noise.
【0079】この方法により、植生指数だけでなく、熱
赤外バンドの輝度温度データについても衛星データから
季節変動プロファイルを抽出できるようになったが、こ
れにより、一週間〜10日毎の時間間隔で全球規模の植
生環境の変動を監視することがはじめて可能となった。
特に、大規模森林火災や乾燥害などの農業・林業災害は
植生の乾燥の程度に影響されるため、植生の状態を迅速
かつ的確に把握する手法の開発は防災面で大きな意味を
持っている。By this method, not only the vegetation index but also the seasonal temperature profile of the brightness temperature data of the thermal infrared band can be extracted from the satellite data. For the first time, it is possible to monitor changes in the global vegetation environment.
In particular, agricultural and forestry disasters such as large-scale forest fires and drought damage are affected by the degree of vegetation dryness, so the development of a method for quickly and accurately ascertaining the vegetation state has great significance in terms of disaster prevention. .
【0080】植生指数や輝度温度の季節変動データか
ら、森林火災危険度を算出できるが、それについては次
節で述べる。また、植生指数の季節変動パターンの分類
によって、広域での土地被覆分類も可能である。The forest fire risk can be calculated from the seasonal variation data of the vegetation index and the brightness temperature, which will be described in the next section. In addition, land cover classification in a wide area is possible by classifying seasonal variation patterns of vegetation index.
【0081】LMFおよび LMF−KFは、画素毎に
処理が独立しているため、並列機に適している。したが
って、MODISなど、従来より解像度の高い高頻度観
測衛星データについても、時間をかけずに処理が行える
利点がある。The LMF and LMF-KF are suitable for parallel machines because the processing is independent for each pixel. Therefore, there is an advantage that high-frequency observation satellite data having a higher resolution than that of the conventional one, such as MODIS, can be processed in a short time.
【0082】7. 輝度温度データを用いた火災検出
つぎに、LMFおよびLMF−KFによって処理された
データの活用、特に輝度温度データを用いた火災検出の
高精度化について言及する。
DMSP/OLSでの雲域の自動推定:7. Fire Detection Using Luminance Temperature Data Next, the utilization of data processed by LMF and LMF-KF, particularly the high accuracy of fire detection using luminance temperature data will be mentioned. Automatic estimation of cloud area in DMSP / OLS:
【0083】DMSP/OLS での雲域除去は、特定
のしきい値を定めるか、NOAA/NGDCから配信さ
れる全球地表面温度データと比較することで行ってい
る。前者の方法では、しきい値の設定に恣意性があった
り、雲除去が適切に行われないなど、精度において問題
がある。一方、配信される全球地表面温度データは、配
信時刻が衛星データと同時でなく1日程度遅れており、
リアルタイム処理には間に合わない。また、地上観測で
の地表面温度と衛星で観測される熱赤外バンドの輝度温
度は必ずしも同じ物理量ではない。The cloud area removal in DMSP / OLS is performed by setting a specific threshold value or by comparing with global surface temperature data distributed from NOAA / NGDC. The former method has a problem in accuracy such that the threshold value is arbitrary and cloud removal is not performed properly. On the other hand, the global surface temperature data that is delivered has a delivery time that is not the same as the satellite data and is delayed by about one day
It is not in time for real-time processing. In addition, the ground surface temperature in ground observation and the brightness temperature in the thermal infrared band observed by satellite are not necessarily the same physical quantity.
【0084】DMSP/OLSの熱赤外輝度温度をLM
F処理しておけば、各時期、各画素での雲の影響のない
値が得られるので、雲の判定は単純にしきい値を設定す
れば足りる。また、観測波長の問題はあるが、次節に示
す方法により、熱によっても火災検出が可能になる可能
性がある。
NOAA/AVHRR での火災検出アルゴリズムの単
純化:The thermal infrared brightness temperature of DMSP / OLS is set to LM.
If the F processing is performed, a value that does not affect the cloud at each pixel at each time can be obtained, so that it suffices to simply set the threshold value to determine the cloud. In addition, although there is a problem with the observation wavelength, the method described in the next section may enable fire detection by heat. Simplification of the fire detection algorithm in NOAA / AVHRR:
【0085】現在、NOAA/AVHRR でのホット
スポット検出は、単純に輝度温度にしきい値を設定する
ものや、バンド間演算をしたもので判断するものが多い
(文献24参照)。この方法では、どうしても、裸地や火
口、人口構造物(石油および製鉄プラント)などの火災
でないものも検出してしまう欠点がある。これらを回避
するため、複雑な火災検出アルゴリズムも開発されてい
るが24)、演算量の点から、リアルタイム処理には適
していない。さらには、土地被覆によって火災検出に用
いるしきい値だけでなく、アルゴリズムも変更する必要
があるとの指摘もなされている。これを実現することを
考えた場合、特定の狭い範囲で火災検出システムを構築
することは可能であっても、国レベルの広域を対象にす
ることは困難である。At present, hot spot detection in NOAA / AVHRR is often made by simply setting a threshold value for the brightness temperature or by performing inter-band calculation.
(See Reference 24). This method has a drawback that it can detect non-fires such as bare ground, craters, and man-made structures (oil and steelmaking plants). In order to avoid these, a complicated fire detection algorithm has been developed, but it is not suitable for real-time processing due to the amount of calculation. Furthermore, it has been pointed out that not only the threshold used for fire detection but also the algorithm must be changed depending on the land cover. Considering to realize this, although it is possible to construct a fire detection system in a specific narrow range, it is difficult to target a wide area at the national level.
【0086】ここで、LMF処理済みの、雲無しの輝度
温度データがあれば、それらは、対象となる画素の土地
被覆や緯度、標高などの影響がすでに含まれているの
で、単純にしきい値を設定するだけでも火災精度の向上
が期待できる。また、次節以下で述べる方法により、画
素内の火災割合を算出することも可能となる。
衛星データによる小規模火災の検出:Here, if there is cloud-free brightness temperature data that has been LMF-processed, since they already include the influences of the land cover, latitude, and elevation of the target pixel, they are simply thresholded. You can expect an improvement in fire accuracy simply by setting. It is also possible to calculate the fire rate within a pixel by the method described in the next section. Small fire detection with satellite data:
【0087】Prince らが、静止気象衛星 GO
ES を用いてブラジルの森林火災監視を行っている
(文献25参照)。一般に、静止気象衛星は高度が非常に
高いため、解像度は低くなってしまうが、画素内の火災
割合を算出することで、NOAA/AVHRR と同程
度の火災検知を可能にしている。ここで、熱赤外バンド
iでの、温度Tの黒体の分光放射輝度をLi(T)、画
素内の火災割合をpとすると、
Li(Tobs)=p・Li(Tfire)+(1−p)・Li(Tb) (15)
の式が成り立つ。ただし、Tobsは観測されたものの
温度、Tfireは火災温度、Tbは画素内の火災のな
い部分の温度である。この式には、未知数がp、T
fire、Tbの3個あるので、熱赤外バンドは3つ必
要である。しかし、Princeらは、Tfire=4
00(K)、Tbを周囲の比較的低温な画素の平均と仮
定して火災割合を算出している。しかし、これでは、火
災のない場合、対象となる画素と、その周囲の画素の温
度が同じであることを仮定している。換言すれば、対象
となる画素と、その周囲の画素の土地被覆や標高が同一
であることになるが、この仮定が成り立つ場所はかなり
限られてくる。ここで、Tbとして、対象となる画素の
同時期のLMF処理された値を使えば、DMSP/OL
Sのように熱赤外バンドが一つしかなくても、火災割合
pを算出することが可能になる。また、NOAA/AV
HRRのように、複数の熱赤外バンドがあれば、火災割
合と火災温度の両方を算出することが可能である。Prince et al., Geostationary Meteorological Satellite GO
Monitoring forest fires in Brazil using ES
(Ref. 25). In general, the geostationary meteorological satellite has a very high altitude and thus has a low resolution. However, by calculating the ratio of fires in pixels, it is possible to detect fires at the same level as NOAA / AVHRR. Let L i (T) be the spectral radiance of the black body at temperature T in the thermal infrared band i, and p be the fire ratio in the pixel, then L i (T obs ) = p · L i (T fire ) + (1−p) · L i (T b ) (15) However, T obs is the observed temperature, T fire is the fire temperature, and T b is the temperature of the non-fire part in the pixel. In this equation, the unknowns are p, T
Since there are three fire and T b , three thermal infrared bands are required. However, Prince et al., T fire = 4
The fire rate is calculated by assuming that 00 (K) and T b are the average of relatively low temperature surrounding pixels. However, this assumes that the temperature of the target pixel and the surrounding pixels are the same when there is no fire. In other words, the target pixel and the surrounding pixels have the same land cover and elevation, but the places where this assumption holds are quite limited. Here, if the value of the target pixel subjected to LMF processing at the same time is used as T b , DMSP / OL
Even if there is only one thermal infrared band like S, the fire rate p can be calculated. In addition, NOAA / AV
If there are a plurality of thermal infrared bands like HRR, it is possible to calculate both the fire rate and the fire temperature.
【0088】NOAA/AVHRRでは、センサーが飽
和してしまう温度が低いため、火災温度の推定は困難で
あるが、MODIS等の新型センサーでは、火災温度に
よる消火優先度の設定なども考えられる。この方法は、
GMS−5の後継機(MTSAT)のデータにも応用が
可能であり、毎時の火災情報提供への応用が期待でき
る。In NOAA / AVHRR, it is difficult to estimate the fire temperature because the temperature at which the sensor saturates is low, but in new sensors such as MODIS, it is possible to set the fire extinguishing priority depending on the fire temperature. This method
It can be applied to the data of the successor to GMS-5 (MTSAT), and can be expected to be applied to the provision of fire information every hour.
【0089】8.付記及び参考文献
本発明の地球観測衛星データのノイズ除去方法又は地球
観測衛星データのノイズ除去装置・システムは、その各
手順をコンピュータに実行させるための地球観測衛星デ
ータのノイズ除去処理プログラム、地球観測衛星データ
のノイズ除去処理プログラムを記録したコンピュータ読
み取り可能な記録媒体、地球観測衛星データのノイズ除
去処理プログラムを含みコンピュータの内部メモリにロ
ード可能なプログラム製品、そのプログラムを含むサー
バ等のコンピュータ、等により提供されることができ
る。以下に、参考文献を示す。8. Additional Notes and References The noise removal method for earth observation satellite data or the noise removal device / system for earth observation satellite data according to the present invention is a noise removal processing program for earth observation satellite data for causing a computer to execute the respective steps, earth observation A computer-readable recording medium recording a satellite data noise removal processing program, a program product containing the earth observation satellite data noise removal processing program that can be loaded into the internal memory of the computer, a computer such as a server containing the program, etc. Can be provided. The following are references.
【0090】1) インドネシア森林火災予防プロジェク
ト、“インドネシア国立公園森林火災調査報告書”、国
際協力事業団、1998
2) 沢田治雄、“大規模森林火災へのリモートセンシン
グの利用”、森林科学、24、22-28、1998
3) 本多嘉明、村井俊治、加藤喜久雄、“世界植生モニ
タリング」、日本写真測量学会誌「写真測量とリモート
センシング」、31、1、4-14、1992
4) Defries,R. S. and J. R. G. Townshed,
“NDVI-derived land cover classifications at
a global scale”, Int. J. Remote Sens.,15,
3567-3586, 1994.
5) Myneni, R. B., C. J. Tucker, G. Asrar
and C. D. Keeling,“Interannual variations in
satellite-sensed vegetation index data from
1981 to 1991”, J. Geophys. Res., 103, D
6, 6145-6160,1998.
6) Los, S. O., C. O. Justice and C. J. T
ucker, “A global 1deg. by 1 deg. NDVI dat
a set for climate studies derived from the
GIMMS continental NDVI Data.”, Int. J. Re
mote Sens., 15, 3493-3518, 1994.
7) GLCC web site : http://edcdaac.usgs.gov/gl
cc/globdoc1_2.html
8) Viovy, N., O. Arino and A. S. Belward,
“The Best Index Slope Extraction (BISE):
A method for reducing noise in NDVI time-se
ries”, , Int. J. Remote Sens., 13, 1585-15
90, 1992.
9) 李 雲慶、大沼一彦、安田嘉純、“時系列植生指数
データのフーリエ解析”、日本写真測量学会誌「写真測
量とリモートセンシング」、31、4、4-14、1992
10) Sellers, P. J., C. J. Tucker, G. J. C
ollatz, S. O. Los,C. O. Justice, D. A. Daz
lich and D. A. Randall, “A global1 deg. b
y 1 deg. NDVI data set for climate studie
s. Part 2: The generation of global fields
of terrestrial biophysical parameters from
the NDVI”, Int. J. Remote Sens., 15, 3519-
3545,1994.
11) Verhoef, W., M. Menenti and S. Azzali,
“A colour composite of NOAA-AVHRR-NDVI bas
ed on time series analysis (1981-1992)”, In
t. J. Remote Sens., 17231-225, 1996.
12) Ricotta, C., G. Avena and A. De Palma,
“Mapping and monitoring net primary produc
tivity with AVHRR NDVI time-series: statistic
al equivalence of cumulative vegetation indic
es”, ISPRS J.Photogramm. Remote Sens., 54,
325-331, 1999.
13) Azzali, S. and M. Menenti, “Mapping ve
getation-soil-climatecomplexes in southern Afri
ca using temporal Fourier analysis of NOAA-A
VHRR NDVI data”, , Int. J. Remote Sens.,
21, 973-996,2000.
14) Roerink, G. J., M. Menenti and W. Verh
oef, “Reconstructingcloudfree NDVI composites
using Fourier analysis of time series”, I
nt. J. Remote Sens., 21, 1911-1917, 2000.
15) 澤田義人、沢田治雄、斎藤英樹、“高頻度観測衛
星による植生の季節変動観測”、日本写真測量学会講演
論文集、2000秋季、209-212、2000
16) 北川源四郎、“FORTRAN77 時系列解析プログラミ
ング”、岩波書店、第5章、1993
17) http://www.affrc.go.jp/Agropedia/
18) 児玉正文、宋 献方、“農林水産リモートセンシ
ングデータベースシステムの構築”、日本リモートセン
シング学会第28回学術講演会論文集、259-260、2000
19) Pathfinder web site:http://daac.gsfc.nasa.g
ov/CAMPAIGN_DOCS/LAND_BIO/GLBDST_main.html
20) http://edcdaac.usgs.gov/1KM/1kmhomepage.htm
l
21) Agbu, P. A. and M. E. James, “NOAA/NA
SA Pathfinder AVHRRLand Data Set User's Manu
al. Goddard Distributed Active Archive Cente
r”, NASA Goddard Space Flight Center, Green
belt, 1994.http://daac.gsfc.nasa.gov/CAMPAIGN_DOC
S/LAND_BIO/AVHRR_GD.pdf
22) 北川源四郎、“FORTRAN77 時系列解析プログラミ
ング”、岩波書店、第9章、1993
23) 西山清、“周波数が既知である成分のみを含む時
変スペクトルの一推定法”、電子情報通信学会論文誌、
J74-A、6、916-918、1991
24) Martin, M. P., P. Ceccato, S. Flasse a
nd I. Downey, “Fire detection and fire gro
wth monitoring using satellite data”,Remote
Sensing of Large Wildfires in the European
MediterraneanBasin, E. Chuvieco (ed.), Splinge
r-Verlag, chapt. 6, 1999.
25) Prins, E. M. and W. P. Menzel, “Geost
ationary satellite detection of biomass burni
ng is South America”, Int. J. RemoteSens.,
13, 2783-2799, 1992.1) Indonesia Forest Fire Prevention Project, “Indonesia National Park Forest Fire Survey Report”, Japan International Cooperation Agency, 1998 2) Harada Sawada, “Use of Remote Sensing for Large-scale Forest Fires”, Forest Science, 24 , 22-28, 1998 3) Yoshiaki Honda, Shunji Murai, Kikuo Kato, "World Vegetation Monitoring", Journal of Photogrammetric Society of Japan "Photogrammetry and Remote Sensing", 31, 1, 4-14, 1992 4) Defries, RS and JRG Townshed,
“NDVI-derived land cover classifications at
a global scale ”, Int. J. Remote Sens., 15,
3567-3586, 1994.5) Myneni, RB, CJ Tucker, G. Asrar
and CD Keeling, “Interannual variations in
satellite-sensed vegetation index data from
1981 to 1991 ”, J. Geophys. Res., 103, D
6, 6145-6160, 1998. 6) Los, SO, CO Justice and CJ T
ucker, “A global 1deg. by 1 deg. NDVI dat
a set for climate studies derived from the
GIMMS continental NDVI Data. ”, Int. J. Re
mote Sens., 15, 3493-3518, 1994.7) GLCC web site: http://edcdaac.usgs.gov/gl
cc / globdoc1_2.html 8) Viovy, N., O. Arino and AS Belward,
“The Best Index Slope Extraction (BISE):
A method for reducing noise in NDVI time-se
ries ”,, Int. J. Remote Sens., 13, 1585-15
90, 1992. 9) Lee Unqing, Kazuhiko Onuma, Yoshizumi Yasuda, “Fourier analysis of time series vegetation index data”, Journal of Photographic Survey of Japan, “Photogrammetry and Remote Sensing”, 31, 4, 4-14, 1992 10 ) Sellers, PJ, CJ Tucker, GJ C
ollatz, SO Los, CO Justice, DA Daz
lich and DA Randall, “A global1 deg. b
y 1 deg.NDVI data set for climate studie
s. Part 2: The generation of global fields
of terrestrial biophysical parameters from
the NDVI ”, Int. J. Remote Sens., 15, 3519-
3545, 1994. 11) Verhoef, W., M. Menenti and S. Azzali,
"A color composite of NOAA-AVHRR-NDVI bas
ed on time series analysis (1981-1992) ”, In
t. J. Remote Sens., 17231-225, 1996.12) Ricotta, C., G. Avena and A. De Palma,
“Mapping and monitoring net primary produc
tivity with AVHRR NDVI time-series: statistic
al equivalence of cumulative vegetation indic
es ”, ISPRS J. Photogramm. Remote Sens., 54,
325-331, 1999.13) Azzali, S. and M. Menenti, “Mapping ve
getation-soil-climatecomplexes in southern Afri
ca using temporal Fourier analysis of NOAA-A
VHRR NDVI data ”,, Int. J. Remote Sens.,
21, 973-996, 2000. 14) Roerink, GJ, M. Menenti and W. Verh
oef, “Reconstructingcloudfree NDVI composites
using Fourier analysis of time series ”, I
nt. J. Remote Sens., 21, 1911-1917, 2000. 15) Yoshito Sawada, Haruo Sawada, Hideki Saito, “Observation of seasonal variation of vegetation by high-frequency observation satellites”, Proc. of Japan Photogrammetric Society, Autumn 2000 , 209-212, 2000 16) Genshiro Kitagawa, “FORTRAN77 Time Series Analysis Programming”, Iwanami Shoten, Chapter 5, 1993 17) http://www.affrc.go.jp/Agropedia/ 18) Masafumi Kodama, Song Song , "Agriculture, Forestry and Fisheries Remote Sensing Database System", Proc. Of the 28th conference of the Remote Sensing Society of Japan, 259-260, 2000 19) Pathfinder web site: http: //daac.gsfc.nasa.g
ov / CAMPAIGN_DOCS / LAND_BIO / GLBDST_main.html 20) http://edcdaac.usgs.gov/1KM/1kmhomepage.htm
l 21) Agbu, PA and ME James, “NOAA / NA
SA Pathfinder AVHRRLand Data Set User's Manu
al. Goddard Distributed Active Archive Cente
r ”, NASA Goddard Space Flight Center, Green
belt, 1994.http: //daac.gsfc.nasa.gov/CAMPAIGN_DOC
S / LAND_BIO / AVHRR_GD.pdf 22) Genshiro Kitagawa, “FORTRAN77 Time Series Analysis Programming,” Iwanami Shoten, Chapter 9, 1993 23) Kiyoshi Nishiyama, “One estimation method for time-varying spectra containing only known frequency components” ", The Institute of Electronics, Information and Communication Engineers,
J74-A, 6, 916-918, 1991 24) Martin, MP, P. Ceccato, S. Flasse a
nd I. Downey, “Fire detection and fire gro
wth monitoring using satellite data ”, Remote
Sensing of Large Wildfires in the European
MediterraneanBasin, E. Chuvieco (ed.), Splinge
r-Verlag, chapt. 6, 1999. 25) Prins, EM and WP Menzel, “Geost
ationary satellite detection of biomass burni
ng is South America ”, Int. J. RemoteSens.,
13, 2783-2799, 1992.
【0091】[0091]
【発明の効果】本発明によると、高頻度観測衛星からの
データの画素毎の時系列的な特徴を保存し、時間分解能
を完全に保持したまま、雲などの影響およびシステムノ
イズの除去をすることができる。また、本発明による
と、時間分解能を保持したまま植生などの変動の様子が
捉えられるため、防災、環境、農業分野での広域現状把
握、および予測が容易となる。さらに、本発明による
と、適切な関数を自動的に選択し、ノイズ除去処理を自
動で行えるものであり、定常業務の一部分とすることが
可能である。さらに本発明によると、水分指数など他の
地球高頻度観測衛星のデータにもそのまま転用および応
用ができる。According to the present invention, the time-series characteristics of each pixel of the data from the high-frequency observation satellite are preserved, and the influence of clouds and the system noise are removed while completely maintaining the time resolution. be able to. Further, according to the present invention, it is possible to grasp the state of changes in vegetation while maintaining the time resolution, and thus it is easy to understand and predict the wide area current situation in the disaster prevention, environment, and agricultural fields. Further, according to the present invention, it is possible to automatically select an appropriate function and automatically perform the noise removal processing, and it is possible to form a part of the routine work. Furthermore, according to the present invention, the data can be diverted and applied as it is to the data of other high frequency earth observation satellites such as the moisture index.
【図1】シミュレーションによるフィルタ特性の評価図FIG. 1 Evaluation diagram of filter characteristics by simulation
【図2】シミュレーション結果の例を示す図FIG. 2 is a diagram showing an example of a simulation result.
【図3】命令列スケジュール処理を実行するための装置
の構成図FIG. 3 is a block diagram of an apparatus for executing an instruction sequence schedule process.
【図4】データ処理のフローチャートFIG. 4 is a flowchart of data processing.
【図5】LMF処理前後のNDVIデータ図FIG. 5: NDVI data diagram before and after LMF processing
【図6】LMF処理前後のデータの相関図FIG. 6 is a correlation diagram of data before and after LMF processing.
【図7】AVHRR/Ch4のデータをLMF処理した
データ図FIG. 7 is a data diagram of LHR processing of AVHRR / Ch4 data.
【図8】PathfinderデータのNDVIをLM
F−KF処理した結果を示す図[FIG. 8] LMVI of NDVI of Pathfinder data
The figure which shows the result of F-KF processing.
【図9】AVHRR/Ch4輝度温度データをLMF−
KF処理した結果を示す図FIG. 9: AVHRR / Ch4 brightness temperature data is LMF-
Diagram showing the result of KF processing
1 処理部 2 入力部 3 出力部 4 記憶部 41 計算パラメータファイル 42 第1の観測データファイル 43 第2の観測データファイル 44 第1の推定観測データファイル 45 第2の推定観測データファイル 46 第1のパラメータファイル 47 第2のパラメータファイル 1 processing unit 2 Input section 3 Output section 4 memory 41 Calculation parameter file 42 First observation data file 43 Second observation data file 44 First estimated observation data file 45 Second estimated observation data file 46 First parameter file 47 Second parameter file
───────────────────────────────────────────────────── フロントページの続き (72)発明者 澤田 義人 茨城県稲敷郡茎崎町高見原1丁目1番179 号 ランドスケープB206 (72)発明者 沢田 治雄 茨城県つくば市春日4丁目13番19号 Fターム(参考) 5B057 AA14 CA08 CA12 CA16 CB08 CB12 CB16 CE02 DC32 ─────────────────────────────────────────────────── ─── Continued front page (72) Inventor Yoshito Sawada 1-1179 Takamihara, Kukizaki-cho, Inashiki-gun, Ibaraki Prefecture No. Landscape B206 (72) Inventor Haruo Sawada 4-13-19 Kasuga, Tsukuba City, Ibaraki Prefecture F term (reference) 5B057 AA14 CA08 CA12 CA16 CB08 CB12 CB16 CE02 DC32
Claims (11)
応して、衛星による植生の季節変動についての時系列の
第1の観測データを記憶した第1の観測データファイル
から、画素毎の第1の観測データを読み込むステップ
と、 処理部は、各画素毎に、ある時期又は時点に対応して、
該時期又は時点の前後所定期間の局所最大値を抽出する
ことにより時系列の第2の観測データを求め、それを各
画素に対応して第2の観測データファイルに記憶するス
テップと、 処理部は、第2の観測データファイルに記憶された第2
の観測データを、長期的非周期項及び季節変化項を含む
次式の関数を用い、各第1の係数ciを最小二乗法又は
他の推定法により決定することでフィッティングし、各
第1の係数ciを第1のパラメータファイルに記憶する
ステップと、 【数1】 [ ここで、 ci:求める第1の係数 t:時期又は時間 N:周期関数のペア数(この場合sinとcosのペア
数) l:周期関数の仮の番号 kl:12ヶ月が回帰の整数倍になる数列(例:1,
2,3,4,6,12) M:1年間のデータ数(シーン数) ] 処理部は、第1の係数ciで定められた前記関数に従
い、第2の観測データから雑音データを除去した第1の
推定観測データを求め、それを第1の推定観測データフ
ァイルに記憶する、又は、出力部に出力若しくは表示す
るステップとを含む地球観測衛星データのノイズ除去方
法。1. A processing unit for each pixel from a first observation data file that stores time-series first observation data regarding seasonal variation of vegetation by satellites corresponding to each pixel in an observation target area. The step of reading the first observation data of, the processing unit, for each pixel, corresponding to a certain time or point,
A step of obtaining time-series second observation data by extracting a local maximum value in a predetermined period before and after the time or point of time, and storing the second observation data in a second observation data file corresponding to each pixel; Is the second stored in the second observation data file.
The observed data of 1 is fitted by determining each first coefficient c i by the method of least squares or another estimation method using the function of the following equation including the long-term aperiodic term and the seasonal variation term, and Storing the coefficients c i of the first parameter file in the first parameter file; [Where, c i : first coefficient to be obtained t: time or time N: number of pairs of periodic function (in this case, number of pairs of sin and cos) l: temporary number of periodic function k l : 12 months of regression Sequences that are integer multiples (eg, 1,
2, 3, 4, 6, 12) M: number of data in one year (number of scenes)] The processing unit removes noise data from the second observation data according to the function defined by the first coefficient c i. Determining the first estimated observation data, storing the first estimated observation data in the first estimated observation data file, or outputting or displaying the first estimated observation data file in the output unit.
の時系列の状態空間モデルを適用し、該状態空間モデル
の各第2の係数ci(t)を推定する処理を実行し、各
第2の係数ci(t)を第2のパラメータファイルに記
憶するステップと、 【数2】 [ ここで、 ci(t):求める第2の係数 wk(t):システム雑音 y(t):観測される値 v(t):平均0、分散がσvの白色雑音(観測雑音) 【数3】 ]処理部は、第2の係数ci(t)で定められた前記関
数に従い、第2の観測データから雑音データを除去した
第2の推定観測データを求め、それを第2の推定観測デ
ータファイルに記憶する、又は出力部に出力、若しくは
表示するステップをさらに含む請求項1に記載の地球観
測衛星データのノイズ除去方法。2. The processing unit applies a time-series state space model of the following equation to the second observation data, and executes a process of estimating each second coefficient c i (t) of the state space model. , Storing each second coefficient c i (t) in a second parameter file; [Where, c i (t): second coefficient to be obtained w k (t): system noise y (t): observed value v (t): average 0, white noise with variance σ v (observation noise ) [Equation 3] ] The processing unit obtains the second estimated observation data by removing the noise data from the second observation data according to the function defined by the second coefficient c i (t), and calculates the second estimated observation data. The method for removing noise from earth observation satellite data according to claim 1, further comprising a step of storing in a file, outputting to an output unit, or displaying.
までのデータについて最大値を抽出するステップと、時
刻kに対応して、時期kから時期k+nのデータについ
て最大値を抽出するステップと、時刻kに対応して、得
られたこれら二つの最大値のうち、値の小さいものを局
所的最大値とするステップと、を含む請求項1又は2に
記載の地球観測衛星データのノイズ除去方法。3. Corresponding to time k, time k-n to time k
To extract the maximum value for the data up to, the step to extract the maximum value for the data from time k to time k + n corresponding to time k, and these two maximum values obtained corresponding to time k The method of removing noise from earth observation satellite data according to claim 1 or 2, further comprising the step of: making the one having a smaller value a local maximum value.
が、設定したしきい値よりも大きい場合には雲やシステ
ムノイズの影響とみなしてその時期の観測データは採用
しないようにして、採用されたデータだけを用いて再
度、係数及び/又は推定観測データを求めることを特徴
とする請求項1又は2に記載の地球観測衛星データのノ
イズ除去方法。4. If the difference between the estimated value and the actual pixel value is larger than the set threshold value at the time k, it is regarded as the influence of clouds or system noise and the observation data at that time is not adopted. The method for removing noise from earth observation satellite data according to claim 1 or 2, wherein the coefficient and / or the estimated observation data is obtained again using only the adopted data.
タ、水分指数、中間赤外線バンド、熱赤外線輝度バン
ド、温度データのいずれか又は複数であることを特徴と
する請求項1又は2に記載の地球観測衛星データのノイ
ズ除去方法。5. The first observation data is any one or a plurality of vegetation index data, moisture index, mid-infrared band, thermal infrared brightness band, and temperature data. Method for removing earth observation satellite data from Japan.
の大気状態に起因するもの、及び/又は、センサ感度等
のシステムに起因するものを含むことを特徴とする請求
項1又は2に記載の地球観測衛星データのノイズ除去方
法。6. The noise data includes data derived from atmospheric conditions such as clouds, fog, rain and haze, and / or data derived from a system such as sensor sensitivity. The method for removing noise from earth observation satellite data according to 2.
複数を組み合わせて出力部に出力又は表示するステップ
をさらに含む請求項1又は2に記載の地球観測衛星デー
タのノイズ除去方法。7. The method for removing noise from earth observation satellite data according to claim 1, wherein the processing unit further comprises a step of combining a plurality of observation data and estimation data and outputting or displaying the combination on an output unit.
応して、衛星による植生の季節変動についての時系列の
第1の観測データを記憶した第1の観測データファイル
から、画素毎の第1の観測データを読み込むステップ
と、 処理部は、各画素毎に、ある時期又は時点に対応して、
該時期又は時点の前後所定期間の局所最大値を抽出する
ことにより時系列の第2の観測データを求め、それを各
画素に対応して第2の観測データファイルに記憶するス
テップと、 処理部は、第2の観測データファイルに記憶された第2
の観測データを、長期的非周期項及び季節変化項を含む
次式の関数を用い、各第1の係数ciを最小二乗法又は
他の推定法により決定することでフィッティングし、各
第1の係数ciを第1のパラメータファイルに記憶する
ステップと、 【数4】 [ ここで、 ci:求める第1の係数 t:時期又は時間 N:周期関数のペア数(この場合sinとcosのペア
数): l:周期関数の仮の番号 kl:12ヶ月が回帰の整数倍になる数列(例:1,
2,3,4,6,12) M:1年間のデータ数(シーン数) ] 処理部は、第1の係数ciで定められた前記関数に従
い、第2の観測データから雑音データを除去した第1の
推定観測データを求め、それを第1の推定観測データフ
ァイルに記憶する、又は出力部に出力若しくは表示する
ステップをコンピュータに実行させるための地球観測衛
星データのノイズ除去処理プログラム。8. The processing unit stores, for each pixel, from a first observation data file in which time-series first observation data regarding seasonal variation of vegetation by a satellite is stored, corresponding to each pixel in the observation target area. The step of reading the first observation data of, the processing unit, for each pixel, corresponding to a certain time or point,
A step of obtaining time-series second observation data by extracting a local maximum value in a predetermined period before and after the time or point of time, and storing the second observation data in a second observation data file corresponding to each pixel; Is the second stored in the second observation data file.
The observed data of 1 is fitted by determining each first coefficient c i by the method of least squares or another estimation method using the function of the following equation including the long-term aperiodic term and the seasonal variation term, and Storing the coefficients c i of the in the first parameter file; [Where, c i : first coefficient to be obtained t: time or time N: number of pairs of periodic function (in this case, number of pairs of sin and cos): l: temporary number of periodic function k l : 12 months regression A sequence that is an integer multiple of (example: 1,
2, 3, 4, 6, 12) M: number of data in one year (number of scenes)] The processing unit removes noise data from the second observation data according to the function defined by the first coefficient c i. A noise removal processing program for earth observation satellite data for causing a computer to perform the step of obtaining the first estimated observation data, storing the first estimated observation data in the first estimated observation data file, or outputting or displaying the first estimated observation data file.
の時系列の状態空間モデルを適用し、該状態空間モデル
の各第2の係数ci(t)を推定する処理を実行し、各
第2の係数ci(t)を第2のパラメータファイルに記
憶するステップと、 【数5】 [ ここで、 ci(t):求める第2の係数 wk(t):システム雑音 y(t):観測される値 v(t):平均0、分散がσvの白色雑音(観測雑音) 【数6】 ] 処理部は、第2の係数ci(t)で定められた前記関数
に従い、第2の観測データから雑音データを除去した第
2の推定観測データを求め、それを第2の推定観測デー
タファイルに記憶する、又は出力部に出力若しくは表示
するステップをさらにコンピュータに実行させるための
請求項13に記載の地球観測衛星データのノイズ除去処
理プログラム。9. A processing unit applies a time-series state space model of the following equation to the second observation data and executes a process of estimating each second coefficient c i (t) of the state space model. , Storing each second coefficient c i (t) in a second parameter file; [Where, c i (t): second coefficient to be obtained w k (t): system noise y (t): observed value v (t): average 0, white noise with variance σ v (observation noise ) [Equation 6] ] The processing unit obtains the second estimated observation data by removing the noise data from the second observation data according to the function defined by the second coefficient c i (t), and obtains the second estimated observation data. The noise removal processing program for earth observation satellite data according to claim 13, which causes a computer to further execute the step of storing in a file or outputting or displaying on an output unit.
対応して、衛星による植生の季節変動についての時系列
の第1の観測データを記憶した第1の観測データファイ
ルから、画素毎の第1の観測データを読み込むステップ
と、 処理部は、各画素毎に、ある時期又は時点に対応して、
該時期又は時点の前後所定期間の局所最大値を抽出する
ことにより時系列の第2の観測データを求め、それを各
画素に対応して第2の観測データファイルに記憶するス
テップと、 処理部は、第2の観測データファイルに記憶された第2
の観測データを、長期的非周期項及び季節変化項を含む
次式の関数を用い、各第1の係数ciを最小二乗法又は
他の推定法により決定することでフィッティングし、各
第1の係数ciを第1のパラメータファイルに記憶する
ステップと、 【数7】 [ ここで、 ci:求める第1の係数 t:時期又は時間 N:周期関数のペア数(この場合sinとcosのペア
数): l:周期関数の仮の番号 kl:12ヶ月が回帰の整数倍になる数列(例:1,
2,3,4,6,12) M:1年間のデータ数(シーン数) ] 処理部は、第1の係数ciで定められた前記関数に従
い、第2の観測データから雑音データを除去した第1の
推定観測データを求め、それを第1の推定観測データフ
ァイルに記憶する、又は出力部に出力、若しくは表示す
るステップをコンピュータに実行させるための地球観測
衛星データのノイズ除去処理プログラムを記録したコン
ピュータ読み取り可能な記録媒体。10. The processing unit stores, for each pixel, from a first observation data file in which time-series first observation data regarding seasonal variation of vegetation by a satellite is stored, corresponding to each pixel in the observation target area. The step of reading the first observation data of, the processing unit, for each pixel, corresponding to a certain time or point,
A step of obtaining time-series second observation data by extracting a local maximum value in a predetermined period before and after the time or point of time, and storing the second observation data in a second observation data file corresponding to each pixel; Is the second stored in the second observation data file.
The observed data of 1 is fitted by determining each first coefficient c i by the method of least squares or another estimation method using the function of the following equation including the long-term aperiodic term and the seasonal variation term, and Storing the coefficient c i of the first parameter file in the first parameter file; [Where, c i : first coefficient to be obtained t: time or time N: number of pairs of periodic function (in this case, number of pairs of sin and cos): l: temporary number of periodic function k l : 12 months regression A sequence that is an integer multiple of (example: 1,
2, 3, 4, 6, 12) M: number of data in one year (number of scenes)] The processing unit removes noise data from the second observation data according to the function defined by the first coefficient c i. A noise removal processing program for earth observation satellite data for causing the computer to execute the step of obtaining the first estimated observation data, storing it in the first estimated observation data file, or outputting or displaying it on the output unit. The recorded computer-readable recording medium.
式の時系列の状態空間モデルを適用し、該状態空間モデ
ルの各第2の係数ci(t)を推定する処理を実行し、
各第2の係数ci(t)を第2のパラメータファイルに
記憶するステップと、 【数8】 [ ここで、 ci(t):求める第2の係数 Wk(t):システム雑音 y(t):観測される値 v(t):平均0、分散がσvの白色雑音(観測雑音) 【数9】 ] 処理部は、第2の係数ci(t)で定められた前記関数
に従い、第2の観測データから雑音データを除去した第
2の推定観測データを求め、それを第2の推定観測デー
タファイルに記憶する、又は出力部に出力、若しくは表
示するステップをさらにコンピュータに実行させるため
の地球観測衛星データのノイズ除去処理プログラムを記
録した請求項15に記載のコンピュータ読み取り可能な
記録媒体。11. A processing unit applies a time-series state space model of the following equation to the second observation data and executes a process of estimating each second coefficient c i (t) of the state space model. ,
Storing each second coefficient c i (t) in a second parameter file; [Where, c i (t): second coefficient Wk (t) to be obtained: system noise y (t): observed value v (t): white noise with mean 0 and variance σ v (observation noise) [Equation 9] ] The processing unit obtains the second estimated observation data by removing the noise data from the second observation data according to the function defined by the second coefficient c i (t), and obtains the second estimated observation data. 16. The computer-readable recording medium according to claim 15, which stores a noise removal processing program for earth observation satellite data for causing a computer to further execute the step of storing in a file or outputting or displaying on an output unit.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002103966A JP4003869B2 (en) | 2002-04-05 | 2002-04-05 | Noise removal processing method for earth observation satellite data, noise removal processing program, and recording medium recording noise removal processing program |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002103966A JP4003869B2 (en) | 2002-04-05 | 2002-04-05 | Noise removal processing method for earth observation satellite data, noise removal processing program, and recording medium recording noise removal processing program |
Publications (3)
Publication Number | Publication Date |
---|---|
JP2003296702A true JP2003296702A (en) | 2003-10-17 |
JP2003296702A5 JP2003296702A5 (en) | 2005-06-16 |
JP4003869B2 JP4003869B2 (en) | 2007-11-07 |
Family
ID=29389479
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2002103966A Expired - Fee Related JP4003869B2 (en) | 2002-04-05 | 2002-04-05 | Noise removal processing method for earth observation satellite data, noise removal processing program, and recording medium recording noise removal processing program |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP4003869B2 (en) |
Cited By (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2008281505A (en) * | 2007-05-14 | 2008-11-20 | Japan Aerospace Exploration Agency | Method of resampling 3-dimensional discrete data, and its device |
JP2011242288A (en) * | 2010-05-19 | 2011-12-01 | Honda Elesys Co Ltd | Electronic scanning radar device, reception wave direction estimation method, and reception wave direction estimation program |
KR101089220B1 (en) | 2010-02-18 | 2011-12-02 | 공주대학교 산학협력단 | Normalized difference vegetation index correction method using spatio-temporal continuity |
JP2012013569A (en) * | 2010-07-01 | 2012-01-19 | Honda Elesys Co Ltd | Electronic scanning type radar apparatus, reception wave direction estimating method, and reception wave direction estimation program |
KR101404430B1 (en) * | 2013-06-11 | 2014-06-10 | 서울시립대학교 산학협력단 | Method for estimation of surface temperature lapse rate Using thermal infrared images |
CN104251662A (en) * | 2013-06-27 | 2014-12-31 | 杭州中科天维科技有限公司 | Ordered point cloud threshold adaptive noise suppression technology |
CN104502241A (en) * | 2014-05-14 | 2015-04-08 | 淮海工学院 | GNSS based haze monitoring and analysis technology |
CN107888878A (en) * | 2017-11-16 | 2018-04-06 | 包云清 | Fire automatic-positioning type photograph platform courses method |
WO2019160161A1 (en) | 2018-02-19 | 2019-08-22 | 株式会社Ihi | Heat source detection device |
JP2020085494A (en) * | 2018-11-16 | 2020-06-04 | ヤフー株式会社 | Information processor, information processing method, and information processing program |
CN111553237A (en) * | 2020-04-23 | 2020-08-18 | 江西理工大学 | LJ1-01 night light data denoising method based on multi-state superposition Gamma distribution |
CN111650102A (en) * | 2020-05-26 | 2020-09-11 | 北京中科锐景科技有限公司 | Haze pollution analysis method, device, medium and equipment based on satellite data |
JP2021533427A (en) * | 2019-07-04 | 2021-12-02 | 浙江大学Zhejiang University | Reconstruction method of geostationary sea color satellite data by empirical orthogonal function decomposition method |
WO2022080418A1 (en) * | 2020-10-14 | 2022-04-21 | 株式会社ブリヂストン | Leaf abscission detection system and leaf abscission detection method |
KR102502154B1 (en) * | 2022-07-08 | 2023-02-22 | 에스케이임업 주식회사 | CO2 Capture Prediction System |
CN117668468A (en) * | 2024-01-31 | 2024-03-08 | 山东省舜天化工集团有限公司 | Intelligent analysis management system for chemical preparation data |
-
2002
- 2002-04-05 JP JP2002103966A patent/JP4003869B2/en not_active Expired - Fee Related
Cited By (30)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2008281505A (en) * | 2007-05-14 | 2008-11-20 | Japan Aerospace Exploration Agency | Method of resampling 3-dimensional discrete data, and its device |
WO2008143077A1 (en) * | 2007-05-14 | 2008-11-27 | Japan Aerospace Exploration Agency | 3-dimensional discrete data re-sampling method and device |
KR101089220B1 (en) | 2010-02-18 | 2011-12-02 | 공주대학교 산학협력단 | Normalized difference vegetation index correction method using spatio-temporal continuity |
JP2011242288A (en) * | 2010-05-19 | 2011-12-01 | Honda Elesys Co Ltd | Electronic scanning radar device, reception wave direction estimation method, and reception wave direction estimation program |
JP2012013569A (en) * | 2010-07-01 | 2012-01-19 | Honda Elesys Co Ltd | Electronic scanning type radar apparatus, reception wave direction estimating method, and reception wave direction estimation program |
US10337925B2 (en) | 2013-06-11 | 2019-07-02 | University of Seoul Cooperation Foundation | Method for estimating land surface temperature lapse rate using infrared image |
KR101404430B1 (en) * | 2013-06-11 | 2014-06-10 | 서울시립대학교 산학협력단 | Method for estimation of surface temperature lapse rate Using thermal infrared images |
WO2014200258A1 (en) * | 2013-06-11 | 2014-12-18 | 서울시립대학교 산학협력단 | Method for estimating surface lapse rate using infrared image |
CN104251662A (en) * | 2013-06-27 | 2014-12-31 | 杭州中科天维科技有限公司 | Ordered point cloud threshold adaptive noise suppression technology |
CN104502241A (en) * | 2014-05-14 | 2015-04-08 | 淮海工学院 | GNSS based haze monitoring and analysis technology |
CN107888878A (en) * | 2017-11-16 | 2018-04-06 | 包云清 | Fire automatic-positioning type photograph platform courses method |
JP6996573B2 (en) | 2018-02-19 | 2022-01-17 | 株式会社Ihi | Heat source detector |
WO2019160161A1 (en) | 2018-02-19 | 2019-08-22 | 株式会社Ihi | Heat source detection device |
JPWO2019160161A1 (en) * | 2018-02-19 | 2020-12-10 | 株式会社Ihi | Heat source detector |
AU2019221160B2 (en) * | 2018-02-19 | 2021-12-16 | Ihi Corporation | Heat source detection device |
US11933673B2 (en) | 2018-02-19 | 2024-03-19 | Ihi Corporation | Heat source detection device |
JP2020085494A (en) * | 2018-11-16 | 2020-06-04 | ヤフー株式会社 | Information processor, information processing method, and information processing program |
JP7036704B2 (en) | 2018-11-16 | 2022-03-15 | ヤフー株式会社 | Information processing equipment, information processing methods, and information processing programs |
JP2021533427A (en) * | 2019-07-04 | 2021-12-02 | 浙江大学Zhejiang University | Reconstruction method of geostationary sea color satellite data by empirical orthogonal function decomposition method |
US11790580B2 (en) | 2019-07-04 | 2023-10-17 | Zhejiang University | Method for reconstructing geostationary ocean color satellite data based on data interpolating empirical orthogonal functions |
JP7004844B2 (en) | 2019-07-04 | 2022-01-21 | 浙江大学 | Reconstruction method of geostationary sea color satellite data by empirical orthogonal function decomposition method |
CN111553237A (en) * | 2020-04-23 | 2020-08-18 | 江西理工大学 | LJ1-01 night light data denoising method based on multi-state superposition Gamma distribution |
CN111553237B (en) * | 2020-04-23 | 2023-11-07 | 江西理工大学 | LJ1-01 night lamplight data denoising method based on polymorphic superposition Gamma distribution |
CN111650102A (en) * | 2020-05-26 | 2020-09-11 | 北京中科锐景科技有限公司 | Haze pollution analysis method, device, medium and equipment based on satellite data |
CN111650102B (en) * | 2020-05-26 | 2023-08-29 | 北京中科锐景科技有限公司 | Haze pollution analysis method, device, medium and equipment based on satellite data |
WO2022080418A1 (en) * | 2020-10-14 | 2022-04-21 | 株式会社ブリヂストン | Leaf abscission detection system and leaf abscission detection method |
KR102502154B1 (en) * | 2022-07-08 | 2023-02-22 | 에스케이임업 주식회사 | CO2 Capture Prediction System |
WO2024010127A1 (en) * | 2022-07-08 | 2024-01-11 | 에스케이임업 주식회사 | System for predicting amount of captured carbon dioxide |
CN117668468A (en) * | 2024-01-31 | 2024-03-08 | 山东省舜天化工集团有限公司 | Intelligent analysis management system for chemical preparation data |
CN117668468B (en) * | 2024-01-31 | 2024-04-26 | 山东省舜天化工集团有限公司 | Intelligent analysis management system for chemical preparation data |
Also Published As
Publication number | Publication date |
---|---|
JP4003869B2 (en) | 2007-11-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Dixon et al. | Satellite prediction of forest flowering phenology | |
McEvoy et al. | The evaporative demand drought index. Part II: CONUS-wide assessment against common drought indicators | |
Sobrino et al. | Toward remote sensing methods for land cover dynamic monitoring: Application to Morocco | |
JP2003296702A (en) | Method for denoising earth observation satellite data, denoising processing program and recording medium recorded with denoising processing program | |
Knipper et al. | Using high-spatiotemporal thermal satellite ET retrievals to monitor water use over California vineyards of different climate, vine variety and trellis design | |
Fernandez-Carrillo et al. | Estimating prescribed fire impacts and post-fire tree survival in eucalyptus forests of Western Australia with L-band SAR data | |
Katagis et al. | Trend analysis of medium-and coarse-resolution time series image data for burned area mapping in a Mediterranean ecosystem | |
De Kauwe et al. | Quantifying land surface temperature variability for two Sahelian mesoscale regions during the wet season | |
Zhang et al. | Understory biomass measurement in a dense plantation forest based on drone-SfM data by a manual low-flying drone under the canopy | |
Lee et al. | A comparative study on generating simulated Landsat NDVI images using data fusion and regression method—the case of the Korean Peninsula | |
Liang et al. | Utilizing digital image processing and two-source energy balance model for the estimation of evapotranspiration of dry edible beans in western Nebraska | |
Santos et al. | Predicting eucalyptus plantation growth and yield using Landsat imagery in Minas Gerais, Brazil | |
Zang et al. | Spatially-explicit mapping annual oil palm heights in peninsular Malaysia combining ICESat-2 and stand age data | |
Kim et al. | Multi-temporal nonlinear regression method for landsat image simulation | |
Clemente et al. | Monitoring post-fire regeneration in Mediterranean ecosystems by employing multitemporal satellite imagery | |
Anurogo et al. | An integrated comparative approach to estimating forest aboveground carbon stock using advanced remote sensing technologies | |
Reyes-Gonzalez | Using remote sensing to estimate crop water use to improve irrigation water management | |
Bhattacharya et al. | A baseline regional evapotranspiration (ET) and change hotspots over Indian sub-tropics using satellite remote sensing data | |
Liu et al. | An improved combined vegetation difference index and burn scar index approach for mapping cropland burned areas using combined data from Landsat 8 multispectral and thermal infrared bands | |
Dineshkumar et al. | Rice monitoring using sentinel-1 data in the *** earth engine platform | |
Praetyo et al. | Mitigation & identification for local aridity, based of vegetation indices combined with spatial statistics & clustering k means | |
Bezerra et al. | STEEP: A remotely-sensed energy balance model for evapotranspiration estimation in seasonally dry tropical forests | |
Verma et al. | Assessment of aboveground biomass in a tropical dry deciduous forest using PRISMA data | |
Santarsiero et al. | Fire severity and vegetation recovery determination using GEE and Sentinel-2: the case of Peschici fire | |
Dubois et al. | Copernicus sentinel-2 data for the determination of groundwater withdrawal in the maghreb region |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A712 Effective date: 20031031 |
|
RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20040129 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040922 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040922 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20070711 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20070724 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20070816 |
|
R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100831 Year of fee payment: 3 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110831 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110831 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120831 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130831 Year of fee payment: 6 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
LAPS | Cancellation because of no payment of annual fees |