Forest tree species group classification method considering harmonic model coefficients and phenological parameters

Abstract

The present invention discloses a forest tree species group classification method considering harmonic model coefficients and phenological parameters, comprising the following steps: S1, image acquisition and preprocessing; S2, establishing an MCCDC model based on multi-source remote sensing data; S3, extracting annual forest phenology and MCCDC coefficients based on a logistic regression equation; S4, identifying tree species groups considering phenological parameters and MCCDC parameters. The present invention uses forest phenological parameters extracted from remote sensing and the constructed harmonic model MCCDC coefficients as input variables for tree species classification modeling, providing unique classification data features different from traditional spectral and texture feature methods, compensating to a certain extent for the influence of the phenomena of "same object with different spectra" and "different objects with the same spectra" on the classification accuracy when only relying on spectral and texture feature inputs, thereby improving the classification accuracy of forest tree species (groups).

Description

Technical Field

[0001] The present invention belongs to the field of forestry remote sensing technology, and specifically relates to a method for taking into account the variation characteristics of the long-term series reflectivity harmonic model simulation coefficients of different ground objects in multi-source remote sensing data, and extracting vegetation phenological parameters based on the same to perform forest tree species (groups) identification and spatial mapping. Background Art

[0002] The distribution of forest tree species (groups) is one of the basic information required by forestry management departments, wildlife protection organizations, etc. The collection of accurate spatial distribution information of forest tree species (groups) is particularly important for the development of sustainable forest management strategies and the achievement of the country's "dual carbon" goals. Traditional forest tree species (groups) identification mainly relies on large-scale manual surveys. Obviously, this has the defects of large time and manpower investment, large terrain restrictions, and poor timeliness. With the continuous development of remote sensing technology, visual interpretation of high-spatial-resolution aerial photos or commercial satellite data to extract tree species distribution information is also a way, but this method is usually costly, has a small scope of application, and has high requirements for the knowledge and experience, professional background and patience of image interpreters. Later, with the advancement of computer algorithms and the increasing popularity of satellite data, the current main method is to couple various machine learning algorithms and deep learning algorithms, with the help of input features such as remote sensing spectral signals, texture features and terrain variables, to automatically classify and identify forest tree species (groups), in order to improve the efficiency and accuracy of forest tree species classification mapping.

[0003] At present, there have been some studies on the identification and classification of forest tree species (groups), including experimental studies based on low spatial resolution remote sensing image data (such as MODIS), medium spatial resolution remote sensing data (such as Landsat, Sentinel-2), high spatial resolution remote sensing data (such as WorldView, QuickBird) and other multispectral remote sensing data. Among them, the classification results based on low spatial resolution data are often rough due to the presence of a large number of mixed pixels, which limits its accuracy and applicability in mapping forest types in fragmented landscapes; the classification results based on medium resolution multispectral remote sensing data are mostly based on the classification of single tree species or the distinction between coniferous, broad-leaved and mixed tree species groups, and often cannot be divided into more detailed tree species; when classifying tree species based on high-resolution remote sensing data, the crown characteristics of different tree species can be better used (expressed by texture measurements derived from high-resolution remote sensing data) to identify tree species, but it often has high time and financial costs, and its practicality on a large regional scale will be limited. Under the influence of the extensive spatiotemporal heterogeneity of surface features, these studies on the classification and identification of forest tree species based on multispectral remote sensing data and its derived texture features and terrain parameters are still widely plagued by the "same object, different spectrum" and "different objects, same spectrum" phenomena, and the reliability of tree species (group) classification results is still far from practical application. In the past 20 years, hyperspectral remote sensing has shown certain application potential in the classification and identification of forest tree species (groups), but it is greatly restricted by cost and data processing complexity. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to address the deficiencies in the prior art and to provide a method for identifying and classifying forest tree species groups by introducing into the classifier, in terms of input features, coefficients of an improved harmonic model for continuous change detection and classification (MCCDC), as well as phenological parameters of major tree species characterized by the harmonic model.

[0005] Technical solution: The forest tree species group classification method considering harmonic model coefficients and phenological parameters of the present invention comprises the following steps:

[0006] S1. Image acquisition and preprocessing: obtain multiple sensor satellite images of the area to be classified within a period of time, and use the least squares method to fit the slope and intercept of the obtained images to perform spectral normalization band by band;

[0007] S2. Establish a MCCDC model based on multi-source remote sensing data. The specific steps include:

[0008] S21, establish a SR time series model, define a pixel-by-pixel reflectance time series model, as shown in formula (1), to express the characteristics of surface reflectance change;

[0009]

[0010] Formula (1) contains a Fourier harmonic model and a long-term trend part, x is the Julian date; T is the cycle day, the value is 365.25; is the predicted value of the surface reflectance of the i-th band on the Julian date x; there are 8 coefficients in this time series model, coefficient a 0,i is the overall mean of clear observations in the i-th band, and the coefficient a 1,i 、a 2,i 、b 1,i 、b 2,i 、c 1,i 、c 2,i It is used to simulate the annual variation of surface reflectance caused by phenology. Its minimum positive period is one year, and the coefficient d 1,i Reflects the trend of gradual changes in surface reflectivity caused by factors such as climate change and soil degradation;

[0011] S22, identifying the mutation point of the time series model, based on the time series model established in step S21, to see whether it satisfies formula (2);

[0012]

[0013] In formula (2), k is the number of bands; n refers to the number of consecutive observations after the mutation point to be judged; ρ(i,x) is the actual observation value of the surface reflectance of the i-th band on the Julian date x; RMSE i It refers to the root mean square error between the model-predicted reflectivity value and the actual remote sensing observed reflectivity value in the i-th band;

[0014] When formula (2) is satisfied, it proves that land cover change has occurred in the area. The changed pixels are sorted in chronological order. First, the first 24 clear pixels at the same location are used to establish an initialized time series model based on formula (1). Then, it is determined whether the difference between the actual observation value and the predicted value of the next 4 pixels in each band exceeds 3 times the root mean square error of the model. If it does not exceed the range, these observations are added to the observations previously involved in the modeling, the parameters are refitted, and the subsequent observations are judged again. If 4 consecutive observations are out of range, the point is defined as the mutation point, that is, the time point when the land cover change occurs. In this process, it is necessary to average the sum of the differences between the observation value and the predicted value of all bands (as shown in formula (2)) to define the mutation;

[0015] S23, synthesizing daily remote sensing images, using the method of steps S21 to S22 to generate time series models of different bands for each pixel, for a single band of each pixel, by giving the value of the variable Julian date x, according to the time series model, generating the SR prediction value of the pixel on a given date, applying this method to all pixels of each day within a period of time, thereby generating a clear synthetic image with a spatial resolution of 30m;

[0016] S3. Extraction of forest annual phenology based on logistic regression equation. The specific method is as follows:

[0017] S31. Based on all the daily composite images obtained, the annual phenological parameters are extracted year by year, and the three vegetation indices EVI, NDVI and LSWI are calculated in turn by equations (3) to (5), and the logistic regression functions are fitted based on the three respectively, so as to determine the estimated values ​​of different phenological parameters;

[0018]

[0019]

[0020]

[0021] In the formula, ρ Blue , Red , NIR and They are the surface reflectance values ​​in the blue, red, near-infrared and short-wave infrared 1 bands;

[0022] S32. Based on the three vegetation indices, the annual SOS, EOS and LOS expressed by each were extracted using the logistic regression function, and the annual vegetation index time series data were modeled using formula (6);

[0023]

[0024] Where t is the observation time, y(t) is the vegetation index value (EVI, or NDVI, or LSWI), a and b are fitting parameters, c+d is the maximum vegetation index value, and d is the initial background vegetation index value;

[0025] After completing the fitting of formula (6), the curvature change rate of the fitting model is calculated by formula (7);

[0026]

[0027] Where y(t)' and y(t)" are the first and second order derivatives of the logistic equation, respectively, and ROC is the rate of change of curvature of the fitted model;

[0028] Criteria for judging phenological parameters: When the curvature change rate reaches the maximum value for the first time, the corresponding number of days is the vegetation growth start period SOS; when the curvature change rate reaches the minimum value for the last time, the corresponding number of days is the vegetation growth end period EOS; the difference between EOS and SOS is the length of the vegetation growth season LOS;

[0029] S4. Considering the phenological parameters and the coefficients of the fitted MCCDC model, the tree species groups were identified by using RF or SVM recognition algorithms to classify the forest tree species groups respectively. The parameters of the MCCDC model and the phenological period parameters were added as input variables to the two machine learning algorithms to draw the spatiotemporal distribution map of the forest tree species groups.

[0030] A further preferred technical solution of the present invention is:

[0031] In step S3, the phenological parameters obtained by the LSWI index are selected as input variables of the recognition algorithm.

[0032] Preferably, in step S4, a SVM recognition algorithm is selected to classify forest tree species groups.

[0033] Beneficial effects: In the process of establishing a forest tree species classifier, the present invention introduces MCCDC model coefficients and phenological parameters as classification basis, which is different from the unique classification basis of traditional spectral and texture feature methods. It compensates to a certain extent for the influence of "same object different spectrum" and "different objects same spectrum" phenomena on classification accuracy when relying solely on spectral, texture and terrain input, thereby improving the classification accuracy of forest tree species (groups); the method of the present invention does not use a classification method that combines medium spatial resolution remote sensing image data with higher resolution image data, but only uses medium spatial resolution remote sensing image data, and uses the MCCDC model to synthesize daily remote sensing images to classify forest tree species; and introduces MCCDC model coefficients and phenological parameters in the classifier to classify tree species. The results show that the method of the present invention is better than the traditional classification method that only uses spectral features and texture features, and the accuracy of forest tree species recognition is improved.

[0034] In step S1, the present invention performs a band-by-band spectral normalization operation on the acquired land satellite image, minimizing the spectral response differences between different sensors, thereby improving the MCCDC model fitting accuracy; in step S22, the present invention uses the subsequent four continuous observations to define the real land cover changes, significantly reducing the risk of the model capturing false land cover changes; the present invention introduces harmonic model coefficients and phenological parameters into the classification model, significantly reducing the impact of "same object different spectrum" and "different objects same spectrum" caused by conventional classification features (spectral signal + texture feature) on classification accuracy. In step S3, the present invention selects the phenological parameters obtained by the LSWI index as the input variables of the recognition algorithm, because it has stronger regional adaptability and the extracted phenological parameters are more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 This is a flow chart of the identification and classification method of the present invention;

[0036] Figure 2 This is a comparison chart of the visual effects of Landsat 5, 7 and 8 four-season synthetic images in the embodiment;

[0037] Figure 3 It is the SOS and EOS accuracy verification map extracted based on EVI, NDVI and LSWI in the embodiment;

[0038] Figure 4 It is the spatial distribution map of forest SOS from 1999 to 2019 in the embodiment;

[0039] Figure 5 This is the spatial distribution map of forest EOS from 1999 to 2019 in the embodiment;

[0040] Figure 6 is the spatial distribution map of forest LOS from 1999 to 2019 in the embodiment;

[0041] Figure 7 Spectral curve characteristic diagram of different forest tree species over time in the embodiment;

[0042] Figure 8 A comparison diagram of phenological periods of different tree species groups in the embodiment;

[0043] Fig. 9 It is a comparison diagram of the tree species group distribution map in 2016 based on RF and SVM in the embodiment. DETAILED DESCRIPTION

[0044] The technical solution of the present invention is described in detail below with reference to the accompanying drawings, but the protection scope of the present invention is not limited to the embodiments.

[0045] Embodiment: A method for classifying forest tree species groups taking into account harmonic model coefficients and phenological parameters comprises the following steps:

[0046] S1, image acquisition and preprocessing;

[0047] S11. Obtain the land satellite images of multiple sensors in the area to be classified for a period of time. The data of this embodiment is downloaded for free from the official website of the United States Geological Survey (USGS) from January 1, 1999 to December 31, 2019. The land satellite images of Landsat TM / ETM+ / OLI with cloud cover less than 80% from January 1, 1999 to December 31, 2019. The downloaded data track number is 131 / 042, with a total of 610 images. There are 152 TM images, 243 ETM+ images, and 119 OLI images. Due to weather reasons, there are relatively few images with cloud cover less than 80% in July in summer. Sentinel-2 data is also downloaded from the European Space Agency. The data includes the atmospheric top reflectance (Top of Atmosphere, TOA) radiation data of Sentinel-2 with cloud cover less than 80% from March to November from 2017 to 2019. The track number is T47RPJ, including 58 2A data and 38 2B data. The basic information of the selected Landsat and Sentinel-2 images is shown in Table 1:

[0048] Sensor Type Satellite Number of images Data Level Start time End Time TM Landsat-5 152 SR 1999.01 2011.09 ETM+ Landsat-7 243 SR 1999.06 2017.05 OLI Landsat-8 119 SR 2013.04 2019.12 MSI Sentinel-2A 58 TOA 2017.04 2019.11 MSI Sentinel-2B 38 TOA 2017.10 2019.11

[0049] Table 1 Basic information of Landsat and Sentinel-2 remote sensing images used

[0050] S12. Preprocess the obtained image to obtain a super-resolution image: For the images provided by USGS, the present invention directly applies for SR images. The TM and ETM+ images use the LEDAPS algorithm for atmospheric correction to obtain SR images, while the OLI images use the LaSRC algorithm for atmospheric correction to obtain SR images. For the downloaded Sentinel-2 TOA image data, the Sen2Cor atmospheric radiation transfer correction code (version 2.8, released on February 21, 2019; http: / / step.esa.int / main / third-party-plugins-2 / sen2cor / ) is used for atmospheric correction.

[0051] S13. In order to minimize the spectral response differences between different sensors, the least squares method is used to fit the slope and intercept, and a spectral normalization operation is performed band by band.

[0052] S2. Establish a MCCDC model based on multi-source remote sensing data. The specific steps include:

[0053] S21. Establish SR time series model:

[0054] Based on Fourier's idea, any complex periodic function can be expressed as a weighted sum of sine and cosine functions with different frequencies. This embodiment first defines a pixel-by-pixel reflectance time series model, as shown in formula (1), to express the characteristics of surface reflectance changes;

[0055]

[0056] Formula (1) contains a Fourier harmonic model and a long-term trend part, x is the Julian date; T is the cycle day, the value is 365.25; is the predicted value of the surface reflectance of the i-th band on the Julian date x; there are 8 coefficients in this time series model, coefficient a 0,i is the overall mean of clear observations in the i-th band, and the coefficient a 1,i 、a 2,i 、b 1,i 、b 2,i 、c 1,i 、c 2,i It is used to simulate the annual variation of surface reflectance caused by phenology. Its minimum positive period is one year, and the coefficient d 1,i Reflects the trend of gradual changes in surface reflectivity caused by factors such as climate change and soil degradation;

[0057] S22, identifying the mutation point of the time series model, based on the time series model established in step S21, to see whether it satisfies formula (2);

[0058]

[0059] In formula (2), ρ(i,x) is the actual observed value of the surface reflectance corresponding to the Julian date x of the i-th band, k is the number of bands (6 in this example), n is the number of consecutive observations after the mutation point to be judged (4 in this example), and RMSE is i It refers to the root mean square error between the model-predicted reflectivity value and the actual remote sensing observed reflectivity value in the i-th band;

[0060] When formula (2) is satisfied, it is considered that the time series model has land cover changes during the entire observation period. In order to accurately determine the specific date when the change occurs, the principle of the CCDC model is used in this embodiment to sort the changed pixels in chronological order. First, the first 24 clear pixels at the same position are used to establish an initialized time series model based on formula (1). The reason for selecting 24 clear observations is that the time series model has 8 coefficients. The number of clear observations is 3 times the number of coefficients of the curve to be fitted, which can make the coefficients estimated by the model more accurate. Then, it is determined whether the difference between the actual observations and the predicted values ​​of the next 4 pixels in each band exceeds 3 times the root mean square error of the model. If it does not exceed the range, these observations are added to the observations previously involved in the modeling, the parameters are refitted, and the subsequent observations are determined again. If 4 consecutive observations are out of range, the point is defined as a mutation point, that is, the time point when the land cover change occurs; in order to make full use of all spectral information to accurately define the breakpoint, the input of the MCCDC model is six bands, so it is necessary to average the sum of the differences between the observations and predicted values ​​of all bands to define the "mutation".

[0061] S23, synthesize daily remote sensing images, use the methods of steps S21 to S23 to generate time series models of different bands for each pixel, and for a single band of each pixel, generate the SR prediction value of the pixel on a given date based on the time series model by giving the value of the variable Julian date X, and apply this method to all pixels of each day from 1999 to 2019, so as to generate a clear synthetic image with a spatial resolution of 30m, as shown in Figure 2 shown.

[0062] S3. Extraction of forest annual phenology based on logistic regression equation. The specific method is as follows:

[0063] S31. Based on all the daily composite images obtained, the annual phenological parameters are extracted year by year, and the three vegetation indices EVI, NDVI and LSWI are calculated in turn by equations (3) to (5), and the logistic regression functions are fitted based on the three respectively, so as to determine the estimated values ​​of different phenological parameters;

[0064]

[0065]

[0066]

[0067] In the formula, ρ Blue , Red , NIR and They are the surface reflectance values ​​in the blue, red, near-infrared and short-wave infrared 1 bands;

[0068] S32. Based on the three vegetation indices, the annual SOS, EOS and LOS expressed by each were extracted using the logistic regression function, and the annual vegetation index time series data were modeled using formula (6);

[0069]

[0070] Where t is the observation time, y(t) is the vegetation index value (EVI, or NDVI, or LSWI), a and b are fitting parameters, c+d is the maximum vegetation index value, and d is the initial background vegetation index value;

[0071] After completing the fitting of formula (6), the curvature change rate of the fitting model is calculated by formula (7);

[0072]

[0073] Where y(t)' and y(t)" are the first and second order derivatives of the logistic equation, respectively, and ROC is the rate of change of curvature of the fitted model;

[0074] Criteria for judging phenological parameters: When the curvature change rate reaches a maximum value for the first time, the corresponding number of days is the start of vegetation growth SOS; when the curvature change rate reaches a minimum value for the last time, the corresponding number of days is the end of vegetation growth EOS; the difference between EOS and SOS is the length of the vegetation growing season LOS.

[0075] By comparing the SOS and EOS extracted from the three remote sensing indices with the field observed phenological data, such as Figure 3 As shown in the figure, it was found that the R2 between the SOS predicted by EVI and the measured SOS was the lowest (0.26), and its RMSE was 19.30 days. The R2 values ​​between the SOS predicted by NDVI and LSWI and the observed values ​​were 0.61 and 0.66, respectively, and the corresponding RMSE was about 11 days. The R2 of EOS obtained based on EVI, NDVI and LSWI were not much different, which were 0.60, 0.62 and 0.64, respectively, and the RMSE were 14.85 days, 12.75 days and 12.04 days, respectively. Therefore, the LSWI index was selected for the subsequent dynamic mapping of forest phenology.

[0076] Figure 4The spatiotemporal variation pattern of SOS in the forests of the study area is shown. In general, SOS in the north of the study area is delayed than that in the south, and the occurrence time of SOS in the east of the study area is later than that in the west. There are mostly yellow and light green areas in the study area, indicating that most of the SOS in the study area are concentrated in the 85th to 135th day. Forests with SOS less than 85 days are mainly concentrated in the westernmost part of the study area. Forests with SOS greater than 135 days are mainly concentrated in the northeast corner of the study area. It can be seen from the SOS map that the area in the northeast corner of the study area greater than 135 days has been gradually increasing from 1999 to 2019, and the red area in the west of the study area has slightly decreased, which indicates that SOS in this part is delayed.

[0077] Figure 5 The spatiotemporal variation pattern of EOS in the forest of the study area is shown. Overall, the EOS in the western part of the study area is earlier than that in the eastern part, and the EOS in the northern part of the study area is earlier than that in the southern part. In addition, the green and blue areas in the study area are mostly present, indicating that most EOS in the study area is later than the 310th day. As can be seen from the figure, from 1999 to 2019, there is a trend of greening in the southern part of the study area, indicating a delay in EOS.

[0078] By comparing the differences between SOS and EOS, the spatial distribution map of LOS can be obtained. Figure 6 It can be seen that the forest LOS in the eastern part of the study area is generally less than 200 days, the LOS in the westernmost area of ​​the study area is concentrated above 250 days, and the LOS in the remaining areas of the study area is mainly concentrated between 200 and 225 days.

[0079] S4. Identification of tree species groups considering phenological parameters and MCCDC parameters:

[0080] The choice of variables for the classification of forest tree species groups is a key factor in improving classification accuracy. A useful feature for forest tree species classification is the phenological parameter. Phenology encompasses very visible processes such as leaf color change in autumn due to leaf senescence (mainly related to the decomposition of chlorophyll) and the gradual greening of leaves and flowering events in spring. Since phenology varies among forest tree species, it is possible to relate forest tree species classification to the phenological cycles of the studied species.

[0081] In this embodiment, the parameters of the MCCDC model and the phenological period parameters (SOS, EOS and LOS) are used for classification and mapping of forest types. Considering the MCCDC model parameters and phenological parameters (SOS, EOS and LOS), a reliable identification framework for forest types (tree species groups) is proposed. The principle of proposing this framework is that based on the established MCCDC parameter discovery, different tree species groups have different model coefficients, and thus respond in different change curves, which provides another unique opportunity for forest tree species classification that is different from traditional spectral and texture feature methods. Figure 7 As shown, the spectral curves of pine in the NIR band are significantly different from those of oak and walnut. The highest SR value of pine is generally lower than 3000, but the highest values ​​of oak and walnut are generally greater than 3000, while the difference between oak and walnut in the NIR band is not much. For SWIR1, the SR value of pine in this band fluctuates between 1600-2100, and the SR value of oak fluctuates between 1000-2000. The peak and trough times of pine and oak are not much different. However, the spectral curve of walnut in the SWIR1 band is obviously opposite to that of the other two categories, which provides a good idea for distinguishing walnut from pine and oak. In addition to these parameters, phenological period can also be used as a good variable to participate in tree species classification, because different tree species often have different phenological periods. This study is based on the tree species distribution map of the second-category forest survey data for two years and the phenological period extracted in Chapter 4 to obtain the average value of each tree species in these two years. As shown Figure 8 As shown in the figure, the SOS of walnut is smaller than that of oak, the SOS of oak is smaller than that of pine, the EOS of pine is later than that of oak and walnut, and finally, the LOS of walnut is longer than that of pine and oak. Therefore, the phenological period is added to the categorical variables together with the MCCDC model parameters.

[0082] For the classification of forest tree species in 2007 and 2016, this embodiment randomly selects regions of interest (ROIs) of training samples of different tree species groups from the second-category forest survey data. It should be noted that according to the small class factor of the second-category forest survey data, the area of ​​pure forest is optimally selected to avoid ROI falling in mixed forest. For the remaining years, due to the lack of second-category forest survey data, with the help of the established annual forest cover map from 1999 to 2019, the area that has maintained a continuous forest state for 21 years is extracted, assuming that the type of forest tree species in the area will not change. The intersection of the continuous forest area and the tree species distribution in 2007 and 2016 is used to find the intersection, and these intersection areas can be used for ROI extraction in the remaining years. As can be seen from Table 2, 350 pine ROIs, 120 oak ROIs, 48 ​​walnut ROIs and 22 other forest tree species ROIs were selected from the second-category survey data in 2007 and 2016, respectively. For the remaining years, 220 pine ROIs, 65 oak ROIs, 26 walnut ROIs and 12 other forest tree species ROIs were selected in turn. These ROIs were combined, 80% of the ROIs were used for classification training, and the remaining ROIs were used for classification validation.

[0083]

[0084]

[0085] Table 2 Number of areas of interest for different forest tree species at different times

[0086] This embodiment performs tree species classification at a more accessible medium spatial resolution scale. In order to obtain more robust classification results, the present invention compares the performance differences of two machine learning algorithms (SVM and RF) in tree species classification at a resolution of 30m. RF and SVM have their own advantages and disadvantages in classification problems, so the present invention attempts to compare the differences between these two classification algorithms in forest tree species group mapping and select the best method to draw the spatiotemporal distribution map of forest tree species groups in the study area.

[0087] Result analysis:

[0088] Classification accuracy of forest tree species groups under different variable combinations and machine learning:

[0089] The present invention first adopts traditional spectral variables, and classifies forest tree species groups based on RF and SVM algorithms respectively (the results are shown in Table 3). The RF classification accuracy results in 2007 show that the highest user accuracy comes from other forest tree species, and the lowest user accuracy comes from pine. The highest producer accuracy appears in non-forest, and the lowest producer accuracy appears in other forest tree species. The overall accuracy of RF classification is 74%, and the Kappa coefficient is 62%. The classification accuracy results of forest tree species groups in 2007 obtained based on the SVM algorithm show that the user accuracy of pine is the lowest, and the highest user accuracy comes from non-forest. The highest producer accuracy comes from non-forest, and the lowest producer accuracy comes from other forest tree species, which is only 28%. The overall accuracy obtained by the SVM algorithm is 80%, and the Kappa coefficient is 72%. For the forest classification accuracy results in 2016, the classification results based on the RF algorithm have improved, and the classification accuracy based on the SVM algorithm has slightly decreased.

[0090]

[0091]

[0092] Table 3 Classification accuracy of forest tree species groups in 2007 and 2016 using RF and SVM models based on spectral variables

[0093] Different from the traditional spectral classification, this embodiment uses the MCCDC coefficient and phenological parameters as input variables, and adds them to the two machine learning algorithms to obtain a new verification accuracy (results are shown in Table 4). Taking the verification results of the forest tree species group in 2007 as an example, for RF classification, the producer accuracy of pine is the highest, which is 92%, and the user accuracy of other forest tree species is the highest, which is 94% (Table 4). On the contrary, the producer accuracy of other forest tree species is the lowest, only 61%. The overall accuracy of all tree species is 86%, and the Kappa coefficient is 81%. For the SVM classification in 2007, the producer accuracy of pine is the highest, which is 92%, and the user accuracy of non-forest is the highest, which is 94%. Similar to the RF classification, the producer accuracy of other forest tree species is also the lowest, which is 74%. The overall accuracy of SVM classification is 89%, and the Kappa coefficient is 83%. For the verification results of forest tree species in 2016, the overall accuracy and Kappa coefficient based on RF classification are slightly increased, and the overall accuracy and Kappa coefficient based on SVM classification are slightly decreased.

[0094]

[0095]

[0096] Table 4 Classification accuracy of forest tree species groups in 2007 and 2016 using RF and SVM models based on MCCDC coefficients and phenological parameters

[0097] Select three small areas to zoom in. Fig. 9 The comparison between the spatial distribution map of RF tree species groups and the spatial distribution map of SVM tree species groups in 2016 is shown. For area (a), the RF model mainly divides the area into pine, but the SVM model divides part of the pine into oak, and there is a small area of ​​other forest tree species at the lower left edge of the pine. For area (b), the oak area obtained by RF is smaller than the oak area obtained by SVM. For area (c), there are a large number of walnuts in this area, and the classification results obtained by RF are basically consistent with those obtained by SVM. Overall, the pine area obtained by the RF model is larger than the pine area obtained by the SVM model, and the oak area obtained by the RF model is smaller than the oak area obtained by the SVM model, that is, the oak area obtained by RF has obvious omissions in space. Therefore, the SVM algorithm based on the MCCDC coefficient and phenological parameters is selected for the classification of forest tree species groups.

[0098] In summary, the present invention uses a classification method based on MCCDC parameters and phenological periods to perform detailed classification of tree species groups at a medium spatial resolution scale. The results show that compared with traditional spectral features and their derivative variables, the classification strategy using MCCDC coefficients and phenological parameters as input variables improves the classification accuracy of forest tree species groups by more than 15%. After determining that MCCDC coefficients and phenological parameters are used as input variables, the SVM and RF algorithms are compared. The classification effect of forest tree species groups by SVM is better than that of RF, and SVM is finally selected as the final tree species group classification result.

[0099] As described above, although the present invention has been shown and described with reference to specific preferred embodiments, it should not be construed as limiting the present invention itself. Various changes in form and details may be made without departing from the spirit and scope of the present invention as defined in the appended claims.

Claims

1. A forest tree species group classification method considering harmonic model coefficients and phenological parameters, Features The steps include: S1. Image acquisition and preprocessing: Obtain multiple sensor satellite images of the area to be classified within a period of time, and use the least squares method to fit the slope and intercept of the acquired images to perform spectral normalization band by band to minimize the spectral response differences between different sensors; S2. Establish a MCCDC model based on multi-source remote sensing data. The specific steps include: S21. Establish a SR time series model and define a pixel-by-pixel reflectance time series model, as shown in formula (1), to express the characteristics of surface reflectance change; （1） Formula (1) contains a Fourier harmonic model and a long-term trend part. is the Julian date; T is the number of days in the cycle, which is 365.25; For the Bands in Julian Date The predicted value of surface reflectivity; there are 8 coefficients in this time series model, coefficient For the The overall mean of the band clear observations, coefficient , , , , , It is used to simulate the annual variation of surface reflectance caused by phenology. Its minimum positive period is one year, and the coefficient Reflects the trend of gradual changes in surface reflectivity caused by factors such as climate change and soil degradation; S22, identifying the mutation point of the time series model, based on the time series model established in step S21, to see whether it satisfies formula (2); （2） In formula (2), is the number of bands; Refers to the number of consecutive observations after the mutation point that needs to be judged; For the Bands in Julian Date Actual observations of surface reflectivity; It refers to the root mean square error between the model-predicted reflectivity value and the actual remote sensing observed reflectivity value in the i-th band; When formula (2) is satisfied, it proves that land cover change has occurred in the area. The changed pixels are sorted in chronological order. First, the first 24 clear pixels at the same location are used to establish an initialized time series model based on formula (1). Then, it is determined whether the difference between the actual observation value and the predicted value of the next 4 pixels in each band exceeds 3 times the root mean square error of the model. If it does not exceed the range, these observations are added to the observations previously involved in the modeling, the parameters are refitted, and the subsequent observations are determined again. If 4 consecutive observations are out of range, the point is defined as the mutation point, that is, the time point when the land cover change occurs. In this process, the sum of the differences between the observation value and the predicted value of all bands needs to be averaged to define the mutation. S23, synthesize daily remote sensing images, use the methods of steps S21 to S22 to generate time series models of different bands for each pixel, and for a single band of each pixel, use the given variable Julian date The SR prediction value of the pixel on a given date is generated based on the time series model. This method is applied to all pixels of each day over a period of time to generate a clear synthetic image with a spatial resolution of 30 m. S3. Extraction of forest annual phenology based on logistic regression equation. The specific method is as follows: S31. Based on all the daily composite images obtained, the annual phenological parameters are extracted year by year. The three vegetation indices EVI, NDVI and LSWI are calculated in turn by equations (3) to (5), and the logistic regression functions are fitted based on the three respectively, so as to determine the estimated values ​​of different phenological parameters. （3） （4） （5） In the formula, , , and They are the surface reflectance values ​​of blue light, red light, near infrared and short-wave infrared 1 bands; S32. Based on the three vegetation indices, the logistic regression function was used to extract the annual SOS, EOS and LOS respectively, and the annual vegetation index time series data were modeled using formula (6); （6） Where t is the observation time, is the vegetation index value EVI, or NDVI, or LSWI, a and b are fitting parameters, the vegetation index value is less than c+d, d is the initial background vegetation index value; After completing the fitting of formula (6), the curvature change rate of the fitting model is calculated by formula (7); （7） In the formula, and are the first and second derivatives of the logistic equation, respectively, and ROC is the rate of change of curvature of the fitted model; Criteria for judging phenological parameters: When the curvature change rate reaches the maximum value for the first time, the corresponding number of days is the vegetation growth start period SOS; when the curvature change rate reaches the minimum value for the last time, the corresponding number of days is the vegetation growth end period EOS; the difference between EOS and SOS is the length of the vegetation growth season LOS; S33, by comparing the SOS and EOS extracted from the three remote sensing indices with the field observed phenological data, the phenological parameters obtained with the LSWI index were selected as the input variables of the recognition algorithm; S4. Considering the phenological parameters and the coefficients of the fitted MCCDC model, the tree species groups were identified by using RF or SVM recognition algorithms to classify the forest tree species groups respectively. The parameters of the MCCDC model and the phenological period parameters were added as input variables to the two machine learning algorithms to draw the spatiotemporal distribution map of the forest tree species groups.

2. The forest tree species group classification method considering harmonic model coefficients and phenological parameters according to claim 1, It is characterized in that In step S4, the SVM recognition algorithm is selected to classify forest tree species groups.

Images