JP2003296702A - Method for denoising earth observation satellite data, denoising processing program and recording medium recorded with denoising processing program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To remove the influence of clouds and haze from data of a high frequency observation satellite while maintaining spatial resolution and while realizing optional time resolution. <P>SOLUTION: A filter for extracting a local maximum value is acted about time series data of a processing object pixel regarding around n time (S101 to S105). A function is assigned to a local maximum value group obtained there to estimate a pixel value (S107 to S123). When the deviation between the estimated value and the actual pixel value is larger than a set threshold, it is considered that the deviation results from the influence of cloud and system noise to prevent the data at that time from being adopted. Only adopted data is used to calculate a local maximum value and to assign a function again. This processing is repeated if necessary until a result is settled. The above processing is performed for the entire image, the influence of the cloud and the system noise can be eliminated. <P>COPYRIGHT: (C)2004,JPO

Description

Detailed Description of the Invention

[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.

[Prior art]

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.

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.

[Problems to be Solved by the Invention]

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.

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.

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.

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] 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.

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.

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.

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]

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.

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.

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.

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

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,

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]

[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]

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]

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.

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]

[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 ₀ is the average, c ₁ 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.

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.

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]

[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).

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.

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.

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]

[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.

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).

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.

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).

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.

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).

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.

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]

[Equation 14] Here, each parameter is the same as equation (1), but c
_There is no ₁ (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]

[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]

[Equation 16] Is a random walk scatter model

[0040]

[Equation 17] And a model assuming that the change rates at time t and time t−1 are equal,

[0041]

[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]

[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]

[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]

[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]

[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.

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.

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]

[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]

[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]

[Equation 25] (However, ε> 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]

[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 ⁻⁴ , 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.

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]

[Equation 27] (However, ε> 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]

[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.

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.

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.

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.

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).

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.

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.

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).

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).

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).

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).

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.

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.)

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.

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.

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.

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.

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.

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. .

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. .

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.

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.

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.

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.

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.

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. .

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.

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.

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:

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.

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:

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.

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:

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.

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.

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.

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]

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.

[Brief description of drawings]

FIG. 1 Evaluation diagram of filter characteristics by simulation

FIG. 2 is a diagram showing an example of a simulation result.

FIG. 3 is a block diagram of an apparatus for executing an instruction sequence schedule process.

FIG. 4 is a flowchart of data processing.

FIG. 5: NDVI data diagram before and after LMF processing

FIG. 6 is a correlation diagram of data before and after LMF processing.

FIG. 7 is a data diagram of LHR processing of AVHRR / Ch4 data.

[FIG. 8] LMVI of NDVI of Pathfinder data
The figure which shows the result of F-KF processing.

FIG. 9: AVHRR / Ch4 brightness temperature data is LMF-
Diagram showing the result of KF processing

[Explanation of symbols]

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

   ─────────────────────────────────────────────────── ─── 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)

[Claims] 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. 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. 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. 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. 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. 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. 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. 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. 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. 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. 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.