Remote sensing image space-time fusion method for complex mountainous area

Abstract

The invention relates to the technical field of remote sensing, in particular to a complex mountainous area-oriented remote sensing image space-time fusion method, which comprises the following steps of: S1, acquiring Landsat and MODIS time sequence data, topographic data and land cover product data, and preprocessing the data; s2, extracting a slope unit: forming a coupling graph based on the data acquired in the step S1, and segmenting the coupling graph to generate the slope unit; and S3, on the basis of introducing the slope unit, simulating reflectivity change caused by a terrain effect in combination with a terrain correction model, constructing reflectivity conversion coefficients of images with different spatial resolutions under a complex terrain, and finally, calculating and generating a fused image according to the weight matrix and the conversion coefficients. The method not only has excellent performance in the aspects of spectrum precision and spatial information retention, but also shows higher flexibility and robustness when processing a complex terrain region, and overcomes the limitation of the ESTARFM in the complex terrain region.

Description

Technical Field

[0001] This invention relates to the field of remote sensing technology, and in particular to a method for spatiotemporal fusion of remote sensing images for complex mountainous areas. Background Technology

[0002] Long-term remote sensing imagery is crucial for monitoring the dynamics of heterogeneous landscapes, particularly in areas such as ecosystem monitoring, land use / land cover change detection, and vegetation / crop monitoring. However, a trade-off always exists between the spatial and temporal resolution of remote sensing imagery. Landsat series imagery, with a spatial resolution of 30 m, is often used for long-term land cover and vegetation monitoring, but its 16-day revisit period makes it difficult to capture short-term surface dynamics. Satellite sensors such as MODIS and AVHRR have shorter revisit periods, allowing for multiple coverages of the Earth within a short timeframe to capture rapid surface changes; however, their insufficient spatial resolution makes them unsuitable for monitoring highly heterogeneous regions. Furthermore, acquiring high spatiotemporal resolution remote sensing imagery is extremely challenging due to the need to balance spatial resolution and revisit period, as well as the challenges posed by cloud contamination and sensor malfunctions.

[0003] Spatiotemporal fusion based on multi-source remote sensing data is an effective way to solve the difficulty of acquiring high spatiotemporal resolution remote sensing images. It combines remote sensing data acquired by multiple sensors on different dates with different spatial resolutions to generate remote sensing images at new time points. Among remote sensing spatiotemporal fusion methods, the spatiotemporal adaptive reflectance fusion model is the most classic. Among various spatiotemporal fusion methods, ESTARFM has high fusion accuracy in heterogeneous regions and has been applied on cloud platforms, enabling rapid online fusion over large areas. However, ESTARFM's performance is limited in areas with complex terrain. ESTARFM uses a spectral thresholding method for similar pixel selection, which relies on the similarity between the spectra of ground objects. Topographic effects can easily lead to mismatches of similar pixels, resulting in uncertainty. Especially in mountainous areas with complex terrain, the shaded and illuminated sides may lead to phenomena such as "same object, different spectrum" or "different object, same spectrum," affecting subsequent calculations. Meanwhile, considering that mountainous areas account for approximately 26.4% of the global land area and possess rich ecosystems, species, and genetic diversity, the acquisition of high spatiotemporal resolution remote sensing imagery in mountainous regions is challenging due to their complex vertical climate and strong spatial heterogeneity, necessitating the spatiotemporal fusion of remote sensing data. Therefore, improving the fusion accuracy of ESTARFM in complex terrain regions is of great significance, but this aspect is rarely addressed in existing technologies and urgently needs improvement. Summary of the Invention

[0004] The features and advantages of the present invention are set forth in part in the description which follows, or may be apparent from the description, or may be learned by practicing the invention.

[0005] To overcome the problems of existing technologies, this invention provides a spatiotemporal fusion method for remote sensing images in complex mountainous areas, specifically including the following steps: S1. Collect Landsat and MODIS time series data, terrain data, and land cover product data, and perform data preprocessing. S2. Extracting slope elements: Based on the data collected in step S1, a coupled graph is formed, and the coupled graph is segmented to generate slope elements; S3. Based on the introduction of the slope unit, the reflectance change caused by the terrain effect is simulated by the terrain correction model, and the reflectance conversion coefficients of images with different spatial resolutions under complex terrain are constructed. Finally, the fused image is generated according to the weight matrix and the conversion coefficients.

[0006] Preferably, the specific operation of step S2 includes: extracting terrain parameters from the DEM, generating a terrain parameter coupling map, and then using simple non-iterative clustering provided by GEE to perform image segmentation on the terrain parameter coupling map to generate various homogeneous and heterogeneous segments, i.e., slope units.

[0007] Preferably, the homogeneity within the slope unit is evaluated by the local variance V, and the heterogeneity between the slope units is evaluated by the autocorrelation index I, wherein the calculation formulas for the local variance V and the autocorrelation index I are shown in formula (1) and formula (2): Formula (1), Formula (2), In the formula, For all in a given partition The first slope unit indivual, It is the first Surface area of ​​each slope unit It is the first The circular square of the slope aspect of each slope unit It is a spatial adjacency matrix. It is the slope direction of the slope unit.

[0008] Preferably, the optimal cutting parameters for slope unit extraction are determined based on the scale of the moving window. Specifically, the operation includes: randomly selecting sample points within a typical study area according to the minimum sample size formula, generating an inscribed rectangle that serves as the circumcircle of the moving window, and calculating the slope units extracted at different superpixel seed position intervals within this circumcircle. To determine the optimal cutting parameters and extract slope elements accordingly; where, The calculation formula is shown in formula (3): Formula (3),

[0009] In the formula, The upper 5th quantile of V. The lower 5th quantile of V. The upper 5th percentile of I, It is the lower 5th percentile of I.

[0010] Preferably, the specific operation steps of step S3 include: S301. Similar pixels are selected based on slope units, and a weight matrix is ​​calculated based on the selected similar pixels. S302. Construct reflectance conversion coefficients for images with different spatial resolutions under complex terrain by combining terrain correction models; S303. Generate the fused image based on the weight matrix and transformation coefficients.

[0011] Preferably, the calculation formula for the similar pixel screening is as shown in formula (4): Formula (4), In the formula, represents the area within the moving window. Location cell and center cell The difference in reflectivity, Spectral threshold, The slope unit where the central pixel is located.

[0012] Preferably, step S302 further includes: known High spatial resolution remote sensing imagery at any given time, The formulas for calculating the reflectance of a flat surface in a high spatial resolution remote sensing image are shown in formulas (5) and (6): Formula (5), Formula (6), In the formula, and for and Time prediction The reflectance of a flat surface in high spatial resolution remote sensing imagery at any given time; for Reflectivity of a flat surface at all times for Reflectivity of a flat surface at all times; When the changes in surface reflectance due to topographic effects are taken into account, then The surface reflectance of the high spatial resolution remote sensing image at a given time is calculated using formulas (7) and (8): Formula (7), Formula (8), In the formula, for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. These are empirical parameters calculated based on the linear relationship between reflectance and solar altitude angle in actual images; for At any given moment, the zenith angle of the sun. The slope.

[0013] Preferably, step S303 further includes: placing Time and prediction High spatial resolution remote sensing images and A regression model is established using low spatial resolution images at different times. The transformation coefficients are calculated by making full use of the information of adjacent similar pixels. The formulas for calculating the transformation coefficients are shown in formulas (9) and (10). Formula (9), Formula (10), In the formula, and for Reflectivity of similar pixels with high spatial resolution at any given time. and Based on Time-of-flight simulation Reflectivity of similar pixels with high spatial resolution at any given time. , , for Reflectance of similar pixels with low spatial resolution at any given time The coefficients of the regression model established for similar pixels in high spatial resolution images and similar pixels in low spatial resolution images. , The conversion coefficient is the center pixel.

[0014] Preferably, step S303 further includes: based on the reflectance conversion coefficients of images with different spatial resolutions obtained in step S302, then based on... Time prediction The formulas for calculating Landsat at time t are shown in formulas (11) and (12): Formula (11), Formula (12), In the formula, for Landsat pixel reflectance at time 1 The weight matrix calculated for the selected similar pixels. This is the transformation coefficient matrix simulated under complex terrain. for arrive The reflectance difference of MODIS pixels at time t. for arrive The error caused by the SCS+C model at any time.

[0015] Preferably, the data preprocessing includes atmospheric correction, geometric correction, cloud removal, and cloud shadow removal.

[0016] The beneficial effects of this invention are as follows: Based on the introduction of slope units, this invention optimizes the screening strategy for similar pixels and combines a terrain correction model to simulate reflectance changes caused by terrain effects. It constructs reflectance conversion coefficients for images with different spatial resolutions under complex terrain conditions to improve the applicability of ESTARFM in complex terrain areas. By introducing slope units and a terrain correction model, the limitations of traditional methods in complex terrain are overcome, resulting in higher homogeneity of pixels within the same slope unit. This effectively avoids the problem of mismatched similar pixels caused by excessive spectral differences due to terrain changes. The similar pixel screening results in this invention are superior to those of the traditional ESTARFM model, with a significant decrease in the standard deviation of reflectance for the screened similar pixels. The conversion coefficients in this invention consider the influence of terrain effects on reflectance changes, making the conversion coefficients more consistent with the actual surface reflectance variation patterns. This invention can more accurately identify similar pixels and obtain conversion coefficients that better conform to the actual surface reflectance variation patterns. It not only performs excellently in terms of spectral accuracy and spatial information preservation but also demonstrates greater flexibility and robustness in handling complex terrain areas, overcoming the limitations of ESTARFM in complex terrain regions. Attached Figure Description

[0017] The present invention will be described in detail below with reference to the accompanying drawings and examples. The advantages and implementation methods of the present invention will become more apparent from this description. The accompanying drawings are for illustrative purposes only and do not constitute any limitation on the present invention. In the accompanying drawings: Figure 1 This is a flowchart illustrating a method for spatiotemporal fusion of remote sensing images in complex mountainous areas, as described in a specific embodiment of the present invention. Figure 2 This is a schematic diagram of the process for extracting slope units in a specific embodiment of the present invention; Figure 3 This is a schematic diagram of similar pixel screening in a specific embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the change in surface reflectance caused by topographic effects in a specific embodiment of the present invention; Figure 5 This is a visual comparison of the fused images of the present invention and ESTARFM in a typical region A in a specific embodiment of the present invention, where a is a magnified view of point a in typical region A, and b is a magnified view of point b in typical region A. Figure 6 The image shown here is a visual comparison of the fused images of the present invention and ESTARFM in a typical region B in a specific embodiment of the present invention. Here, a is a magnified view of point a in the typical region B, and b is a magnified view of point b in the typical region B. Figure 7 The APA diagram shows the performance of the present invention and the ESTARFM algorithm in typical region A of a specific embodiment of the present invention. Figure 8 The APA diagram shows the performance of the present invention and the ESTARFM algorithm in typical region B of a specific embodiment of the present invention. Detailed Implementation

[0018] The embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and should not be construed as limiting the scope of the invention.

[0019] like Figure 1 As shown, this invention provides a spatiotemporal fusion method for remote sensing images in complex mountainous areas, specifically including the following steps: S1. Collect Landsat and MODIS time series data, terrain data, and land cover product data, and preprocess the data. The high spatial resolution imagery uses Landsat time-series data, specifically Landsat-5 (1984–2012), Landsat-7 (1999–2021), and Landsat-8 (2013–Present) Surface Reflectance Tier 1 products from the Landsat Collection 2 data released by the U.S. Geological Survey (USGS) in 2021. This data has a spatial resolution of 30 m and a revisit period of 16 days. The data underwent geometric and atmospheric corrections. Landsat 5 and Landsat 7 products underwent atmospheric correction using the Landsat Ecosystem Disturbance Adaptive Processing System algorithm (LEDAPS) software, while Landsat 8 underwent correction using the Land Surface Reflectance Code (LaSRC). Quality Assessment (QA) bands generated using the FMASK method were used to remove pixels affected by clouds, cloud shadows, etc. These preprocessing operations were performed using the Google Earth Engine cloud platform.

[0020] The low spatial resolution imagery uses MODIS time series data, provided by NASALP DAAC of the U.S. Geological Survey's EROS Center. This includes 500m surface reflectance (MOD09A1.061 Terra SurfaceReflectance 8-Day Global 500m), 1000m surface temperature (MOD11A2.061 Terra LandSurface Temperature and Emissivity 8-Day Global 1km), and 500m vegetation index (MOD13A1.061 Terra Vegetation Indices 16-Day Global 500m). The product itself has already undergone atmospheric and geometric correction; only cloud and cloud shadow removal is needed.

[0021] After completing the Landsat and MODIS preprocessing, the MODIS image is reprojected and resampled to the resolution and range corresponding to the Landsat image. The resampling and registration functions provided by GEE are used, with bilinear interpolation as the resampling method.

[0022] The terrain data used was the SRTM (Shuttle Radar Topography Mission) digital elevation model (DEM) with a sampling interval of 1 arcsecond in the Earth's conformal coordinate system. The SRTM V3 product (SRTM Plus) was provided by NASA JPL, with a spatial resolution of 30 m. The data was sourced from https: / / gdex.cr.usgs.gov / gdex / . The SRTM DEM was primarily used to extract slope elements.

[0023] The land cover product used was GLC_FCS30D (data source: https: / / zenodo.org / records / 8239305), with a resolution of 30 m and an overall accuracy of 80.88% (±0.27%). The GLC_FCS30D product was reclassified into cultivated land, forest land, grassland, water bodies, construction land, and unused land.

[0024] S2. Extracting slope elements: Based on the data collected in step S1, a coupled graph is formed, and the coupled graph is segmented to generate slope elements; the specific operations are as follows: like Figure 2 As shown, topographic parameters, including slope, aspect, and topographic location index (TPI), are extracted from the DEM to generate a topographic parameter coupling map. This coupling map can effectively visualize and assess the slope variation from the ridgeline to the downhill valley near the drainage network, reflecting the contours of the slope units. Then, Simple Non-Iterative Clustering (SNIC) provided by GEE is used to segment the topographic parameter coupling map, generating various homogeneous and heterogeneous segments, i.e., slope units.

[0025] In the above-mentioned ramp unit extraction process, the parameter affecting the SNIC segmentation effect is mainly the superpixel seed position spacing (initial seed points uniformly distributed in the image, each seed point will become the center of a superpixel, denoted as "SIZE"). Specifically, SIZE determines the fineness of the ramp unit segmentation; the smaller the SIZE, the more initial seed points in the image, and the finer the ramp unit segmentation. In order to quickly evaluate the segmentation quality of ramp units of different regions and different SIZEs, this application uses local variance V to evaluate the homogeneity within ramp units and autocorrelation index I to evaluate the heterogeneity between ramp units. The calculation formulas for V and I are shown in Formula 1 and Formula 2: Formula 1, Formula 2, In the formula, For all in a given partition The first slope unit indivual, It is the first The surface area of ​​each slope unit. It is the first The circular square of the slope aspect of each slope unit. It is a spatial adjacency matrix. It is the slope direction of the slope unit.

[0026] To avoid the effects of smaller V and I values ​​caused by finer or coarser cuts of the slope element, the following approach is adopted: The objective function quantifies the effect of slope element cutting. If... The smaller the value, the better the slope element cutting effect is considered to be. The calculation formula is shown in Formula 3: Formula 3, In the formula, The upper 5th quantile of V. The lower 5th quantile of V. The upper 5th percentile of I, It is the lower 5th percentile of I.

[0027] To ensure that the extracted slope units have scale consistency in subsequent remote sensing spatiotemporal fusion, the selection of parameters during slope unit extraction should be based on the scale of the moving window.

[0028] To determine the optimal cutting parameters for extracting ramp elements within a moving window, sample points were randomly selected from a typical study area based on the minimum sample size formula. An inscribed rectangle was then generated as the circumcircle of the moving window. The slope elements extracted within this circumcircle at different sizes were calculated. To determine the optimal cutting parameters and extract slope elements accordingly.

[0029] The formula for the minimum sample size is shown in Formula 5: Formula 17, In the formula, Confidence level corresponding value, Standard deviation, This is within acceptable error limits.

[0030] Slope units extracted using optimal cutting parameters and those extracted using a GIS-based hydrological model method were visually compared and evaluated. The evaluation method involved overlaying the slope units extracted by both methods onto the original DEM image of the study area. The specific process for extracting slope units using the GIS-based hydrological model method included: First, subtracting the original DEM value from the sum of its maximum and minimum values ​​to obtain a reverse DEM, thus transforming the original valley lines into ridge lines and vice versa. Second, filling was performed on both the original and reverse DEMs. Third, the flow direction was calculated using the eight-direction method, and the cumulative discharge was obtained based on the flow direction. Fourth, the cumulative discharge values ​​exceeding a certain threshold in each grid unit were selected to generate a river network. Finally, the watershed maps of the original and reverse DEMs were obtained, which, after conversion to vectors, became the slope unit outlines.

[0031] S3. Based on the introduction of slope units, and combined with a terrain correction model to simulate reflectance changes caused by terrain effects, reflectance conversion coefficients for images with different spatial resolutions under complex terrain are constructed. Finally, a fused image is generated based on the weight matrix and conversion coefficients, thereby obtaining an improved ESTARFM algorithm suitable for complex terrain, namely CT-ESTARFM, to improve the fusion accuracy of ESTARFM in complex terrain areas. Specifically, the steps are as follows: S301. Similar pixels are selected based on slope units, and a weight matrix is ​​calculated based on the selected similar pixels: pixels within the moving window that have the same land cover type as the central pixel are called "similar pixels." In areas with complex topography, pixels within the same slope unit are considered to have high homogeneity and are approximately of the same land cover type. Therefore, when a pixel satisfies both the spectral threshold and is within the same slope unit as the central pixel, it is considered a similar pixel, such as... Figure 3 As shown, the formula for calculating similar pixels is shown in Formula 4: Formula 4, In the formula, Within the moving window Location cell and center cell The difference in reflectivity, The spectral threshold is given by Equation 4. The slope unit where the central pixel is located.

[0032] Calculate the weights of all similar pixels. Weight It is composed of the spectral differences and spatial similarity of similar pixels, and its calculation formula is as follows: Formula 18: Formula 18, in, Calculated using Formula 19: Formula 19, In the formula, Due to spectral differences, For distance index; where spectral difference Reflectivity of low spatial resolution pixels and high spatial resolution pixel reflectivity The calculation is shown in Formula 20: Formula 20, Distance Index The distance from the center pixel to similar pixels is calculated using the formula shown in Formula 21: Formula 21, S302. Construct reflectance conversion coefficients for images with different spatial resolutions under complex terrain by combining the terrain correction model. For example... Figure 4 As shown, and The endpoint of a certain time interval. Time for and At some point between times, since it is known High spatial resolution remote sensing imagery at any given time, The formulas for calculating the reflectance of a flat surface in a high spatial resolution remote sensing image are shown in Formulas 5 and 6: Formula 5, Formula 6, In the formula, and for and Time prediction The reflectance of a flat surface in high spatial resolution remote sensing imagery at any given time; for Reflectivity of a flat surface at all times for Reflectivity of a flat surface at any given time; because when there is no topographic effect... If the reflectance of a flat surface changes linearly at any given time, then... Calculated using Formula 16: Formula 16, In the formula, for Reflectivity of a flat surface at all times.

[0033] Therefore, when the change in surface reflectance due to topographic effects is taken into account, then The surface reflectance of the high spatial resolution remote sensing image at any given time is: Formula 7, Formula 8, In the formula, for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. These are empirical parameters calculated based on the linear relationship between reflectivity and solar altitude angle in actual images. for At any given moment, the zenith angle of the sun. The slope is where, It can be calculated using formulas 13 and 14: Formula 13, Formula 14, In the formula, The zenith angle of the sun. For slope, Slope direction, The solar altitude angle, This is the solar azimuth angle.

[0034] To avoid overly fragmented patches, a morphological dilation operator is used to expand the slope cell boundaries, combined with stable linear regression for modeling. For slope cells consisting of fewer pixels, the dilation operator can thicken the slope cell boundaries, incorporating adjacent pixels for regression modeling. This avoids situations where too few pixels prevent modeling or outliers significantly impact the data structure, thus altering the entire regression model.

[0035] use Time and prediction High spatial resolution remote sensing images and A regression model is established using low spatial resolution images at different times, and the transformation coefficients are calculated by making full use of the information of adjacent similar pixels. The formulas for calculating the transformation coefficients are shown in Formulas 9 and 10: Formula 9, Formula 10, In the formula, and for Reflectivity of similar pixels with high spatial resolution at any given time. and Based on Time simulation Reflectivity of similar pixels with high spatial resolution at any given time. , , for Reflectance of similar pixels with low spatial resolution at any given time The coefficients of the regression model established for similar pixels in high spatial resolution images and similar pixels in low spatial resolution images. , These are the transformation coefficients for the center pixel. When a model cannot be established... .

[0036] S303. Calculate and generate the fused image based on the weight matrix and transformation coefficients: Based on the reflectance transformation coefficients of the images with different spatial resolutions obtained in step S302, then... Time prediction The formulas for calculating Landsat at time t are shown in Formulas 11 and 12: Formula 11, Formula 12, In the formula, for Landsat pixel reflectance at time 1 The weight matrix calculated for the selected similar pixels. This is the transformation coefficient matrix simulated under complex terrain. for arrive The reflectance difference of MODIS pixels at time t. for arrive The error caused by the SCS+C model at any time.

[0037] Due to the prediction time Between and At time t, the errors in the two predicted directions satisfy: ,but The Landsat pixel reflectance at time t is shown in Equation 17: Formula 17, Finally, a fused image is generated.

[0038] For the constructed fused imagery of a certain river basin from 2001 to 2023, the root mean square error (RMSE), average difference (AD), and similarity coefficient (r) were used to evaluate the fusion results of a remote sensing image spatiotemporal fusion method for complex mountainous areas proposed in this invention. The basin exhibits a complete range of topographical types, including plains, plateaus, hills, gently undulating mountains, moderately undulating mountains, undulating mountains, and extremely undulating mountains. Among these, gently undulating mountains are the dominant terrain, accounting for approximately 59.78%, followed by hills, accounting for approximately 32.44%. Based on the topographical classification results and land cover conditions, two typical areas were selected as the experimental area (typical area A) and the test area (typical area B) to verify the robustness of this invention. Figure 5 As shown, for ease of distinction, the spatiotemporal fusion method for remote sensing images in complex mountainous areas in this invention is abbreviated as CT-ESTARFM. Images generated using this method have a visually closer color tone to real images, accurately predicting the surface reflectance of ridges and adjacent slopes while preserving clear texture details. In contrast, the ESTARFM-fused image exhibits outliers in water bodies and their adjacent mountain areas, overestimates green band reflectance in shady slopes and adjacent sunny slopes, and loses spectral information and spatial details in some shaded areas. Figure 6 As shown, in the typical study area B, the tonal values ​​of the images generated by both CT-ESTARFM and ESTARFM are close to the real images. However, in terms of spatial information representation, ESTARFM overestimates (sharpenes), while CT-ESTARFM underestimates (smooths). Figures 7-8As shown, CT-ESTARFM outperforms ESTARFM. In typical study area A, the performance evaluation result of CT-ESTARFM is "average," close to "good," while the performance evaluation result of ESTARFM is "ineffective." This indicates that the image fused by CT-ESTARFM has better spectral information, while ESTARFM cannot effectively retrieve spectral or spatial information. The average RMSE of the CT-ESTARFM fusion result is 0.019, a reduction of 52.80% compared to ESTARFM. For different bands, CT-ESTARFM performs better in the visible light band, with reduced RMSE of reflectance in the blue, green, and red light bands by 56.77%, 72.85%, and 61.89%, respectively. In terms of representing spatial information, CT-ESTARFM and ESTARFM are similar, with average edges of -0.16 and -0.10, respectively. In typical study area B, except for the near-infrared band, the performance evaluation of CT-ESTARFM is "average," with an average RMSE of 0.016, a reduction of 32.66% compared to ESTARFM. For different wavelengths, CT-ESTARFM still performs better in the visible light band, with reduced RMSE of reflectance in the blue, green, and red bands by 25.64%, 0.32%, and 70.12%, respectively. In characterizing spatial information performance, ESTARFM has an average edge of 0.02, while CT-ESTARFM has an average edge of -0.11. Furthermore, in typical study area A, the AD values ​​of CT-ESTARFM and ESTARFM are close. Specifically, the AD value for shortwave infrared 2 reflectance shows a significant improvement, decreasing from overestimated to between -0.001 and 0.001. The AD values ​​for blue, near-infrared, and shortwave infrared 1 reflectance are underestimated, while the AD values ​​for green and red reflectance are between -0.001 and 0.001. In typical study area B, except for the near-infrared and shortwave infrared 1 reflectance, which are overestimated, the AD values ​​for reflectance in all other bands are underestimated.

[0039] In summary, CT-ESTARFM outperforms the traditional ESTARFM model in similar pixel selection for both typical study areas A and B, with a significant decrease in the standard deviation of reflectance for the selected similar pixels. Regarding the simulation of conversion coefficients, CT-ESTARFM considers the influence of terrain effects on reflectance variations, making the simulated conversion coefficients more consistent with actual surface reflectance changes. The RMSE of the CT-ESTARFM conversion coefficients shows a significant decrease, indicating stronger adaptability to terrain effects. In terms of fusion performance, CT-ESTARFM demonstrates advantages in complex terrain areas. Therefore, CT-ESTARFM exhibits superior fusion accuracy across different bands compared to ESTARFM, while effectively preserving the spatial information of the fused image, demonstrating strong robustness and adaptability.

[0040] The preferred embodiments of the present invention have been described above with reference to the accompanying drawings. Those skilled in the art can implement the present invention in various modifications without departing from its scope and spirit. For example, a feature shown or described in one embodiment can be used in another embodiment to obtain yet another embodiment. The above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. All equivalent changes made based on the description and drawings of the present invention are included within the scope of the present invention.

Claims

1. A method for spatiotemporal fusion of remote sensing images in complex mountainous areas, characterized in that, Specifically, the following steps are included: S1. Collect Landsat and MODIS time series data, terrain data, and land cover product data, and perform data preprocessing. S2. Extracting slope elements: Based on the data collected in step S1, a coupled graph is formed, and the coupled graph is segmented to generate slope elements; S3. Based on the introduction of the slope unit, the reflectance change caused by the terrain effect is simulated by the terrain correction model, and the reflectance conversion coefficients of images with different spatial resolutions under complex terrain are constructed. Finally, the fused image is generated according to the weight matrix and the conversion coefficients.

2. The method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 1, characterized in that, The specific operations of step S2 include: extracting terrain parameters from the DEM, generating a terrain parameter coupling map, and then using simple non-iterative clustering provided by GEE to perform image segmentation on the terrain parameter coupling map to generate various homogeneous and heterogeneous segments, i.e., slope units.

3. The method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 2, characterized in that, The homogeneity within the slope unit is evaluated by the local variance V, and the heterogeneity between the slope units is evaluated by the autocorrelation index I. The calculation formulas for the local variance V and the autocorrelation index I are shown in formula (1) and formula (2): Official (1) Official (2) In the formula, For all in a given partition The first slope unit indivual, It is the first Surface area of ​​each slope unit It is the first The circular square of the slope aspect of each slope unit It is a spatial adjacency matrix. It is the slope direction of the slope unit.

4. The method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 2, characterized in that, The optimal cutting parameters for slope unit extraction are determined based on the scale of the moving window. Specifically, the process includes: randomly selecting sample points within a typical study area according to the minimum sample size formula, generating an inscribed rectangle that serves as the circumcircle of the moving window, and calculating the slope units extracted at different superpixel seed position intervals within this circumcircle. To determine the optimal cutting parameters and extract slope elements accordingly; where, The calculation formula is shown in formula (3): Official (3) In the formula, for The upper 5th percentile, for The lower 5th percentile, for The upper 5th percentile, for The lower 5th percentile.

5. The method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 1, characterized in that, The specific steps of step S3 include: S301. Similar pixels are selected based on slope units, and a weight matrix is ​​calculated based on the selected similar pixels. S302. Construct reflectance conversion coefficients for images with different spatial resolutions under complex terrain by combining terrain correction models; S303. Generate the fused image based on the weight matrix and transformation coefficients.

6. The method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 5, characterized in that, The calculation formula for the similar pixel selection is shown in formula (4): Official (4) In the formula, represents the area within the moving window. Location cell and center cell The difference in reflectivity, Spectral threshold, The slope unit where the central pixel is located.

7. A method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 5, characterized in that, Step S302 further includes: known High spatial resolution remote sensing imagery at any given time, The formulas for calculating the reflectance of a flat surface in a high spatial resolution remote sensing image are shown in formulas (5) and (6): Official (5) Official (6) In the formula, and for and Time prediction The reflectance of a flat surface in high spatial resolution remote sensing imagery at any given time; for Reflectivity of a flat surface at all times for Reflectivity of a flat surface at all times; When the changes in surface reflectance due to topographic effects are taken into account, then The surface reflectance of the high spatial resolution remote sensing image at a given time is calculated using formulas (7) and (8): Official (7) Official (8) In the formula, for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. for Time prediction Surface reflectance of high spatial resolution remote sensing imagery at any given time. These are empirical parameters calculated based on the linear relationship between reflectance and solar altitude angle in actual images; for At any given moment, the zenith angle of the sun. The slope.

8. A method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 7, characterized in that, Step S302 further includes: placing Time and prediction High spatial resolution remote sensing images and A regression model is established using low spatial resolution images at different times. The transformation coefficients are calculated by making full use of the information of adjacent similar pixels. The formulas for calculating the transformation coefficients are shown in formulas (9) and (10). Official (9) Official (10) In the formula, and for Reflectivity of similar pixels with high spatial resolution at any given time. and Based on Time-of-flight simulation Reflectivity of similar pixels with high spatial resolution at any given time. , , for Reflectance of similar pixels with low spatial resolution at any given time The coefficients of the regression model established for similar pixels in high spatial resolution images and similar pixels in low spatial resolution images. , The conversion coefficient is the center pixel.

9. A method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 7, characterized in that, Step S303 further includes: based on the reflectance conversion coefficients of images with different spatial resolutions obtained in step S302, then based on... Time prediction The formulas for calculating Landsat at time t are shown in formulas (11) and (12): Official (11) Official (12) In the formula, for Landsat pixel reflectance at time 1 The weight matrix calculated for the selected similar pixels. This is the transformation coefficient matrix simulated under complex terrain. for arrive The difference in reflectance of MODIS pixels at time t. for arrive The error caused by the SCS+C model at any time.

10. A method for spatiotemporal fusion of remote sensing images for complex mountainous areas according to claim 1, characterized in that, The data preprocessing includes atmospheric correction, geometric correction, cloud removal, and cloud and shadow removal.

Images