Wetland landscape pattern monitoring method considering sub-pixel information

Abstract

The invention discloses a wetland landscape pattern monitoring method considering sub-pixel information, and the method comprises the steps: fusing a first remote sensing image with a high spatial resolution and a second remote sensing image with a high temporal resolution into a fused image, decomposing mixed pixels in the fused image into sub-pixels, and carrying out the decomposition of the sub-pixels based on the decomposition result of the mixed pixels. And on the sub-pixel level, the landscape pattern index, the landscape pattern space dynamic attitude and the landscape pattern time sequence change characteristics are calculated for the fused image, so that the fine granularity of the landscape pattern change of the remote sensing image is accurately monitored.

Description

Technical Field

[0001] The present invention relates to the technical field of humidity landscape evolution monitoring, and in particular to a wetland landscape pattern monitoring method considering sub-pixel information. Background Art

[0002] The use of satellite remote sensing technology to dynamically monitor changes in wetland landscape patterns has been widely used. However, remote sensing images obtained using different technologies usually have different spatial resolutions and temporal resolutions. Using low-resolution remote sensing images to monitor the dynamic situation of changes in wetland landscape patterns, the low-resolution remote sensing images themselves will affect the accuracy of monitoring the dynamic changes in wetland landscape patterns. Therefore, in the prior art, some means are used to first improve the resolution of remote sensing images. However, most of the existing methods for improving the resolution of remote sensing images perform image processing such as denoising on a single type of remote sensing images, which has a poor effect on improving the resolution. Therefore, providing an effective remote sensing image resolution improvement method is the first technical problem that the present invention needs to solve. In addition, after obtaining high-resolution remote sensing images, how to accurately identify areas where the landscape pattern has changed dynamically from the remote sensing images is the second technical problem that the present invention needs to solve. Summary of the invention

[0003] The present invention aims to improve the accuracy of wetland landscape evolution monitoring and provides a wetland landscape pattern monitoring method taking sub-pixel information into consideration.

[0004] To achieve this object, the present invention adopts the following technical solutions:

[0005] A wetland landscape pattern monitoring method considering sub-pixel information is provided, comprising the steps of:

[0006] S1, after data preprocessing, a first remote sensing image with high spatial resolution and a second remote sensing image with high temporal resolution are fused into a fused image with high spatiotemporal resolution;

[0007] S2, decomposing the mixed pixel in the fused image to obtain the proportion and position of the sub-pixel inside the mixed pixel;

[0008] S3, according to the decomposition result of the mixed pixels in the fused image, calculate the landscape pattern index, and / or the spatial dynamic degree of the landscape pattern, and / or the temporal change characteristics of the landscape pattern of the fused image, and complete the monitoring of the landscape pattern change.

[0009] Preferably, the first remote sensing image is a Landsat remote sensing image, and the second remote sensing image is a MODIS remote sensing image.

[0010] Preferably, in step S1, the data preprocessing method for the first remote sensing image and / or the second remote sensing image includes any one or more of radiation correction, geometric correction, atmospheric correction, reprojection, and resampling;

[0011] The method for performing the radiation correction on the first remote sensing image and / or the second remote sensing image is expressed by the following formula (1):

[0012]

[0013] In formula (1), L λ is the solar radiation intensity at the top of the atmosphere;

[0014] Q col is the DN value recorded for the pixel;

[0015] Q colmax is the Q recorded for all pixels in the remote sensing image to be corrected col The maximum value of

[0016] LMAX λ Indicates the upper limit of the remote sensor spectral radiation value;

[0017] LMIN λ Indicates the lower limit of the remote sensor spectral radiation value;

[0018] LMAX λ 、LMIN λ It is a constant value recorded in the remote sensing parameter file of the remote sensor.

[0019] Preferably, the method for geometrically correcting the remote sensing image to be corrected is expressed by the following formulas (2)-(3):

[0020]

[0021] In formulas (2)-(3), X and Y represent the coordinate values ​​of a pixel in the remote sensing image to be corrected on the X-axis and Y-axis respectively;

[0022] F x 、F y Represents the resampling correction distortion function associated with the X-axis and the Y-axis respectively;

[0023] (μ,v) represents the actual geographic coordinates of the pixel to be corrected on the X-axis and Y-axis in the corrected image space;

[0024] μ p v j Represents the product of the horizontal and vertical coordinate orders of the actual geographic coordinates;

[0025] p represents the order of the horizontal axis;

[0026] j represents the order of the ordinate;

[0027] n represents the order of the polynomial;

[0028] a pj and b pj is the coefficient.

[0029] Preferably, the method for performing atmospheric correction on the remote sensing image to be corrected comprises the steps of:

[0030] A1, construct an atmospheric correction coefficient lookup table that characterizes the relationship between surface reflectance and apparent reflectance at the top of the atmosphere under different atmospheric conditions, observation conditions, and environmental conditions;

[0031] A2, look up the atmospheric correction coefficient lookup table to obtain the atmospheric correction coefficients Q corresponding to the atmospheric conditions, observation conditions, and environmental conditions of the study area. a , Q b , Q c ;

[0032] A3, atmospheric correction coefficient Q obtained by looking up the table a , Q b , Q c , correcting the surface reflectivity of the remote sensing image to be corrected taken for the study area.

[0033] Preferably, in step A3, the surface reflectivity of the remote sensing image to be corrected is corrected by the following formulas (4)-(5):

[0034]

[0035] in:

[0036] P=Q a LQ b (5)

[0037] In formulas (4)-(5), P is an intermediate variable;

[0038] L is the radiant brightness value after calibration.

[0039] Preferably, in step S1, the method for fusing the first remote sensing image and the second remote sensing image comprises the steps of:

[0040] B1, calculate the reflectivity of the central pixel in the moving window. The calculation method is expressed by the following formula (6):

[0041]

[0042] In formula (6), Represents the reflectivity of Landsat of the central pixel in the moving window;

[0043] C(x i ,y i ,t n ,B) represents the reflectivity of MODIS of the ith similar pixel;

[0044] w i represents the weight of the i-th similar pixel in the first remote sensing image and the second remote sensing image;

[0045] v i represents the conversion coefficient of the i-th similar pixel;

[0046] N represents the number of types of features;

[0047] Represents time as t m The reflectivity of Landsat at the center pixel in the moving window;

[0048] (x i ,y i ) represents the position of the i-th similar pixel;

[0049] t n ,t m They represent different times in the same area.

[0050] B represents the band of the image;

[0051] Indicates the center pixel position;

[0052] B2, searching similar pixels in the first remote sensing image and the second remote sensing image with the moving window, and then calculating the spatial weight and temporal weight between similar pixels, wherein the calculation method of the spatial weight between similar pixels is expressed by the following formula (8):

[0053] D i =(1-R i )×d i

[0054] in:

[0055] Where D i is the spatial weight, R i is the spectral similarity weight, that is, the correlation coefficient between the ith spectral similarity pixel and its corresponding low spatial resolution pixel, and its value range is (-1,1); d i is the geographical distance between similar pixels and the central pixel, Indicates the center pixel position, (x i ,yi ) represents the position of the i-th similar pixel.

[0056] The calculation method of the time weight between similar pixels is expressed by the following formula (7):

[0057]

[0058] In formula (7), T k Indicates time t k With the prediction time t p The time weight for calculating the degree of change in reflectivity of MODIS images between them;

[0059] C(x i ,y i ,t k , B) indicates the i-th similar pixel at t k MODIS reflectivity over time;

[0060] C(x i ,y i ,t p , B) indicates the i-th similar pixel at t p MODIS reflectivity over time;

[0061] t k represents the kth moment, t p represents the pth moment;

[0062] B3, t is obtained by the following formula (9): [ The fused image at time:

[0063]

[0064] In formula (9), Indicates t [ The central pixel value of the fused image at each moment;

[0065] Indicates t p The central pixel value of the fused image calculated by formula (6) at the moment;

[0066] Indicates t [ The central pixel value of the fused image is calculated by formula (6) at the moment.

[0067] Preferably, in step S2, the method for decomposing the mixed pixels in the fused image comprises the steps of:

[0068] C1, extracting pure land class end members in the fused image, the extraction method is: calculating the spectral index of each land class pixel in the fused image, and identifying the pixels with spectral index greater than or equal to or less than or equal to a preset spectral index threshold as pure land class end members of the corresponding landscape land class;

[0069] C2, identifying the pixels at the junction of pure end elements of different landscape types that are smaller than the spectral index threshold of the corresponding landscape type as the mixed pixels, and then further identifying sub-pixels from the mixed pixels, and calculating the area proportions and locations of the sub-pixels of different landscape types in the fused image.

[0070] Preferably, in step S3, the landscape pattern index includes any one or more of the patch area, landscape area ratio, and patch number on the fused image.

[0071] Preferably, in step S3, the spatial dynamics of the landscape pattern is calculated by the following formulas (12)-(13):

[0072]

[0073] In formula (12), K i represents the spatial dynamics of the landscape pattern of the i-th landscape type;

[0074] S it2 , S it1 Respectively represent the area of ​​the i-th landscape type in the fused image at time t2 and time t1; where t2 is later than t1;

[0075]

[0076] In formula (13), S represents the comprehensive spatial dynamics of the landscape pattern of the fused image;

[0077] △S i-j It represents the sum of the landscape conversion areas from the i-th landscape type to the j-th landscape type;

[0078] S i represents the total area of ​​the i-th landscape type on the fused image;

[0079] t is the research period;

[0080] n represents the number of landscape types in the fused image;

[0081] In step S3, the temporal variation characteristics of the landscape pattern are calculated by the following formula (14):

[0082]

[0083] In formula (14): X t ,Y t are the latitude and longitude coordinates of the distribution center of a certain landscape type in the fused image of year t on the X-axis and Y-axis respectively; k , Y k are the longitude and latitude coordinates of the geometric center of the kth evaluation unit on the X-axis and Y-axis respectively; S ki is the area of ​​land used for designated functions (such as wetland types, etc.) of the kth evaluation unit; n is the number of landscape types.

[0084] The present invention fuses a first remote sensing image with high spatial resolution and a second remote sensing image with high temporal resolution into a fused image, decomposes the mixed pixels in the fused image into sub-pixels, and based on the decomposition results of the mixed pixels, at the sub-pixel level, calculates the landscape pattern index, the spatial dynamics of the landscape pattern, and the temporal change characteristics of the landscape pattern for the fused image, thereby achieving fine-grained and accurate monitoring of the landscape pattern changes of the remote sensing images. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments of the present invention. Obviously, the drawings described below are only some embodiments of the present invention, and for ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0086] Figure 1 is a diagram of implementation steps of a wetland landscape pattern method considering sub-pixel information provided by an embodiment of the present invention;

[0087] Figure 2 This is an example of searching for similar pixels on MODIS and Landsat remote sensing images. DETAILED DESCRIPTION

[0088] The technical solution of the present invention is further described below with reference to the accompanying drawings and through specific implementation methods.

[0089] Among them, the drawings are only used for illustrative explanations, and they only represent schematic diagrams rather than actual pictures, and should not be understood as limitations on this patent; in order to better illustrate the embodiments of the present invention, some parts of the drawings may be omitted, enlarged or reduced, and do not represent the size of the actual product; for those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted.

[0090] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right", "inner", "outer", etc. appear, the orientation or position relationship indicated is based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, the terms describing the position relationship in the drawings are only used for illustrative purposes and cannot be understood as a limitation on this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0091] In the description of the present invention, unless otherwise clearly specified and limited, if the term "connection" or the like appears to indicate the connection relationship between components, the term should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, it can be the internal connection of two components or the interaction relationship between two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0092] The wetland landscape pattern monitoring method considering sub-pixel information provided by the embodiment of the present invention includes the following three technical steps:

[0093] Step 1: Multivariate data fusion

[0094] Step 2: Mixed pixel decomposition

[0095] Step 3: Wetland landscape pattern analysis and monitoring

[0096] Step one includes:

[0097] 1. Preprocess high spatial resolution remote sensing images to obtain ten-meter-level surface reflectivity; preprocess high temporal resolution remote sensing images to obtain daily surface reflectivity;

[0098] 2. Fuse data of different temporal and spatial resolutions to obtain high spatial resolution spatiotemporal seamless surface reflectance.

[0099] Step 2 includes:

[0100] 1. Classification of wetland cover types;

[0101] 2. Automatically extract pure end members based on the combination of spectral characteristic index and wetland surface cover type;

[0102] 3. Analyze and obtain the proportion and position of sub-pixels inside the mixed pixel.

[0103] Step three includes:

[0104] 1. Calculation of landscape pattern index at sub-pixel scale;

[0105] 2. Spatial structure analysis of wetland landscape types;

[0106] 3. Time series analysis of the overall wetland landscape level.

[0107] The following describes how to implement each technical point in each step.

[0108] The technical process of step one is:

[0109] Take Landsat (high spatial resolution) and MODIS (high temporal resolution) remote sensing images as examples. The Landsat series satellites of a certain country provide remote sensing image data with high spatial resolution, spectral resolution and radiation resolution for monitoring the earth's resources and environment. The remote sensing image data provided in this embodiment is preferably downloaded from the country's public remote sensing image data download website, with a spatial resolution of 30 meters and a temporal resolution of 16 days. MODIS is an important remote sensing instrument of a certain country's National Space Administration, installed on a satellite. MODIS has a high temporal resolution, and satellites equipped with MODIS can provide daily remote sensing image data covering the world. The remote sensing image provided in this embodiment preferably uses MOD09GQ data from the MODIS series data products, downloaded from the country's public website. The MOD09GQ remote sensing data set contains two bands, red light and near-red light, with a spatial resolution of 250m and a temporal resolution of 1 day.

[0110] 1. Data preprocessing for Landsat remote sensing images with high spatial resolution and MODIS remote sensing images with high temporal resolution

[0111] The data preprocessing method for Landsat remote sensing images and / or MODIS remote sensing images includes any one or more of radiation correction, geometric correction, atmospheric correction, reprojection, and resampling. In this embodiment, radiation correction, geometric correction, and atmospheric correction are performed on Landsat remote sensing images to ensure the quality of remote sensing data. MODIS remote sensing images are reprojected and resampled to match their spatial resolution with that of Landsat remote sensing images.

[0112] 1) Radiation correction

[0113] Affected by atmospheric absorption, scattering, sensor response characteristics and other random factors, the imaging process may cause blurring and distortion of remote sensing images, and reduction of image resolution and contrast. In order to eliminate the impact of these external or internal factors on the quality of remote sensing images, this embodiment restores remote sensing images through radiation correction (eliminating various distortion processes attached to the radiation brightness in remote sensing images). The principle is: convert the DN value of the remote sensing image to be corrected (the electromagnetic radiation reflected or emitted by the ground object is received by the satellite sensor, and is converted into a grayscale value, i.e., the DN value, after photoelectric conversion, analog-to-digital conversion and other internal coefficient processing) into the true reflectance value of the corresponding pixel. The conversion method is expressed by the following formula (1):

[0114]

[0115] In formula (1), L λ is the solar radiation intensity at the top of the atmosphere;

[0116] Q col is the DN value recorded for the pixel; Q col The minimum and maximum values ​​are 0 (the DN value of the corresponding pixel = 0) and 255 (the DN value of the corresponding pixel = 255).

[0117] Q colmax Q col The maximum value of

[0118] LMAX λ Indicates the upper limit of the remote sensor spectral radiation value;

[0119] LMIN λ Indicates the lower limit of the remote sensor spectral radiation value;

[0120] LMAX λ 、LMIN λ It is a constant value recorded in the remote sensing parameter file of the remote sensor.

[0121] 2) Geometric correction

[0122] In this embodiment, the remote sensing image to be corrected is corrected by the method expressed by the following formulas (2)-(3):

[0123]

[0124] In formulas (2)-(3), X and Y represent the coordinate values ​​of a pixel in the remote sensing image to be corrected on the X-axis and Y-axis respectively;

[0125] F x 、F y Represents the resampling correction distortion function associated with the X-axis and the Y-axis respectively;

[0126] (μ, ν) represents the actual geographic coordinates of the pixel to be corrected on the X-axis and Y-axis in the correction image space;

[0127] μ p v j The product of the horizontal and vertical coordinates of the actual geographic coordinates;

[0128] p represents the order of the horizontal axis;

[0129] j represents the order of the ordinate;

[0130] n represents the order of the polynomial;

[0131] a pj and b pj are coefficients, and are the coefficients of the polynomial solved by the least squares method; the following is a pj and b pj The process of solving the problem using the least squares method is described in detail:

[0132] The coefficients of the polynomial are solved by establishing a set of equations using ground control points. More than 10 control points are selected and the known coordinates of these ground control points are substituted into formulas (2) and (3) to solve the coefficients of the polynomial.

[0133] 3) Atmospheric correction

[0134] In the process of inverting the real reflectivity / emissivity of the surface from the top of the atmosphere radiation information obtained by the remote sensor, in order to quantitatively invert the surface parameters (including surface reflectivity, emissivity, etc.), the influence of the atmosphere on the image quality must be eliminated from the remote sensing data to obtain the real surface radiation information. In this embodiment, the atmospheric correction process used to eliminate the influence of the atmosphere on the quality of remote sensing images is provided as follows:

[0135] First, with the help of the atmospheric radiation transfer model, by presetting different atmospheric conditions (such as aerosol optical thickness, single scattering albedo and phase function, etc.), observation conditions (such as observation altitude angle, observation zenith angle, observation time, etc.) and environmental conditions (such as solar altitude angle, solar zenith angle and altitude, etc.), when the sensing parameters of the remote sensor are known, the relationship between the surface reflectivity and the apparent reflectivity of the top of the atmosphere is calculated, and a certain step size is established to obtain the atmospheric correction coefficient Q a , Q b , Q c The lookup table.

[0136] For example, with the help of the 6S model, the preset atmospheric conditions are: when the solar zenith angle is 0°~17°, the step length is set to 3°, when the solar zenith angle is 24°~60°, the step length is set to 12°, for a total of 11 values; the satellite zenith angle is 0°~75° with a step length of 15°, for a total of 6 values; the aerosol optical thickness is 0~1.95, with irregular step length, for a total of 11 values; the surface altitude is 0~1.9km, with irregular step length, for a total of 9 values; and the first 7 bands of MODIS. The above parameter combinations are brought into the 6S radiation model to obtain a lookup table of atmospheric correction coefficients under different parameters.

[0137] Then, using the obtained parameter information of atmospheric conditions, observation conditions and environmental conditions in the study area, the atmospheric correction coefficient Q corresponding to the atmospheric conditions, observation conditions and environmental conditions in the study area is obtained by looking up the table. a , Q b , Q c .

[0138] The following are examples of table lookup methods:

[0139] Find the atmospheric correction coefficient based on the values ​​of different parameters.

[0140] Finally, the atmospheric correction coefficient Q is obtained by looking up the table a , Q b , Q c Calculate the surface reflectivity of the study area. The calculation formula is expressed as follows:

[0141]

[0142] in:

[0143] P=Q a LQ b (5)

[0144] In formulas (4)-(5), P is an intermediate variable;

[0145] L is the radiant brightness value after calibration.

[0146] P is calculated as follows:

[0147] It is calculated by formula (5), where Q a , Q b By looking up the table, we can get that L can be calculated by formula (1) to obtain the method of reprojecting the remote sensing image:

[0148] In this embodiment, remote sensing image processing tools (such as ArcGIS, QGIS, ENVI, etc.) are used to reproject the MODIS remote sensing image to a coordinate system consistent with the Landsat remote sensing image.

[0149] The technical principle of reprojection is briefly described as follows:

[0150] Reprojection converts spatial data from one projection coordinate system to another. The reprojection process is to ensure that spatial data from different sources or different projections can be overlaid, analyzed, and displayed in the same geographic coordinate system.

[0151] In this embodiment, the MODIS remote sensing image also needs to be resampled from a spatial resolution of 500 m to a spatial resolution of 30 m, so that the spatial resolution of the MODIS remote sensing image is the same as that of the Landsat remote sensing image.

[0152] The method of resampling MODIS remote sensing images from 500m spatial resolution to 30m spatial resolution is briefly described as follows:

[0153] Generate higher resolution images through interpolation methods (bilinear interpolation or cubic convolution interpolation). Redistribute lower resolution pixel values ​​into more and smaller pixels, thereby improving spatial resolution.

[0154] 2. Fusion of pre-processed Landsat remote sensing images and MODIS remote sensing images

[0155] 1) Calculate the reflectivity of the central pixel in the moving window

[0156] In this embodiment, by determining the size of the moving window and the number of similar pixels, a linear conversion coefficient between the Landsat remote sensing image and the MODIS remote sensing image is calculated based on the reflectivity data of the two.

[0157] For example, after registering MODIS remote sensing images with Landsat remote sensing images, the reflectivity deviation between MODIS remote sensing images and Landsat remote sensing images is generally caused only by system deviation (such as deviation caused by different sensors). In order to eliminate this system deviation, the moving window method of this embodiment searches for similar pixels on the MODIS remote sensing image and the Landsat remote sensing image within the moving window, and then incorporates the similar pixel information into the calculation of the pixel reflectivity of the predicted fused image.

[0158] In this embodiment, similar pixels refer to neighboring pixels with the same ground object type as the central pixel. Figure 2 , a brief description of the method of searching for similar pixels in MODIS remote sensing images and Landsat remote sensing images using the moving window method is given below:

[0159] like Figure 2As shown in the figure, assuming that the moving window size is 3x3, the middle pixel is the central pixel. Finding similar pixels is to search all the pixels in the window to find the pixels with the same ground object type as the central pixel.

[0160] The reflectivity of the central pixel in the moving window is calculated by the following formula (6):

[0161]

[0162] In formula (6), Represents the reflectivity of Landsat of the central pixel in the moving window;

[0163] C(x i ,y i ,t n ,B) represents the reflectivity of MODIS of the ith similar pixel;

[0164] w i represents the weight of the i-th similar pixel in the first remote sensing image and the second remote sensing image, and represents the contribution of the i-th similar pixel to the reflectivity of the central pixel in the moving window, which is jointly determined by the spectral similarity of the Landsat remote sensing image and the MODIS remote sensing image at the position of the similar pixel and the distance between the similar pixel and the central pixel;

[0165] v i Represents the conversion coefficient of the ith similar pixel. The conversion coefficient is calculated based on image data at multiple moments. Usually, a linear regression model is used to establish the relationship between the spectral value and the time change. The formula is generally expressed as:

[0166] R(t)=αR HR (t1)+βR HR (t2)+γ

[0167] Among them, R(t) is the high-resolution reflectivity value at the prediction time, R HR (t1) and R HR (t2) are the reflectance values ​​of the high-resolution images at two known moments, α, β, γ are the conversion coefficients to be determined, and describe the changes in the spectrum at different moments.

[0168] N represents the number of types of features;

[0169] Represents time as t m The reflectivity of Landsat at the center pixel in the moving window;

[0170] (x i ,y i ) represents the position of the i-th similar pixel;

[0171] t n ,t m They represent different times in the same area.

[0172] B represents the band of the image;

[0173] Indicates the center pixel position;

[0174] 2) Calculation of spatiotemporal weight

[0175] Search similar pixels with the same feature type as the central pixel in the moving window and calculate their spatial weights. The calculation method of the spatial weights between similar pixels is briefly described as follows:

[0176] D i =(1-R i )×d i

[0177] in;

[0178] Where D i is the spatial weight, R i is the spectral similarity weight, that is, the correlation coefficient between the ith spectral similarity pixel and its corresponding low spatial resolution pixel, and its value range is (-1,1); R i The calculation method of is: through the spectral reflectance difference between the target pixel and the neighboring pixel, it is calculated, and the Euclidean distance is usually used to measure the similarity.

[0179] d i is the geographical distance between similar pixels and the central pixel, Indicates the center pixel position, (x i ,y i ) represents the position of the i-th similar pixel. In addition, the time weight between similar pixels is calculated based on the change in reflectivity of the MODIS remote sensing image before and after the prediction time. The calculation method is expressed by the following formula (7):

[0180]

[0181] In formula (7), T k Indicates time t k With the prediction time t p The temporal weight for calculating the degree of change in reflectivity of MODIS images between them;

[0182] C(x i ,y i ,t k , B) represents the i-th similar pixel t k MODIS reflectivity over time;

[0183] C(x i ,y i ,t p , B) represents the i-th similar pixel t p MODIS reflectivity over time;

[0184] t k represents time. Finally, t is obtained by the following formula (9): p The fused image is the result of fusing Landsat remote sensing images and MODIS remote sensing images at all times.

[0185]

[0186] In formula (9), Indicates t [ The central pixel value of the fused image at each moment;

[0187] represents the central pixel value of the fused image calculated by equation (6) at the time;

[0188] represents the central pixel value of the fused image calculated by formula (6) at the time. The technical process of step 2 is:

[0189] 1. Constructing a landscape classification system

[0190] According to the actual conditions of the study area, an appropriate landscape classification system was constructed, including different levels such as natural wetlands, artificial wetlands, and non-wetlands. Each level covers soil, water bodies, vegetation, and artificial surfaces.

[0191] 2. Land type end member extraction

[0192] Through landscape classification and spectral index calculation, pure land class end members are automatically selected according to the spectral index calculation results of different surface types. For example, dense vegetation pixels with NDVI (spectral index) > 0.7 are used as vegetation land class end members, and soil pixels with NDVI < 0.2 are used as soil land class end members.

[0193] 3. Based on the results of pure end-member extraction, the area proportions of sub-pixels of different surface types are extracted from the mixed pixels at the junction of different landscape types. At the same time, according to the spatial distribution, the positions of various sub-pixels are determined by the nearest neighbor principle. Sub-pixels refer to the more subtle ground feature information contained in a pixel in a remote sensing image.

[0194] Specifically, after the extraction of pure end members is completed, the pixels at the junction of pure end members of different land types are mixed pixels. For mixed pixels, the area proportions of sub-pixels of different surface types are calculated. The calculation method is: there are multiple types of land features in mixed pixels. For each type of land feature, their proportions in this pixel are calculated separately.

[0195] In addition, the position of the sub-pixel is determined based on the adjacent pure end members. The method of determining the position of the sub-pixel is as follows:

[0196] Generally, it is believed that mixed pixels are located at the junction of pure soil and pure vegetation, and that the position of soil in mixed pixels is close to pure soil pixels, and the position of vegetation is close to pure vegetation pixels.

[0197] The technical process of step three is:

[0198] 1. Calculation of landscape pattern index

[0199] Based on the mixed pixel decomposition results, the landscape pattern index is calculated, such as the patch area (CA), landscape area ratio (PLAND), number of patches (NP), etc. on the fused image.

[0200] Patch area (CA) refers to the area of ​​patches within a certain type of landscape. A patch refers to a geographical area or plot with certain relatively uniform characteristics in the landscape. A patch can be a forest, grassland, water body, etc. in a natural landscape, or it can be a land use type formed by human activities, such as farmland, urban buildings, etc.

[0201] The plaque area was calculated using the following formula (10):

[0202]

[0203] In formula (10), A i A is the area of ​​the i-th patch in the fused image, and n is the number of patches. i The calculation method is briefly described as follows:

[0204] Assuming that the area of ​​forest patches in the fused image is calculated, then A i It represents the area of ​​a forest patch in the fused image. Its calculation only requires the sum of the forest areas of all sub-pixels in the patch.

[0205] The landscape area ratio (PLAND) represents the ratio of the area of ​​the patch type to the total area of ​​the landscape, and is in the range of (0, 100). The closer the value is to 0, the sparser the patches in the landscape. When it is equal to 100, it means that the landscape is composed of one type of patches. The calculation of the landscape area ratio is expressed by the following formula (11):

[0206]

[0207] Among them, CA represents the area of ​​a patch of a certain land type, and A is the total area of ​​the entire landscape.

[0208] The patch number (NP) indicates the total number of patches in a certain patch or landscape area. The larger the NP, the higher the landscape fragmentation.

[0209] 2. Analysis of spatial dynamics of landscape pattern

[0210] In this embodiment, the landscape pattern dynamics of single landscape types and comprehensive landscape types at the sub-pixel scale are calculated to reveal the magnitude and speed of changes in landscape types during the study period.

[0211] In this embodiment, the spatial dynamics of a single type of landscape pattern is used to represent the quantitative changes of a certain type of landscape in a specific period of time, and it can quantitatively describe its change range and speed.

[0212] The dynamic degree of a single landscape type refers to the change in the quantity of a certain landscape type in the study area during the study period, and its calculation expression is:

[0213]

[0214] In formula (12), K i represents the spatial dynamics of the landscape pattern of the i-th landscape type;

[0215] S it2 , S it1 They represent the area of ​​the i-th landscape type in the fused image at time t2 and time t1 respectively; t2 is later than t1.

[0216] The following example shows K i The calculation method is explained in detail:

[0217] Assuming that we calculate the dynamics of a single landscape type in the two study periods t1 and t2 of the study area, we need to calculate the area of ​​a certain landscape type in the fused image of t1 and the area of ​​the fused image of t2. According to formula (12), K can be calculated: i

[0218] The spatial dynamics of the comprehensive landscape pattern of the fused image is calculated using the following formula (13):

[0219]

[0220] In formula (13), S represents the comprehensive spatial dynamics of the landscape pattern of the fused image;

[0221] △S i-j It represents the sum of the landscape conversion areas from the i-th landscape type to the j-th landscape type;

[0222] S i represents the total area of ​​the i-th landscape type on the fused image;

[0223] t is the research period;

[0224] n represents the number of landscape types in the fused image.

[0225] The following example illustrates the calculation method of the spatial dynamics S of the comprehensive landscape pattern of fused images:

[0226] Calculate the sum of the areas of all landscape types converted to other landscapes in the fused image. According to formula (13), the comprehensive landscape dynamics can be finally calculated.

[0227] △S i-j The calculation method is briefly described as follows:

[0228] The landscape type of the same area will change over time. For example, the original vegetation will turn into soil after a period of time. i-j It is to calculate how much vegetation in this area has been converted into soil.

[0229] 3. Analysis of temporal changes in landscape pattern

[0230] The moving window method and landscape pattern index were used to analyze the changing characteristics of landscape patterns at different time scales. The center of gravity of different landscape types was calculated, and the migration direction and speed of the center of gravity during the monitoring period were analyzed to explain the overall characteristics of different types of landscapes in terms of spatial changes.

[0231] In this embodiment, the temporal variation characteristics of the landscape pattern of the fused image are expressed by the following formula (14):

[0232]

[0233] In formula (14): X t ,Y t are the latitude and longitude coordinates of the distribution center of a certain landscape type in the fused image of year t on the X-axis and Y-axis respectively; k , Y k are the longitude and latitude coordinates of the geometric center on the X-axis and Y-axis of the k-th evaluation unit (for the study area, we usually divide it into different evaluation units according to the characteristics of the data. For example, the same area can be divided according to patches or grid units. Patches and grids are what we call evaluation units. The study area is usually divided into many blocks, the first block, the second block...the k-th block); S ki is the area of ​​land used for designated functions (such as wetland types, etc.) of the kth evaluation unit; n is the number of landscape types.

[0234] In summary, the present invention provides a wetland landscape pattern monitoring method considering sub-pixels, such as Figure 1 As shown, the steps include:

[0235] S1, after data preprocessing, a first remote sensing image with high spatial resolution and a second remote sensing image with high temporal resolution are fused into a fused image with high spatiotemporal resolution;

[0236] S2, decomposing the mixed pixels in the fused image to obtain the proportion and position of the sub-pixels inside the mixed pixels;

[0237] S3, based on the decomposition results of the mixed pixels in the fused image, calculate the landscape pattern index, and / or the spatial dynamics of the landscape pattern, and / or the temporal change characteristics of the landscape pattern of the fused image to complete the monitoring of the landscape pattern change.

[0238] In summary, the present invention fuses a first remote sensing image with high spatial resolution and a second remote sensing image with high temporal resolution into a fused image, decomposes the mixed pixels in the fused image into sub-pixels, and based on the decomposition results of the mixed pixels, at the sub-pixel level, calculates the landscape pattern index, the spatial dynamics of the landscape pattern, and the temporal change characteristics of the landscape pattern for the fused image, thereby achieving fine-grained and accurate monitoring of the landscape pattern changes of the remote sensing images.

[0239] It should be noted that the above specific implementations are only preferred embodiments of the present invention and the technical principles used. Those skilled in the art should understand that various modifications, equivalent substitutions, changes, etc. can be made to the present invention. However, as long as these changes do not deviate from the spirit of the present invention, they should be within the scope of protection of the present invention. In addition, some terms used in the specification and claims of this application are not restrictive, but are only for the convenience of description.

Claims

1. A wetland landscape pattern monitoring method considering sub-pixel information, characterized in that: Includes steps: S1, after data preprocessing, a first remote sensing image with high spatial resolution and a second remote sensing image with high temporal resolution are fused into a fused image with high spatiotemporal resolution; S2, decomposing the mixed pixel in the fused image to obtain the proportion and position of the sub-pixel inside the mixed pixel; S3, according to the decomposition result of the mixed pixels in the fused image, calculate the landscape pattern index, and / or the spatial dynamic degree of the landscape pattern, and / or the temporal change characteristics of the landscape pattern of the fused image, and complete the monitoring of the landscape pattern change.

2. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 1 is characterized in that: The first remote sensing image is a Landsat remote sensing image, and the second remote sensing image is a MODIS remote sensing image.

3. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 1 is characterized in that: In step S1, the data preprocessing method for the first remote sensing image and / or the second remote sensing image includes any one or more of radiation correction, geometric correction, atmospheric correction, reprojection, and resampling; The method for performing the radiation correction on the first remote sensing image and / or the second remote sensing image is expressed by the following formula (1): In formula (1), L λ is the solar radiation intensity at the top of the atmosphere; Q col is the DN value recorded for the pixel; Q col max is the Q recorded for all pixels in the remote sensing image to be corrected col The maximum value of LMAX λ Indicates the upper limit of the remote sensor spectral radiation value; LMIN λ Indicates the lower limit of the remote sensor spectral radiation value; LMAX λ 、LMIN λ It is a constant value recorded in the remote sensing parameter file of the remote sensor.

4. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 3 is characterized in that: The method for geometrically correcting the remote sensing image to be corrected is expressed by the following formulas (2)-(3): In formulas (2)-(3), X and Y represent the coordinate values ​​of a pixel in the remote sensing image to be corrected on the X-axis and Y-axis respectively; F x 、F y Represents the resampling correction distortion function associated with the X-axis and the Y-axis respectively; (μ,v) represents the actual geographic coordinates of the pixel to be corrected on the X-axis and Y-axis in the corrected image space; μ p v j Represents the product of the horizontal and vertical coordinate orders of the actual geographic coordinates; p represents the order of the horizontal axis; j represents the order of the ordinate; n represents the order of the polynomial; a pj and b pj is the coefficient.

5. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 3 is characterized in that: The method for performing atmospheric correction on the remote sensing image to be corrected comprises the steps of: A1, construct an atmospheric correction coefficient lookup table that characterizes the relationship between surface reflectance and apparent reflectance at the top of the atmosphere under different atmospheric conditions, observation conditions, and environmental conditions; A2, look up the atmospheric correction coefficient lookup table to obtain the atmospheric correction coefficients Q corresponding to the atmospheric conditions, observation conditions, and environmental conditions of the study area. a , Q b , Q c ; A3, atmospheric correction coefficient Q obtained by looking up the table a , Q b , Q c , correcting the surface reflectivity of the remote sensing image to be corrected taken for the study area.

6. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 5 is characterized in that: In step A3, the surface reflectivity of the remote sensing image to be corrected is corrected by the following formulas (4)-(5): in: P=Q a LQ b (5) In formulas (4)-(5), P is an intermediate variable; L is the radiant brightness value after calibration.

7. The wetland landscape pattern monitoring method considering sub-pixel information according to any one of claims 1 to 6, characterized in that: In step S1, the method for fusing the first remote sensing image and the second remote sensing image comprises the steps of: B1, calculate the reflectivity of the central pixel in the moving window. The calculation method is expressed by the following formula (6): In formula (6), Represents the reflectivity of Landsat of the central pixel in the moving window; C(x i ,y i , t n , B) represents the reflectivity of MODIS of the i-th similar pixel; w i represents the weight of the i-th similar pixel in the first remote sensing image and the second remote sensing image; v i represents the conversion coefficient of the i-th similar pixel; N represents the number of types of features; Represents time as t m The reflectivity of Landsat at the center pixel in the moving window; (x i ,y i ) represents the position of the i-th similar pixel; t n ,t m They represent different times in the same area. B represents the band of the image; Indicates the center pixel position; B2, searching similar pixels in the first remote sensing image and the second remote sensing image with the moving window, and then calculating the spatial weight and temporal weight between similar pixels, wherein the calculation method of the spatial weight between similar pixels is expressed by the following formula (8): D i =(1-R i )×d i in: Where D i is the spatial weight, R i is the spectral similarity weight, that is, the correlation coefficient between the ith spectrally similar pixel and its corresponding low spatial resolution pixel, and its value range is (-1, 1); d i is the geographical distance between similar pixels and the central pixel, represents the center pixel position, (x i ,y i ) represents the position of the i-th similar pixel; The calculation method of the time weight between similar pixels is expressed by the following formula (7): In formula (7), T k Indicates time t k With the prediction time t p The time weight for calculating the degree of change in reflectivity of MODIS images between them; C(x i ,y i , t k , B) indicates the i-th similar pixel at t k MODIS reflectivity over time; C(x i ,y i , t p , B) indicates the i-th similar pixel at t p MODIS reflectivity over time; t k represents the kth moment, t p represents the pth moment; B3, t is obtained by the following formula (9): p The fused image at time: In formula (9), Indicates t p The central pixel value of the fused image at each moment; Indicates t p The central pixel value of the fused image calculated by formula (6) at the moment; Indicates t p The central pixel value of the fused image is calculated by formula (6) at the moment.

8. The wetland landscape pattern monitoring method considering sub-pixel information according to any one of claims 1 to 6, characterized in that: In step S2, the method for decomposing the mixed pixels in the fused image comprises the steps of: C1, extracting pure land class end members in the fused image, the extraction method is: calculating the spectral index of each land class pixel in the fused image, and identifying the pixels whose spectral index is greater than or equal to or less than or equal to a preset spectral index threshold as pure land class end members of the corresponding landscape land class; C2, identifying the pixels at the junction of pure end elements of different landscape types that are smaller than the spectral index threshold of the corresponding landscape type as the mixed pixels, and then further identifying sub-pixels from the mixed pixels, and calculating the area proportions and locations of the sub-pixels of different landscape types in the fused image.

9. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 1, characterized in that: In step S3, the landscape pattern index includes any one or more of the patch area, landscape area ratio, and patch number on the fused image.

10. The wetland landscape pattern monitoring method considering sub-pixel information according to claim 1, characterized in that: In step S3, the spatial dynamics of the landscape pattern is calculated using the following formulas (12)-(13): In formula (12), K i represents the spatial dynamics of the landscape pattern of the i-th landscape type; S it2 , S it1 Respectively represent the area of ​​the i-th landscape type in the fused image at time t2 and time t1; where t2 is later than t1; In formula (13), S represents the comprehensive spatial dynamics of the landscape pattern of the fused image; △S i-j It represents the sum of the landscape conversion areas from the i-th landscape type to the j-th landscape type; S i represents the total area of ​​the i-th landscape type on the fused image; t is the research period; n represents the number of landscape types in the fused image; In step S3, the temporal variation characteristics of the landscape pattern are calculated by the following formula (14): In formula (14): X t ,Y t are the latitude and longitude coordinates of the distribution center of a certain landscape type in the fused image of year t on the X-axis and Y-axis respectively; k , Y k are the longitude and latitude coordinates of the geometric center of the kth evaluation unit on the X-axis and Y-axis respectively; S ki is the area of ​​land used for designated functions (such as wetland types, etc.) of the kth evaluation unit; n is the number of landscape types.

Images