A rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting

Abstract

The present application relates to a kind of based on two-way dynamic harmonic fitting rubber forest age remote sensing estimation method, comprising: S1, based on two-way dynamic harmonic fitting rubber forest change detection;S2, on the basis of the detection result of step S1, data missing pixel is interpolated, and rubber forest planting year graph is obtained;S3, in combination with rubber forest planting year graph and rubber forest spatial distribution graph, rubber forest age calculation is carried out.The method is conducive to improving the accuracy of rubber forest age estimation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of remote sensing, and particularly relates to a rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting. BACKGROUND

[0002] In the prior art, there are two methods for estimating the age of a rubber forest using remote sensing data. One method is to classify rubber forests year by year, and then use the classification results of each year to estimate the age of the rubber forest. For details, please refer to the document "Obtaining rubber plantation age information from very dense Landsat TM & ETM+ time series data and pixel-based image compositing". The other method is to use a change detection algorithm to detect the planting year of the rubber forest by breakpoint detection to estimate the age of the rubber forest. The mainstream change detection algorithms currently used are CCDC and Landtrendr. For details, please refer to the documents "30m Map of Young Forest Age in China" and "A global map of planting years of plantations". The method of classifying rubber forests year by year has errors in each year's classification. Using the classification results of multiple years to estimate the age of the rubber forest can easily cause cumulative errors, resulting in a large error in the final estimation result. Compared with the above method, the change detection method for estimating the age of the rubber forest has improved accuracy. However, the Landtrendr algorithm generates only one value for each pixel per year, which can easily overlook changes within a year. In contrast to Landtrendr, CCDC can use all available data and use dynamic RMSE for breakpoint detection, making it more suitable for estimating the age of rubber forests. In the past, the CCDC algorithm has been used to estimate the age of rubber forests by adjusting the model parameters to improve accuracy. This method is relatively simple and effective, but the adjustment of the parameters has the problems of difficulty in controlling the regulation scale and unclear interaction mechanism between the parameters. Therefore, how to more effectively improve the performance of the change detection algorithm in age extraction and achieve more accurate age estimation is a problem that needs to be solved. SUMMARY

[0003] The present application relates to the technical field of remote sensing, and particularly relates to a rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting.

[0004] To achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows: a rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting, comprising:

[0005] S1, performing change detection on the rubber forest based on bidirectional dynamic harmonic fitting;

[0006] S2, on the basis of the detection result in step S1, the missing data pixels are interpolated to obtain a rubber plantation year graph;

[0007] S3, combining the rubber plantation year graph and the rubber forest spatial distribution graph, the rubber forest age is calculated.

[0008] Further, the implementation method of step S1 is: respectively performing forward stacking and reverse stacking on the long-time sequence growth season Landsat image data removing the rubber forest phenological phase, respectively performing forward harmonic fitting and reverse harmonic fitting according to time to realize forward detection and reverse detection, and detecting the starting time of the nearest stable section to the detection target year, i.e. the rubber plantation time, pixel by pixel; respectively calculating the RMSE of forward detection and reverse detection, and taking the detection result with smaller RMSE as the detection result of the pixel.

[0009] Further, a plurality of detection bands are used to respectively perform forward harmonic fitting and reverse harmonic fitting; for each forward or reverse harmonic fitting of a detection band, the absolute value of the difference between all observation values in each year and the corresponding harmonic fitting prediction value is averaged, and then three-fold normalization processing is performed to obtain the judgment reference value of each year;

[0010] For forward harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of forward harmonic fitting; then the nearest breakpoint to the target year among all the breakpoints obtained by forward harmonic fitting is determined as the rubber plantation time obtained by forward detection, and the detection bands corresponding to the k1 judgment reference values greater than 1 in the rubber plantation time of forward detection are recorded;

[0011] For reverse harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of reverse harmonic fitting; then the nearest breakpoint to the target year among all the breakpoints obtained by reverse harmonic fitting is determined as the rubber plantation time obtained by reverse detection, and the detection bands corresponding to the k2 judgment reference values greater than 1 in the rubber plantation time of reverse detection are recorded;

[0012] It is judged whether there are the same detection bands in the k1 detection bands and the k2 detection bands, if there are the same detection bands, these same detection bands are taken, and then the RMSE of each detection band from the nearest stable section of the rubber plantation time detected by the detection band is calculated, and the rubber plantation time corresponding to the detection band with the smallest RMSE is taken as the detection result; if there are no same detection bands, the RMSE of each detection band from the nearest stable section of the rubber plantation time detected by the detection band is directly calculated, and the rubber plantation time corresponding to the detection band with the smallest RMSE is taken as the detection result.

[0013] Further, the plurality of detection wavebands includes 8 detection wavebands: GREEN, RED, NIR, SWIR1, SWIR2, NDVI, EVI and NBR; k is 3.

[0014] Further, the calculation method of the harmonic fitting prediction value is as follows:

[0015]

[0016] In the formula, denotes the harmonic fitting prediction value of the i th Landsat waveband on the Julian date x; x denotes the Julian date, i denotes the Landsat waveband, T denotes the number of days per year, N denotes the number of years of Landsat, a 0,i denotes the total system coefficient of the i th Landsat waveband, a 1,i ,b 1,i denotes the intra-annual variation coefficient of the i th Landsat waveband, a 2,i ,b 2,i denotes the inter-annual variation coefficient of the i th Landsat waveband.

[0017] Further, in the step S2, for the pixels with less available Landsat data, resulting in that the number of observation values does not reach the harmonic fitting threshold, the change detection is performed using the annual available Landsat data, the data of the pixels with missing data is interpolated, and then the detection result is rounded to obtain the rubber plantation year map.

[0018] Further, the change detection using the annual available Landsat data is specifically forward detection.

[0019] Further, in the step S3, the rubber plantation year map obtained is masked using the target year rubber forest spatial distribution map obtained by classification, and then rubber forest age calculation is performed to finally obtain the rubber forest age spatial distribution map.

[0020] Further, the rubber forest age is calculated by the following formula:

[0021] Age rubber = Year target -Year planting +1

[0022] In the formula, Age rubber denotes the rubber forest age, Year target denotes the target year, and Year planting denotes the rubber plantation year.

[0023] Compared with the prior art, the present application has the following beneficial effects: the present application provides a rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting, overcomes the problems of false detection and missed detection in one-way dynamic harmonic fitting change detection, improves the accuracy of rubber forest age estimation, and has strong practicability and broad application prospect. BRIEF DESCRIPTION OF DRAWINGS

[0024] Figure 1 is a method implementation flowchart of an embodiment of the present application;

[0025] Figure 2 is a schematic diagram of forward harmonic fitting using the EVI detection band in an embodiment of the present application;

[0026] Figure 3 is a schematic diagram of reverse harmonic fitting using the EVI detection band in an embodiment of the present application. DETAILED DESCRIPTION

[0027] The present application will be further described below in conjunction with the drawings and embodiments.

[0028] It should be noted that the following detailed description is exemplary and is intended to provide further explanation of the present application. Unless otherwise indicated, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the present application pertains.

[0029] It should be noted that the terms used herein are only for the purpose of describing specific embodiments, and are not intended to limit the exemplary embodiments according to the present application. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form, and in addition, it should be understood that when the terms "comprise" and / or "include" are used in the specification, there is a feature, step, operation, device, component and / or combination thereof.

[0030] As shown in Figure 1 The present embodiment provides a rubber forest age remote sensing estimation method based on bidirectional dynamic harmonic fitting, which comprises:

[0031] S1, rubber forest change detection based on bidirectional dynamic harmonic fitting.

[0032] S2, based on the detection results of step S1, the missing data pixels are interpolated to obtain a rubber forest planting year map.

[0033] S3, combining the rubber forest planting year map and the rubber forest spatial distribution map, rubber forest age calculation is performed.

[0034] The data sources used in the present embodiment include:

[0035] The long time series Landsat5 / 7 / 8 Level2, Collection2, Tier1 Surface Reflectance dataset is used as the data source. The images are pre-processed, including shadow, cloud, cloud shadow and snow mask, to obtain higher quality images.

[0036] In reality, because the felling and planting of rubber forests is an event that has already occurred, its time is determined, and from long time series images, there will be no deviation whether observing forward or backward. When detecting rubber forests forward, it is a normal growth cycle, and when detecting backward, it is a reverse growth process. This means that when detecting forward, missing or detecting errors of breakpoints have a chance to be detected correctly when detecting backward. Ideally, when we detect forward, the time breakpoint that changes occurs is not detected (i.e. missed), and then the reverse detection can detect the breakpoint; or when detecting forward detects an error breakpoint (i.e. false detection), the reverse detection can detect a breakpoint that is relatively closer to the actual disturbance time. Of course, this optimization is based on the premise that the reverse detection is completely correct. In fact, like forward detection, reverse detection also has false detection and missed detection, so when using the bidirectional dynamic harmonic fitting algorithm for rubber forest change detection, the results of forward and backward detection need to be selected, and the optimal detection result is used as the final detection result.

[0037] Therefore, in step S1 of the embodiment, the long time series growth season (May-October) Landsat image data of the rubber forest phenology are respectively forward stacked and backward stacked, and forward harmonic fitting and backward harmonic fitting are respectively performed according to time to realize forward detection and backward detection, and the starting time of the nearest stable section to the detection target year is detected pixel by pixel, i.e. the rubber forest planting time; the RMSE of forward detection and backward detection is calculated, and the detection result with smaller RMSE is taken as the detection result of the pixel.

[0038] In order to further improve the accuracy of the detection result, the method uses multiple detection bands to perform forward harmonic fitting and backward harmonic fitting respectively. The specific implementation method is as follows.

[0039] 1) For each forward or backward harmonic fitting of a detection band, the absolute value of the difference between all observed values in each year and the corresponding harmonic fitting prediction value is averaged, and then three normalization processing (divided by 3 times of RMSE) is performed to obtain the judgment reference value of each year.

[0040] In the embodiment,

[0041] 2) For forward harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of forward harmonic fitting; then the breakpoint closest to the target year among all the breakpoints obtained by forward harmonic fitting is determined as the rubber plantation time detected by forward detection, and the k1 detection bands with the judgment reference value greater than 1 corresponding to the rubber plantation time detected by forward detection are recorded.

[0042] For reverse harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of reverse harmonic fitting; then the breakpoint closest to the target year among all the breakpoints obtained by reverse harmonic fitting is determined as the rubber plantation time detected by reverse detection, and the k2 detection bands with the judgment reference value greater than 1 corresponding to the rubber plantation time detected by reverse detection are recorded.

[0043] 3) Determine whether there are the same detection bands in the k1 detection bands and the k2 detection bands, if there are the same detection bands, take these same detection bands, and then calculate the RMSE of each detection band from the stable section closest to the rubber plantation time detected by the detection band, and take the rubber plantation time corresponding to the detection band with the minimum RMSE as the detection result; if there are no the same detection bands, directly calculate the RMSE of each detection band from the stable section closest to the rubber plantation time detected by the detection band, and take the rubber plantation time corresponding to the detection band with the minimum RMSE as the detection result.

[0044] In the embodiment, the plurality of detection bands includes 8 detection bands: GREEN, RED, NIR, SWIR1, SWIR2, NDVI, EVI and NBR. k is 3.

[0045] Taking the EVI detection band as an example, the forward harmonic fitting is as shown in Figure 2 , and the reverse harmonic fitting is as shown in Figure 3 .

[0046] The calculation method of the harmonic fitting prediction value is:

[0047]

[0048] In the formula, represents the harmonic (RIRLS) fitting prediction value of the i th Landsat band at Julian date x; x represents the Julian date, i represents the Landsat band, T represents the number of days per year (T = 365), N represents the number of Landsat years, a 0,i represents the total system constant of the i th Landsat band, a 1,i ,b 1,i represents the annual variation coefficient of the i th Landsat band, a2,i 2,i denotes the interannual variability coefficient of the i-th Landsat band.

[0049] The difference between forward detection and reverse detection lies in the processing of input data. The above formula only detects the intra-annual and inter-annual variation of time series in the order of input data by harmonic fitting. If the input data is forward stacked from 1988->2022, the data processing order is 1988->2022, which is forward detection; if the input data is reverse stacked from 2022->1988, the data processing order is 2022->1988, which is reverse detection.

[0050] Due to the problems of optical image such as cloud interference or no pixel, the growth season (May-October) of rubber forest phenology is in the rainy season, at this time, the available Landsat data is less, which will lead to the observation value of part of the pixel cannot reach the threshold, and then the harmonic fitting cannot be carried out, and a large number of pixels appear missing detection.

[0051] Therefore, in step S2 of the embodiment, for the pixels whose observation value number cannot reach the harmonic fitting threshold due to less available Landsat data, forward detection is carried out using the available Landsat data of the whole year, data interpolation is carried out for the missing data pixels, and then the detection result is rounded to obtain the rubber forest planting year graph. In the embodiment, the rounding method is downward rounding (i.e. directly removing the decimal part), because the detection result is the planting year, for the decimal year, no matter whether the rubber tree planting year is 2000.1 year or 2000.9 year, we think it is planted in 2000.

[0052] In step S3 of the embodiment, the target year rubber forest spatial distribution graph obtained by classification is used to mask the obtained rubber forest planting year graph, and then the following formula is used to calculate the rubber forest age:

[0053] Age rubber =Year target -Year planting +1

[0054] In the formula, Age rubber represents the rubber forest age, Year target represents the target year, and Year planting represents the rubber forest planting year.

[0055] Finally, the rubber forest age spatial distribution graph is obtained.

[0056] ​Those skilled in the art will appreciate that embodiments of the application can be readily used as software, hardware, or a combination of software and hardware. In one

[0057] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing each of the flowchart blocks or the functions indicated in the blocks

[0058] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instructions which implement the flowchart Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing each of the flowchart blocks or the functions indicated in the blocks

[0059] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing each of the flowchart blocks or the functions indicated in the blocks

[0060] The above descriptions are only preferred embodiments of the present application and are not intended in any way to limit the present application and other forms. Any skilled in the art can use the disclosed technology to make changes or modifications as equivalent embodiments. However, any simple modification, equivalent change and modification of the above embodiments without departing from the technical solution of the present application, according to the technical essence of the present application, still belongs to the protection scope of the technical solution of the present application.

Claims

1. A method for estimating the age of rubber forest based on bidirectional dynamic harmonic fitting, characterized in that, The method comprises the following steps: S1, detecting the change of rubber forest based on bidirectional dynamic harmonic fitting; S2, based on the detection result of step S1, interpolating the missing data pixels to obtain a rubber forest planting year graph; S3, combining the rubber forest planting year graph and the rubber forest spatial distribution graph to calculate the age of the rubber forest; The implementation method of step S1 is as follows: respectively performing forward stacking and reverse stacking on long-time sequence growth season Landsat image data excluding the rubber forest phenological period, respectively performing forward harmonic fitting and reverse harmonic fitting according to time to realize forward detection and reverse detection, and detecting, pixel by pixel, the starting time of the nearest stable section to the detection target year, i.e. the rubber forest planting time; respectively calculating the RMSE of forward detection and reverse detection, and taking the detection result with smaller RMSE as the detection result of the pixel; A plurality of detection bands are used to respectively perform forward harmonic fitting and reverse harmonic fitting; for each forward or reverse harmonic fitting of a detection band, the absolute value of the difference between all observation values in each year and the corresponding harmonic fitting prediction value is averaged, and then three normalization processing is performed to obtain the judgment reference value of each year; For forward harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of forward harmonic fitting; then the nearest breakpoint to the target year among all the breakpoints obtained by forward harmonic fitting is determined as the rubber forest planting time obtained by forward detection, and the k1 detection bands with the judgment reference value greater than 1 corresponding to the rubber forest planting time of forward detection are recorded; For reverse harmonic fitting, if the judgment reference value of more than k detection bands is greater than 1 in a certain year, it is judged that the year is a breakpoint of reverse harmonic fitting; then the nearest breakpoint to the target year among all the breakpoints obtained by reverse harmonic fitting is determined as the rubber forest planting time obtained by reverse detection, and the k2 detection bands with the judgment reference value greater than 1 corresponding to the rubber forest planting time of reverse detection are recorded; It is judged whether there are the same detection bands in the k1 detection bands and the k2 detection bands, if there are the same detection bands, the same detection bands are taken, and then the RMSE of each detection band from the nearest stable section to the rubber forest planting time detected by the detection band is calculated, and the rubber forest planting time corresponding to the detection band with the minimum RMSE is taken as the detection result; if there are no same detection bands, the RMSE of each detection band from the nearest stable section to the rubber forest planting time detected by the detection band is calculated, and the rubber forest planting time corresponding to the detection band with the minimum RMSE is taken as the detection result.

2. The method according to claim 1, wherein, The plurality of detection bands include 8 detection bands: GREEN, RED, NIR, SWIR1, SWIR2, NDVI, EVI and NBR; k is 3.

3. The method according to claim 1, wherein, The calculation method of the harmonic fitting prediction value is as follows: wherein, represents the harmonic fitted prediction value of the i-th Landsat band on Julian date x; x represents the Julian date, i represents the Landsat band, T represents the number of days per year, N represents the number of years of Landsat, a 0,i represents the total system coefficient of the i-th Landsat band, a 1,i and b 1,i represents the intra-annual variation coefficient of the i-th Landsat band, a 2,i and b 2,i represents the inter-annual variation coefficient of the i-th Landsat band.

4. The method according to claim 1, wherein, In the step S2, for the pixels whose observation value number does not reach the harmonic fitting threshold due to less Landsat data, change detection is performed using the whole year available Landsat data, data interpolation is performed for the pixels with missing data, then the detection result is rounded to obtain the rubber plantation year map.

5. The method according to claim 4, wherein, The change detection using the whole year available Landsat data is specifically forward detection.

6. The method according to claim 1, wherein, In the step S3, the rubber plantation year map obtained is masked using the target year rubber forest spatial distribution map obtained by classification, then rubber forest age calculation is performed, and finally a rubber forest age spatial distribution map is obtained.

7. The method according to claim 6, wherein, The rubber forest age is calculated by using the following formula: Age rubber = Year target -Year planting +1 In the formula, Age rubber represents the age of the rubber forest, Year target represents the target year, Year planting represents the planting year of the rubber forest.

Images