CN113850139B

CN113850139B - Multi-source remote sensing-based forest annual phenological monitoring method

Info

Publication number: CN113850139B
Application number: CN202110986333.6A
Authority: CN
Inventors: 李明诗; 张亚丽; 孙敏
Original assignee: Nanjing Forestry University
Current assignee: Nanjing Forestry University
Priority date: 2021-08-25
Filing date: 2021-08-25
Publication date: 2022-04-08
Anticipated expiration: 2041-08-25
Also published as: CN113850139A

Abstract

The invention discloses a forest annual phenological monitoring method based on multi-source remote sensing, which comprises the steps of firstly, collecting all satellite remote sensing images with available cloud amount lower than 80%, and then correcting an integration method of different satellite remote sensing images to improve the space and spectrum matching degrees of different sensors; then, an improved continuous change detection and classification model is applied to generate a daily vegetation index curve; and finally, based on the daily synthetic image, adopting a logistic regression model to check the enhanced vegetation index, the normalized vegetation index and the surface water body index to extract the optimal forest annual SOS. The invention improves the integration method of different satellite data and increases the observation frequency; providing an MCCDC model, taking radiation difference into consideration, optimizing a model algorithm, shortening calculation time while ensuring precision, and finally generating a daily clear cloud-free remote sensing image; and 3, introducing a planted vegetation index to estimate the forest annual SOS, and estimating the difference of different indexes in estimating the forest SOS.

Description

Multi-source remote sensing-based forest annual phenological monitoring method

Technical Field

The invention belongs to the technical field of forest phenological monitoring, and particularly relates to a forest annual phenological monitoring method based on multi-source remote sensing.

Background

In recent decades, remote sensing gradually becomes an effective means for monitoring forest climate dynamics, aiming at the limitation that a ground climate observation network is difficult to monitor and map on a large scale. Some low-spatial-resolution sensors are widely applied to the extraction of the phenological information due to the fact that the low-spatial-resolution sensors have high time resolution, such as an Advanced Very High Resolution Radiometer (AVHRR) and a medium resolution imaging spectrometer (MODIS), have high time resolution, can acquire remote sensing data with the spatial resolution of 500 meters to 1100 meters, have the advantage of wide-range observation, and domestic and foreign scholars provide numerous model algorithms for estimating the start date of forest growing season (SOS) on the basis. However, when the low spatial resolution sensor is observed in a region-oriented manner, in a region with strong surface heterogeneity and high degree of fragmentation, a large number of mixed pixels are included in an image, so that a serious estimation error is caused; medium spatial resolution sensors such as Landsat (NASA terrestrial satellite in usa) have a spatial resolution of 30 meters and a long time span, which can overcome the above limitations, however, its 16-day revisit period and common cloud interference make it difficult for Landsat to accurately estimate forest year SOS.

In recent years, two satellites (Sentinel-2A and Sentinel-2B) respectively emitting Sentinel-2 in European space bureau are similar to Landsat spectral resolution, which provides a feasible scheme for improving the time resolution of vegetation phenological observation. The combination of Landsat and Sentinel-2 provides an average observation interval of 2.9 days, greatly increasing the ability to identify forest SOS on a regional scale, and NASA scientists currently provide a suite of integrated products of Landsat and Sentinel-2, but studies have shown that the integrated products are not applicable to all regions, and the integrated products are not shareable in most regions of China, thus requiring re-modification of the integration process to create new integrated products. In addition, for areas that are persistent cloudy, these observations are far from adequate for accurate assessment of forest SOS.

At present, most remote sensing phenological researches are carried out on the basis of a sensor with high time frequency and coarse spatial resolution, and are not suitable for accurate depiction of forest phenological variables in regional scales. Sensors with medium spatial resolution such as Landsat are often used to characterize forest average SOS over long time series, and extracting annual SOS becomes a big challenge. Increasing observation frequency or synthesizing a noiseless image becomes an effective method for remotely sensing and monitoring vegetation year forest SOS with medium spatial resolution.

Melaas et al, which increase the observation frequency to detect vegetation phenology in areas where multiple Landsat orbit numbers overlap, are generally not suitable for large-scale phenology estimation. Bolton et al found that integrating Landsat 8 and Sentinel-2 time series improved the accuracy of the landscape scale evaluation vegetation phenology model compared to MODIS-based monitoring evaluation.

The synthetic noiseless image solves the problem of cloud and shadow interference in a cloudy and rainy area, and is one of feasible strategies for improving the estimation accuracy of the forest year SOS. In addition, some researchers have tried to create high-temporal-resolution images from MODIS using image fusion techniques and existing Landsat images, for example, Bhandari et al have synthesized unacquired or inappropriate Landsat images using an enhanced spatio-temporal adaptive reflectivity fusion model (ESTARFM) algorithm, and finally generated vegetation index time series with 8-day intervals to evaluate forest phenology. However, this algorithm inevitably generates uncertainties in the image fusion process, and these uncertainties are particularly common in highly heterogeneous landscapes.

In summary, the following problems exist for forest annual phenological monitoring:

(1) at present, Landsat and Sentinel-2 integrated products developed by scientists of the United states space agency are not completely provided in China, so that the improvement of the integration method of Landsat and Sentinel-2 becomes one of the problems to be solved in the current climate monitoring.

(2) The premise of accurate estimation of the CCDC model is that a large number of clear Landsat observed values are required, if continuous rainy and cloudy weather occurs, the problem of over-fitting or under-fitting of the model can occur, the model has no universality, and the CCDC model does not consider radiation difference among different sensors.

(3) The CCDC model utilizes the first 12 clear observation values to establish an initial time sequence to simulate an SR motion track, tracks the SR change period, compares the model predicted value with the actual observation value, judges whether the land coverage changes or not by setting a threshold value, iteratively updates model parameters, and predicts the Landsat Surface Reflectivity (SR) of any expected date. Therefore, the CCDC model judges each observation value pixel by pixel, iteratively updates model parameters, has large calculation amount and low working efficiency, and is time-consuming for comparing the actual observation value with the model prediction value pixel by pixel aiming at a land block with unchanged ground surface coverage type in the Landsat scene.

Disclosure of Invention

The technical problems solved by the invention are as follows: the invention provides a forest inter-year phenological monitoring method based on multi-source remote sensing, which has universality, high calculation efficiency and accurate prediction result.

The technical scheme is as follows: in order to solve the technical problems, the technical scheme adopted by the invention is as follows:

a forest annual phenological monitoring method based on multi-source remote sensing comprises the steps of firstly, collecting all satellite remote sensing images with available cloud amount lower than 80%, and then correcting an integration method of different satellite remote sensing images to improve space and spectrum matching degrees of different sensors; then, an improved continuous change detection and classification model is applied to generate a daily vegetation index curve; and finally, based on the daily synthetic image, adopting a logistic regression model to check the enhanced vegetation index, the normalized vegetation index and the surface water body index to extract the optimal forest annual SOS. The method comprises the following specific steps:

s1, acquiring data and preprocessing the data;

s2, establishing a continuous change detection and classification model, and generating a daily vegetation index curve by using the continuous change detection and classification model;

s3: inputting the fitting coefficient of the established continuous change detection and classification model of each pixel into a support vector machine classifier to perform forest coverage mapping;

s4: and extracting the information of the start date of the annual forest growing season.

Further, in step S1, the data includes land satellite images of Landsat ETM +/OLI, Sentinel-2 images, forest resource second-class survey data, and phenological measurement data;

further, in step S1, the preprocessing of the data includes atmospheric correction, cloud and shadow masks, image registration, and image spectrum normalization processing.

Further, in step S2, the continuous change detection and classification model considers different change types of the pixels, fits the model on the basis of obtaining sufficiently clear pixels, gives a judgment standard to determine whether there is a land cover change in the observation period, and if so, judges each observation value one by one, identifies the specific date on which the land cover change occurs, and finally generates a daily remote sensing image through the given model.

Further, in step S3, a radial basis function is selected as a kernel function, a cost parameter C is defined, and a parameter γ is set as a reciprocal of the number of input features, so as to obtain a plurality of coefficients as input features for the support vector machine classification.

Further, in step S4, daily earth surface reflectance data is generated for the cloud and shadow areas, and fused with clear Landsat and Sentinel-2 observation values, a normalized vegetation index, an enhanced vegetation index and earth surface water are calculated, a logistic regression model is used to analyze 3 a pixel-by-pixel fluctuation mode of the planted vegetation index in the whole observation period, and a date corresponding to the maximum value of the change rate of the fitting model at the rising stage is used as the forest annual phenology.

Has the advantages that: compared with the prior art, the invention has the following advantages:

(1) according to the method, the Sentinel-2 data and Landsat series image data are integrated by adopting a linear model (slope and slope) to normalize and minimize the radiation difference of multi-source data, a Landsat image with higher time resolution is generated, the period is about 2.9 days, and the Landsat image is more accurate and reliable than a Landsat image synthesized by the conventional ESTARFM algorithm.

(2) The CCDC model is improved, on the basis of the CCDC model, multi-source data are introduced to increase the observation frequency, and a proper judgment standard is provided to determine whether the land coverage changes in the observation period. If yes, each observation is judged one by one, and the occurrence date of the land cover change is positioned, so that the calculation time is shortened while the accuracy is ensured, and the calculation efficiency is improved.

(3) Most of cloud and shadow noises are eliminated through integration of Landsat and Sentinel-2 and synthesis of clear and non-cloud images, so that the phenological monitoring is no longer affected by extraction of forest phenological information in a multi-cloud rainy region, and the phenological monitoring method is more universal in most regions. In addition, the method can expand the regional scale forest SOS mapping product to 30m spatial scale and has the direct potential of expanding the regional scale SOS mapping product to 10 m observation scale.

(4) The method comprises the steps of extracting forest annual SOS by using three vegetation indexes, namely a normalized vegetation index (NDVI), an Enhanced Vegetation Index (EVI) and a surface water body index (LSWI), and obtaining more accurate phenological information through preferential comparison.

Drawings

FIG. 1 is a flow chart of a forest annual phenological monitoring method based on multi-source remote sensing;

FIG. 2 is a schematic diagram of the MCCDC model;

FIG. 3 is a visual comparison of a composite image (right) with an actual image (left);

in fig. 4, (a) is a forest coverage map from 2013 to 2019; (b) forest area changes from 2013 to 2019;

FIG. 5 is a validation result of three vegetation index estimates SOS;

Detailed Description

The present invention will be further illustrated by the following specific examples, which are carried out on the premise of the technical scheme of the present invention, and it should be understood that these examples are only for illustrating the present invention and are not intended to limit the scope of the present invention.

The embodiment discloses a forest annual phenological monitoring method based on multi-source remote sensing, and a flow chart of the method is shown in figure 1: firstly, collecting all Landsat and Sentinel-2 images with the available cloud amount lower than 80%, and then correcting the integration method of Landsat and Sentinel-2 to improve the spatial and spectral matching degree of different sensors. A Modified Continuous Change Detection and Classification (MCCDC) model was then applied to generate a daily vegetation index curve at 30m spatial resolution. And finally, based on the daily synthetic image, adopting a logistic regression model to check an Enhanced Vegetation Index (EVI), a normalized vegetation index (NDVI) and a surface water body index (LSWI) to extract the optimal forest annual SOS. The method specifically comprises the following steps:

s1, acquiring satellite remote sensing data and preprocessing the data;

firstly, data are obtained, and the data comprise Landsat ETM +/OLI land satellite images, Sentinel-2 images, forest resource two-class survey data and phenological measurement data.

(1) Land satellite images of Landsat ETM +/OLI with cloud content less than 80% from 1/2013 to 12/30/2019 are downloaded from a United States Geological Survey (USGS) website, and the downloaded data track number is 131/42, and 188 scenes are counted. The ETM + image is 89 shots and the OLI image is 99 shots. The basic information of the selected Landsat images is shown in table 1, wherein Landsat ETM +/OLI are both surface reflectance products;

(2) the method comprises the steps of downloading Sentinel-2 atmospheric apparent reflectance (TOA) radiance data with cloud amount less than 80% in 3-11 months in 2017 and 2019 from an European aviation administration (ESA) website, wherein the used spectral bands correspond to Landsat (track numbers T47RPJ, Sentinel-2A58 scene and Sentinel-2B38 scene), and the basic information is shown in Table 1;

(3) a forestry department provides 2016 research data of forest resources of a research area for model precision verification;

(4) and 40 phenological observation plots, wherein 28 observation data are obtained from field worker records, and 12 observation data are obtained by field investigation of an author and are used for model precision verification.

Table 1 basic information of remote sensing image

Then, preprocessing the acquired data, wherein the data preprocessing comprises the following steps:

(1) atmospheric correction: the SR data of Landsat is provided directly by USGS, Landsat ETM + uses the LEDAPS algorithm, and OLI images are corrected for atmospheric pressure using the LaSRC algorithm. The Sentinel-2 data is provided by the ESA, and the TOA is calibrated to SR products using a Sen2cor radiation tool.

(2) Cloud and shadow masks: clouds and snow have high brightness in the visible band, while cloud shadows are low brightness and are easily confused with other types of land cover, causing interference. Therefore, the invention uses Fmask 4.0 software to remove the cloud and the corresponding shadow interference source from the Landsat and the Sentinel-2 images.

(3) Image registration: landsat and Sentiniel-2 data do not match in spatial position and require image registration. And (3) manually selecting a plurality of ground control points in the overlapped area of Landsat and Sentinel-2 by taking Sentinel-2 as a reference image, fitting a first-order polynomial to perform batch registration, wherein the geometric error of the registration is less than 0.5 pixel.

(4) Image spectrum normalization processing of Landsat and Sentiniel-2: landsat and Sentiniel-2 sensors have different spectral responses at specific wavelengths, and therefore band-to-band adjustments are needed to eliminate differences in radiation due to sensor differences.

The method adopts a general coefficient and c factor method to carry out BRDF correction, then identifies a plurality of pixel samples on image pairs with the same or adjacent acquisition date (ensuring that the image pairs do not have land coverage change), then fits a linear model (slope and intercept) to carry out wave band-by-wave radiation normalization operation, and resamples the Sentinel-2 image to 30m by adopting a nearest neighbor point resampling method.

The invention normalizes and minimizes the radiation difference of multi-source data by modifying the integration method of Sentinel-2 data and Landsat data and adopting a linear model (slope and slant range) to generate a more consistent synthetic image.

the MCCDC model considers different change types of the pixels, fits the model (formula 1) on the basis of obtaining enough clear pixels, and gives a judgment standard to determine whether the land coverage changes in an observation period. If yes, judging each observation value one by one, identifying the specific date of land cover change, and finally generating a daily remote sensing image through a given model.

Firstly, sequencing all clear pixels in an observation period at the same position of each wave band according to image acquisition time, fitting a complete time sequence model based on a formula (1), and then calculating a model R²The value is obtained. R²A value of greater than or equal to 0.6 can be used to determine that the land cover is unchanged; if the land cover changes, the fitted curve will deviate from the actual observed value to a great extent, resulting in R²Low in value (<0.6). If the land cover change occurs at the beginning or end of the entire observation period, R²It is also very large, so the invention sets that if the absolute value of the difference between the 4 consecutive predicted values and the actual observed value exceeds 3 times the model RMSE, then it is also assumed that a land cover change has occurred. In summary, if any condition of equation (2) is satisfied, the pixel position occurs in the whole observation periodAnd (4) changing the land coverage.

After a land cover change event is determined according to the formula (2), the occurrence of the land cover change event on a specific change date needs to be accurately extracted. For the pixel with the land cover change, the method is based on the basic principle of a CCDC model, the first 24 clear pixels at the same position are fitted with an initialized time sequence model, then whether the difference between the subsequent 4 continuous observed values and the predicted value of the pixel exceeds 3 times of the Root Mean Square Error (RMSE) of the model or not is judged, and if the difference exceeds the range, the 1 st occurrence time is defined as the occurrence time of the land cover change.

In the formula: x is julian date, T is the number of days in a year, a_0,i，a_1,i，a_2,i，b_1,i，b_2,i，c_1,i，c_2,i，d_1,iThe model coefficients to be fitted are obtained; i is the band number, k is the total number of bands, here the value 6,

is the SR predicted value of the ith wave band in the julian date x; ρ (i, x) is the actual SR observed value of the ith band in julian date x; RMSE_iIs the root mean square error, R, of the predicted value and the actual observed value of the ith band²Is the decision coefficient of the model; expressing the degree of reliability of the original observed data, R, represented by the fitted model²The larger the value, the better the performance of the model to fit the original observed data.

Figure 2 shows the NIR and SWIR1 band fit curves for 2 picture elements at different locations and it can be seen that MCCDC fits well to the regular change in SR when considering the occurrence of land cover change, where the fit curves in (a) and (b) indicate that the land cover type of the first picture element has not changed and the fit curves in (c) and (d) indicate that the land cover of the second picture element has changed during the observation.

S3: drawing forest coverage

The MCCDC model fitting coefficient of each pixel is input into a Support Vector Machine (SVM) classifier to perform forest coverage mapping, a Radial Basis Function (RBF) is selected as a kernel function, a cost parameter C is defined as 10, a parameter gamma is usually set as the reciprocal of the number of input features, and 48 coefficients (8 coefficients in each wave band and 6 wave bands) are used as the input features of SVM classification.

S4 extraction of information on start date of annual forest growing season (SOS)

Generating daily earth surface reflectivity data aiming at cloud and shadow areas, fusing the daily earth surface reflectivity data with clear Landsat and Sentinel-2 observation values, calculating a normalized vegetation index (NDVI), an Enhanced Vegetation Index (EVI) and an earth surface water body (LSWI), analyzing a pixel-by-pixel fluctuation mode of a planting index in the whole observation period by adopting a logistic regression model 3, and taking the date corresponding to the maximum value of the change Rate (ROC) of a fitting model in the rising stage as a forest SOS.

Where ρ is_Blue、ρ_Red、ρ_NIRAnd ρ_swir1Surface reflectance values for blue, red, near infrared and short wave infrared 1 bands, respectively. t is the observation time, y (t) is the vegetation index value, a and b are the fitting parameters, the sum of c + d is the maximum vegetation index value, d is the initial background vegetation index value, and y (t)' are the first and second derivatives of the logistic regression equation, respectively.

S5 evaluation of accuracy

By calculating R²And the RMSE method evaluates the precision of the synthesized image and the precision of the forest SOS, and adopts a confusion matrix to calculate the user precision, the producer precision and the overall precision of the forest and non-forest to evaluate the precision of forest coverage mapping.

The method has the following advantages: (1) improving the integration method of Landsat and Sentinel-2 and increasing the observation frequency; (2) an MCCDC model is provided, radiation difference between Landsat 7ETM + and Landsat 8OLI is taken into consideration, a model algorithm is optimized, the calculation time is shortened while the precision is guaranteed, and finally a daily clear cloud-free remote sensing image is generated; (3) and 3, introducing a planted vegetation index to estimate the forest annual SOS, and estimating the difference of different indexes in estimating the forest SOS.

The forest annual phenological model of the embodiment is subjected to result analysis

(1) Integration and composite Landsat image evaluation of multi-source data

Table 2 shows the slope, intercept and R of the SR match of Landsat 8 with Landsat 7 and Sentiniel-2². First, R of Blue band between Landsat 8 and Sentinel-2²The lowest, value is 0.83. R between SWIR2 bands²Maximum, the value is 0.94. Landsat 8 and Landsat 7 are also R for Blue band²R between the SWIR2 bands with the lowest value²The highest value, the value is 0.95. In general, all R²Values were all greater than 0.80, demonstrating successful band-to-band spectral radiance normalization between Landsat series and Sentinel-2.

TABLE 2 band-by-band spectral normalization transform coefficients for Landsat and Sentinel-2

Table 3 shows R of the synthesized band of Landsat 7/8 and the corresponding real band²R in the blue and green bands²The values are lower than the other bands. Except for red band, each band R in 2013 in summer²All values are lower than in other seasons. Winter R²R of highest, red, near infrared and SWIR1/2 band²All values are greater than 0.90. Each wave band R in spring and winter of 2019²The values are not very different, both higher than the other two seasons (table 3).

TABLE 3 bands R of synthetic and real images²Value of

By visually comparing the original image with the composite image in 2013 winter and 2019 spring (fig. 3), the composite image is very similar to the original image for different land cover types, and more importantly, the composite image completely overcomes the effects of noise, such as clouds, shadows and SLC-off image gaps of Landsat 7. The method can well depict some characteristic ground objects such as lakes, roads, forests and buildings.

(2) Dynamic change of forest coverage

By comparing 2016 forest resource secondary survey data with 2016 forest coverage classification maps, the result shows that the user precision of the forest is 89%, the producer precision is 88%, and the overall precision is 88%. On the basis, the Landsat image data is used for verifying the 2019 forest coverage classification map, and the result shows that: the user precision, the producer precision and the overall precision of the forest and the non-forest are all larger than 88%. Overall, these accuracies indicate that the forest coverage maps obtained by the prior art methods are authentic (table 4). The forest area dropped from 4015.69 square kilometers in 2013 to 4010.17 square kilometers in 2019 (fig. 4).

TABLE 4 forest coverage mapping accuracy evaluation

(3) Forest SOS extraction

The method collects forest SOS information from 40 samples in 2013, 2018 and 2019, compares SOS extracted based on remote sensing with field observation data, and obtains R of forest SOS estimated by different vegetation indexes²And RMSE value (figure 5), and carrying out forest SOS distribution mapping to obtain forest climate dynamic change information (figure 6).

As seen in FIG. 5, R of EVI predicted SOS and measured SOS²The lowest, value was 0.26 and the highest RMSE was 19.30 days. NDVI and LSWI predicted SOS versus observed R²Values were 0.61 and 0.66, respectively, corresponding to an RMSE of about 11 days.

Fig. 6(a) shows the SOS distribution extracted from the three indices in 2013 and 2019. Fig. 6(b) and (c) show annual SOS dynamics for two different sites from 2013 to 2019. NDVI and LSWI extract SOS with similar spatial and temporal distribution, most SOS mainly from 80-150 days. Specifically, the western SOS of the study area extracted from NDVI and LSWI was dominated by 80-120 days, and the SOS of the bottom and southeast corner was dominated by 120-150 days. The north SOS was later than 150 days (fig. 6 (a)). In general, the south SOS is earlier than the north SOS, and the west SOS is earlier than the east SOS. For some study areas, the maximum SOS difference for the three indices did not exceed one week (fig. 6 (b)). However, in other small parts of the study area, the SOS of different indices varies greatly. When the EVI-based SOS was less than 40 days, the SOS from NDVI and LSWI was later than 90 days, with a difference of up to two months (fig. 6 (c)).

The improved model is more complex (more parameters) than the existing model, and not only the trend term but also the annual variation are considered; according to the method, the Sentinel-2 data is introduced, the density of observation points can be obviously increased, the phenomenon that the accuracy of a CCDC model is low in perennial cloudy and rainy areas is made up, and the model is more universal; the MCCDC model provided by the invention makes a trade-off between model precision and calculation efficiency, and reduces the calculation time while ensuring certain precision. Traditional CCDC model calculation may take several days, but the MCCDC model only needs tens of hours to complete; the CCDC model was built without taking into account the sensor differences of Landsat5, 7, and 8. In fact, the reflectivity of different features may vary. The difference in reflectivity between them is eliminated by linear regression herein, resulting in a more accurate composite image.

The foregoing is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, various modifications and decorations can be made without departing from the principle of the present invention, and these modifications and decorations should also be regarded as the protection scope of the present invention.

Claims

1. A forest annual phenological monitoring method based on multi-source remote sensing is characterized by comprising the following steps:

s1, acquiring data and preprocessing the data;

collecting all Landsat and Sentinel-2 satellite remote sensing images with available cloud amount lower than 80%, and then correcting the integration method of different satellite remote sensing images to improve the space and spectrum matching degree of different sensors;

performing BRDF correction by adopting a general coefficient and c factor method, identifying a plurality of pixel samples on an image pair with the same or adjacent acquisition date to ensure that the image pair has no land coverage change, fitting a linear model to perform wave band-by-band radiation normalization operation, and resampling the Sentinel-2 image to 30m by adopting a nearest point resampling method;

the continuous change detection and classification model considers different change types of pixels, the model is fitted on the basis of obtaining sufficiently clear pixels, a judgment standard is given to determine whether the land coverage change exists in an observation period, if yes, each observation value is judged one by one, specific dates of the land coverage change are identified, and finally, daily remote sensing images are generated through the given model;

firstly, sequencing all clear pixels in an observation period at the same position of each wave band according to image acquisition time, fitting a complete time sequence model based on a formula (1), and then calculating a model R²Value R²A value of greater than or equal to 0.6 can be used to determine that the land cover is unchanged; if the land cover changes, the fitted curve will deviate from the actual observed value to a great extent, resulting in R²The value is low, if the land cover change occurs at the beginning or end of the whole observation period, R²It is also very large; setting that if the absolute value of the difference value between the continuous 4 predicted values and the actual observed value exceeds 3 times of the model RMSE, assuming that the land coverage change occurs; therefore, if any condition of the formula (2) is satisfied, the pixel position is changed in land coverage in the whole observation period;

after a land cover change event is determined according to the formula (2), the occurrence condition of the land cover change event on a specific change date needs to be accurately extracted, for a pixel with the land cover change, the first 24 clear pixels at the same position are fitted with an initialized time sequence model, then whether the difference between the subsequent 4 continuous observed values and the predicted value of the pixel exceeds 3 times of the root mean square error of the model or not is judged, and if the difference exceeds the range, the 1 st occurrence time is defined as the occurrence time of the land cover change;

is the SR predicted value of the ith wave band in the julian date x; ρ (i, x) is the actual SR observed value of the ith band in julian date x; RMSE_iIs the root mean square error, R, of the predicted value and the actual observed value of the ith band²Is the decision coefficient of the model; expressing the degree of reliability of the original observed data, R, represented by the fitted model²The larger the value, the better the performance of the model to fit the original observed data;

s4: extracting information of the beginning date of the growth season of the annual forest;

and (3) based on the daily synthetic image, adopting a logistic regression model to check the enhanced vegetation index, the normalized vegetation index and the surface water body index to extract the optimal forest annual SOS.

2. The forest annual phenological monitoring method based on multi-source remote sensing according to claim 1, characterized in that: in step S1, the data includes Landsat ETM +/OLI land satellite images, Sentinel-2 images, forest resource survey data and phenological measurement data.

3. The forest annual phenological monitoring method based on multi-source remote sensing according to claim 1, characterized in that: in step S1, the preprocessing of the data includes atmospheric correction, cloud and shadow masks, image registration, and image spectral normalization processing.

4. The forest annual phenological monitoring method based on multi-source remote sensing according to claim 1, characterized in that: in step S3, a radial basis function is selected as a kernel function, a cost parameter C is defined, and a parameter γ is set as the reciprocal of the number of input features, so as to obtain a plurality of coefficients as input features for support vector machine classification.

5. The forest annual phenological monitoring method based on multi-source remote sensing according to claim 1, characterized in that: in step S4, daily earth surface reflectance data is generated for cloud and shadow areas, and fused with clear Landsat and Sentinel-2 observation values, a normalized vegetation index, an enhanced vegetation index and an earth surface water body are calculated, a pixel-by-pixel fluctuation mode of the planting index in the whole observation period is analyzed 3 by a logistic regression model, and a date corresponding to the maximum value of the change rate of the fitting model in the rising stage is used as a forest annual phenology.