Method for detecting greenhouse gas emissions from grazing natural grasslands

Abstract

The present application relates to the technical field of grassland ecological monitoring, and specifically discloses a detection method for greenhouse gas emission of grazing natural grassland, which comprises the following steps: arranging a ground monitoring array, continuously collecting greenhouse gas flux, precipitation intensity and soil temperature and humidity data, and obtaining multi-temporal, multi-spectral and thermal infrared images; extracting flux peak value and time phase images before and after precipitation according to precipitation events, synthesizing remote sensing feature maps with the same phase as the flux peak value through time sequence synthesis; correlatively storing the synthesized remote sensing feature maps and the flux peak value to form a pulse synchronous paired data set, and then correcting the benchmark flux estimation rule; applying the updated rule to the whole grassland remote sensing image, calculating the flux pixel by pixel and generating an emission spatial distribution map; the present application solves the problem that unmanned aerial vehicles cannot capture the short-time emission peak value after precipitation due to meteorological restrictions, and realizes dynamic detection of greenhouse gas emission of grazing natural grassland.

Description

Technical Field

[0001] This invention relates to the field of grassland ecological monitoring technology, specifically to a method for detecting greenhouse gas emissions from grazing natural grasslands. Background Technology

[0002] Grazing grasslands are important sources and sinks of greenhouse gases in terrestrial ecosystems, and their emissions of carbon dioxide, methane, and nitrous oxide are significantly influenced by grazing intensity, soil moisture changes, and precipitation events. In recent years, ground-based sensor networks and UAV remote sensing technologies have been increasingly applied to the detection of greenhouse gas emissions from grasslands. Ground-based sensor networks can continuously and frequently collect flux data and environmental parameters at specific locations, while UAVs equipped with multispectral and thermal infrared cameras can acquire grassland ecological parameters (such as vegetation cover, biomass, and soil moisture) over a wide area with high spatial resolution. Combining these two technologies and using regression or machine learning models, flux estimation can be achieved from point scale to area scale.

[0003] In monitoring grazing grasslands, limitations imposed by meteorological and airspace conditions on UAV remote sensing prevent the acquisition of synchronous remote sensing images within the peak window of greenhouse gas pulse emissions (6 to 12 hours after precipitation). This results in a long-term temporal misalignment between the short-term high emission fluxes captured by ground sensors and the corresponding remote sensing features. Existing technologies cannot construct paired remote sensing-flux samples at the pulse peak state, thus hindering the correction of the nonlinear mapping relationship between transient emissions triggered by precipitation events and grassland spectral characteristics using conventional regression or machine learning methods. This invention effectively solves the problem of pulse emission detection failure caused by missing temporal sampling by synthesizing remote sensing feature maps contemporaneous with flux peaks and using these maps to form a pulse synchronous paired dataset to update flux estimation rules. Summary of the Invention

[0004] The purpose of this invention is to provide a method for detecting greenhouse gas emissions from grazing natural grasslands, in order to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions:

[0006] A method for detecting greenhouse gas emissions from grazing natural grasslands includes the following steps:

[0007] S1: Deploy ground monitoring arrays on grazing natural grasslands to continuously collect data on greenhouse gas fluxes, precipitation intensity, and soil temperature and humidity. At the same time, acquire multi-temporal multispectral and thermal infrared images of the grasslands to form ground flux sequences and remote sensing image sequences.

[0008] S2: Identify precipitation events based on precipitation intensity, extract the peak flux in the first time period after precipitation from the ground flux sequence, and extract the first time phase image before precipitation and the second time phase image after precipitation from the remote sensing image sequence; through time-series synthesis, generate a synthetic remote sensing feature map that is in the same phase as the peak flux using the first time phase image, the second time phase image, and the high-resolution flux sequence before and after the peak flux.

[0009] S3: Associate and store the synthetic remote sensing feature map with the flux peak to form a pulse synchronization paired dataset;

[0010] S4: Based on the pulse synchronization paired dataset, the mapping relationship between remote sensing features and flux in the baseline flux estimation rule is corrected to obtain the updated flux estimation rule;

[0011] S5: Apply the updated flux estimation rules to the whole grassland remote sensing image, calculate greenhouse gas flux pixel by pixel, and generate a spatial distribution map of greenhouse gas emissions.

[0012] As a further aspect of the present invention: the extraction of flux peak values ​​from the ground flux sequence within the first time period after precipitation, and the extraction of the first temporal image before precipitation and the second temporal image after precipitation from the remote sensing image sequence, specifically includes:

[0013] The soil moisture change rate after the falling edge of precipitation intensity is extracted from the surface flux sequence. When the change rate turns from positive to negative and the absolute value exceeds the standard deviation of historical background fluctuations, the flux value corresponding to the turning point is taken as the flux peak.

[0014] The first phase image is the image frame with the smallest spatial variance of surface temperature in the most recent thermal infrared image before the precipitation begins.

[0015] The second temporal image is the image frame corresponding to the moment when the slope of the vegetation red edge returns to the pre-precipitation level in the first multispectral image after the precipitation ends.

[0016] As a further aspect of the present invention: the process of generating the synthetic remote sensing feature map is as follows:

[0017] The first and second temporal images were subjected to non-downsampled contour wave transforms to decompose the low-frequency subband coefficients and high-frequency directional subband coefficients at multiple scales.

[0018] The flux ratio between the peak time and the previous time in the high-resolution flux sequence before and after the flux peak is used as the low-frequency weighting factor, and the flux ratio between the peak time and the next time is used as the directional weighting factor for the high-frequency directional subband coefficient.

[0019] The low-frequency weighting factor is multiplied by the sum of the low-frequency sub-band coefficients of the first and second temporal images, and the directional weighting factor is multiplied by the difference between the high-frequency directional sub-band coefficients of the corresponding directions of the first and second temporal images. The product result is then subjected to non-downsampled contour wave inverse transform to reconstruct the synthetic remote sensing feature map.

[0020] As a further aspect of the present invention: the non-subsampled contour wave transformation specifically includes:

[0021] Calculate the local grayscale variance maps of the first and second temporal images, use the absolute value of the difference between the two local grayscale variance maps as the weight map, and then determine the number of decomposition scale levels of the non-subsampled contour wave transform based on the energy distribution of the weight map. The number of directional sub-bands at each scale level is determined by the number of peaks in the polar coordinate sector after the Fourier transform of the weight map at the corresponding level.

[0022] Based on the number of decomposition scale levels and the number of directional subbands at each level, non-subsampled Laplacian pyramid decomposition and non-subsampled directional filter bank decomposition are performed on each image to output low-frequency subband coefficients and high-frequency directional subband coefficients at multiple scales.

[0023] As a further aspect of the present invention: S3 specifically includes:

[0024] A two-dimensional hash map is constructed using the spatial pixel locations of the synthesized remote sensing feature maps as row indices and the occurrence time of flux peaks as column indices.

[0025] The reflectance value vector of each pixel in the synthetic remote sensing feature map and the flux peak value at the corresponding pixel location are stored as key-value pairs in a two-dimensional hash mapping table.

[0026] Using the geographical coordinates and time window of the precipitation event as the query key, all associated key-value pairs are extracted and concatenated in chronological order to form a pulse synchronization paired dataset.

[0027] As a further aspect of the present invention: the construction of the two-dimensional hash mapping table specifically includes:

[0028] For each spatial pixel location, calculate its geospatial grid code of latitude and longitude coordinates, and use the geospatial grid code as the row key value;

[0029] For each flux peak occurrence time, it is converted into a continuous time sequence number from the reference time and mapped to a column key value using a prime number hash function;

[0030] The row key and column key are concatenated to form a composite key. The reflectance vector of the synthetic remote sensing feature map at the corresponding pixel position and the corresponding flux peak are used as the composite value and stored in the hash bucket to complete the construction of the two-dimensional hash mapping table.

[0031] As a further aspect of the present invention: S4 specifically includes:

[0032] Extract the pixel reflectance vector and its corresponding flux peak of each synthetic remote sensing feature map from the pulse synchronization pairing dataset, and calculate the residual between the flux peak and the predicted flux output by the baseline flux estimation rule.

[0033] Using the Mahalanobis distance between pixel reflectance vectors in feature space as weights, the residuals are inversely distributed to the mapping coefficients of each feature component in the baseline rule according to the distance. This process is iterated until the sum of squared residuals of all paired datasets converges, thus obtaining the updated mapping relationship.

[0034] As a further aspect of the present invention: the generation of the spatial distribution map of greenhouse gas emissions specifically includes:

[0035] The whole grassland remote sensing image is divided into multiple topographically homogeneous patches according to the topographic slope and aspect. The mode vector of the pixel reflectance vector in each patch is used as the representative feature of the patch.

[0036] The representative features are substituted into the updated flux estimation rules to calculate the representative flux of the patch. Then, the representative flux is linearly inversely distributed according to the spectral angular distance between the reflectance vector of each pixel and the representative features to obtain the greenhouse gas flux of each pixel.

[0037] The flux values ​​of all pixels are filled into the blank raster according to the geographic coordinates of the original image, and after median filtering and smoothing, a spatial distribution map of greenhouse gas emissions is output.

[0038] As a further aspect of the present invention: the calculation of the representative flux of the patch specifically includes:

[0039] The reflectance values ​​of each band in the representative features are multiplied one by one with the mapping coefficients of the corresponding bands in the updated flux estimation rules, and then summed to obtain the preliminary flux.

[0040] Then, retrieve the three paired samples with the smallest spectral angular distance to the representative feature from the pulse synchronization paired dataset, and calculate the average residual between the peak flux of the corresponding pixel of the synthesized remote sensing feature map in these three paired samples and the initial flux.

[0041] The residual average is multiplied by a weighting coefficient that decreases as the spectral angular distance increases, and then superimposed onto the initial flux to obtain the final representative flux.

[0042] The beneficial effects of this invention are:

[0043] (1) This invention compensates for the insufficiency of UAVs being unable to acquire actual remote sensing images within a short window after precipitation due to meteorological conditions such as cloud cover and airspace restrictions by synthesizing remote sensing feature maps that are in the same phase as the flux peak. This allows the original missing remote sensing features within 6 to 12 hours after precipitation to be reconstructed, thereby achieving effective detection of greenhouse gas pulse emissions caused by extreme precipitation events in grazing natural grasslands. This avoids the problem of key emission processes being ignored due to time sampling misalignment in traditional methods.

[0044] (2) This invention uses pulse synchronization paired dataset to correct the remote sensing features and flux mapping relationship in the baseline flux estimation rule, so that the updated flux estimation rule can better reflect the actual correlation between grassland spectral features and greenhouse gas flux under short-term high emission conditions after precipitation events. The resulting spatial distribution map can distinguish emission differences under different terrain and vegetation conditions, providing a more accurate monitoring basis for grazing management based on actual emission dynamics. Attached Figure Description

[0045] The invention will now be further described with reference to the accompanying drawings.

[0046] Figure 1 This is a flowchart of the method of the present invention;

[0047] Figure 2 This is a flowchart of the process of generating the synthetic remote sensing feature map in this invention. Detailed Implementation

[0048] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0049] Please see Figure 1 As shown, this invention provides a method for detecting greenhouse gas emissions from grazing natural grasslands, comprising the following steps:

[0050] S1: Deploy ground monitoring arrays on grazing natural grasslands to continuously collect data on greenhouse gas fluxes, precipitation intensity, and soil temperature and humidity. At the same time, acquire multi-temporal multispectral and thermal infrared images of the grasslands to form ground flux sequences and remote sensing image sequences.

[0051] S2: Identify precipitation events based on precipitation intensity, extract the peak flux in the first time period after precipitation from the ground flux sequence, and extract the first time phase image before precipitation and the second time phase image after precipitation from the remote sensing image sequence; through time-series synthesis, generate a synthetic remote sensing feature map that is in the same phase as the peak flux using the first time phase image, the second time phase image, and the high-resolution flux sequence before and after the peak flux.

[0052] S3: Associate and store the synthetic remote sensing feature map with the flux peak to form a pulse synchronization paired dataset;

[0053] S4: Based on the pulse synchronization paired dataset, the mapping relationship between remote sensing features and flux in the baseline flux estimation rule is corrected to obtain the updated flux estimation rule;

[0054] S5: Apply the updated flux estimation rules to the whole grassland remote sensing image, calculate greenhouse gas flux pixel by pixel, and generate a spatial distribution map of greenhouse gas emissions.

[0055] In S1, a ground-based monitoring array is deployed on grazing natural grasslands to continuously collect data on greenhouse gas fluxes, precipitation intensity, and soil temperature and humidity. Simultaneously, multi-temporal multispectral and thermal infrared images of the grasslands are acquired, forming ground flux sequences and remote sensing image sequences, specifically including:

[0056] A typical monitoring area of ​​1 square kilometer was selected within grazing natural grassland, and a ground monitoring array was deployed with a grid spacing of 200 meters × 200 meters. The array consists of 9 monitoring points, each with a multi-parameter sensor integrated column. Within each integrated column, a non-dispersive infrared absorption gas sensor is used to detect the concentrations of carbon dioxide, methane, and nitrous oxide in real time. Combined with a built-in micro-pump and a sealed gas chamber, near-surface gas samples are collected at a frequency of once per minute, and greenhouse gas fluxes are calculated based on the rate of change in gas concentration. A tipping bucket rain gauge is used to record precipitation intensity with a resolution of 0.2 millimeters. A thermistor and a frequency domain reflectometer are used to measure soil temperature (at depths of 5 cm and 15 cm) and soil volumetric water content, respectively. All sensor data are uploaded to a data aggregation terminal via a low-power wide-area network at a frequency of once per minute to form a ground flux sequence.

[0057] Simultaneously, a quadcopter drone equipped with a multispectral camera and a thermal infrared camera was used. The multispectral camera includes five bands: blue, green, red, red-edge, and near-infrared, while the thermal infrared camera operates in the 8–14 micrometer band. The drone flew every seven days according to a pre-planned route, under rain-free and cloud-free weather conditions, at an altitude of 80 meters, achieving a ground resolution of 10 centimeters. Each flight acquired multispectral and thermal infrared images covering the entire monitoring area, recording vegetation reflectance spectral information and surface radiation temperature, respectively. Images acquired from different flight dates were arranged chronologically to form a remote sensing image sequence. The ground flux sequence and the remote sensing image sequence were mapped one-to-one using GPS coordinates and acquisition timestamps, forming the foundation dataset for subsequent steps.

[0058] Please see Figure 2 As shown, in S2, precipitation events are identified based on precipitation intensity. The peak flux value in the first time period after precipitation is extracted from the surface flux sequence, and the first-phase image before precipitation and the second-phase image after precipitation are extracted from the remote sensing image sequence. Through temporal synthesis, a synthetic remote sensing feature map synchronous with the flux peak is generated using the first-phase image, the second-phase image, and the high-resolution flux sequence before and after the flux peak. Specifically, this includes:

[0059] Precipitation intensity data is continuously monitored from the surface flux sequence. The moment when the precipitation intensity changes from greater than 0 mm / h to 0 mm / h is defined as the end of precipitation. Starting from the end of precipitation, the time period from the moment when the precipitation intensity began to be greater than 0 mm / h is defined as the precipitation duration. After the precipitation duration ends, the soil moisture change rate is extracted from the surface flux sequence. The change rate is calculated as: the current soil volumetric water content minus the soil volumetric water content one minute prior, divided by the sampling interval (1 minute). The soil moisture change rate is evaluated point by point within the first 6 hours after the precipitation ends. When the change rate changes from a positive value (indicating continuous soil water absorption) to a negative value (indicating the beginning of soil water loss), and the absolute value of this negative change rate exceeds the historical background fluctuation standard deviation, the moment corresponding to this inflection point is recorded as the peak moment, and the greenhouse gas flux value measured by the surface sensor at this moment is taken as the flux peak value. The historical background fluctuation standard deviation is obtained by selecting three consecutive days of data with no precipitation and stable grazing intensity in the past 30 days, and calculating the standard deviation of the soil moisture change rate per minute during this period as the historical background fluctuation standard deviation.

[0060] From the remote sensing image sequence, locate the most recent thermal infrared image before the start of precipitation (i.e., the first minute before precipitation intensity changes from 0 to greater than 0). Calculate the spatial variance of surface temperature for this thermal infrared image: divide the entire image into a 10m × 10m grid, calculate the average pixel temperature within each grid, and then use the average values ​​of all grids as a sample to calculate the overall variance. If multiple thermal infrared images exist before the start of precipitation (e.g., every 7 days), select the frame with the smallest spatial variance of surface temperature as the first temporal image. This image represents the most uniform spatial distribution of surface temperature before precipitation, which helps reduce spatial heterogeneity interference during subsequent image synthesis.

[0061] From the remote sensing image sequence, the first multispectral image appearing after the end of precipitation (i.e., the first minute after the precipitation intensity changes from greater than 0 to 0) is located. The vegetation red edge slope is calculated for this multispectral image, defined as the difference between the near-infrared and red band reflectance divided by the difference in the center wavelengths of the two bands (near-infrared center wavelength 800 nm, red band center wavelength 660 nm). Simultaneously, the red edge slope is calculated from the multispectral image of the most recent rainless day before the precipitation began, serving as the pre-precipitation level. The red edge slopes of each multispectral image after the precipitation ended are compared frame by frame. When the red edge slope value of a certain frame recovers to 95% to 105% of the pre-precipitation level, that frame is used as the second temporal image.

[0062] The first and second temporal images were subjected to non-downsampled contourlet transforms. Before the transform, the local gray-level variance map of each image was calculated: a 9-pixel × 9-pixel window was taken with each pixel as the center, the variance of the gray-level values ​​of all pixels within the window was calculated, and this variance was assigned to the center pixel. This process was repeated across the entire image to obtain the local gray-level variance map. The local gray-level variance map of the first temporal image was subtracted pixel by pixel from the local gray-level variance map of the second temporal image, and the absolute value was taken to obtain the weight map. A two-dimensional Fourier transform was performed on the weight map, and the transformed spectrum map was converted to polar coordinate sectors (15-degree angular intervals per sector). The cumulative sum of the spectral amplitudes within each sector was counted, and the number of sectors containing the peak of the cumulative sum is the number of directional sub-bands at that scale level. The method for determining the number of decomposition scale levels is as follows: Starting from the finest scale, calculate the energy of the weight map at the corresponding scale (i.e., the sum of squares of all pixel gray levels at that scale). Decomposition stops when the energy at a certain scale is less than 5% of the energy at the finest scale. The number of levels before that scale is the number of decomposition scale levels. Based on the determined number of decomposition scale levels and the number of directional subbands corresponding to each level, perform non-downsampled Laplacian pyramid decomposition (to obtain low-frequency subband coefficients at different scales) and non-downsampled directional filter bank decomposition (to obtain high-frequency directional subband coefficients in each direction at each scale) on the first and second temporal images, respectively. The output results include: low-frequency subband coefficients at multiple scales and high-frequency directional subband coefficients in multiple directions at each scale for the first temporal image, as well as two sets of coefficients for the second temporal image.

[0063] High-resolution flux sequences before and after the flux peak are extracted. These sequences include flux values ​​one minute before the peak, flux values ​​at the peak, and flux values ​​one minute after the peak. The ratio of the flux value at the peak to the flux value at the previous moment is used as a low-frequency weighting factor. The ratio of the flux value at the peak to the flux value at the next moment is used as a direction weighting factor. For each scale, the corresponding low-frequency subband coefficients of the first and second temporal images are added pixel-by-pixel and multiplied by the low-frequency weighting factor to obtain the synthesized low-frequency subband coefficients. For each direction at each scale, the high-frequency directional subband coefficients of the second temporal image are subtracted from the corresponding high-frequency directional subband coefficients of the first temporal image and multiplied by the direction weighting factor to obtain the synthesized high-frequency directional subband coefficients. The synthesized low-frequency subband coefficients and synthesized high-frequency directional subband coefficients at various scales and directions are input into the inverse non-downsampled contour wave transform (the inverse transform process involves first combining the directional subband coefficients using a non-downsampled directional filter, and then inversely reconstructing the Laplace pyramid). This reconstructs an image of the same size as the first and second temporal images, which is the synthesized remote sensing feature map. The spatial resolution of this feature map is consistent with the original remote sensing image, and its pixel values ​​reflect the simulated remote sensing characteristics within 6 to 12 hours after precipitation (i.e., the time when the flux peak occurs), including equivalent multispectral reflectance and surface radiant temperature.

[0064] In S3, the synthetic remote sensing feature maps are associated and stored with flux peaks to form a pulse synchronization paired dataset, specifically including:

[0065] A two-dimensional hash map is constructed using the spatial pixel locations of the synthesized remote sensing feature maps as row indices and the occurrence time of flux peaks as column indices. The specific construction method is as follows:

[0066] First, for each pixel location in the synthetic remote sensing feature map, obtain the corresponding longitude and latitude coordinates (provided by the georegistration file of the UAV imagery). Convert these coordinates to geospatial grid codes using the following method: recursively divide the global map into 32 layers using a quartering method, with each layer's grid size being half that of the previous layer. Take the 20th grid layer as the coding baseline layer; its side length is approximately 0.3 meters. Calculate the row and column numbers of the pixel's latitude and longitude coordinates within the 20th grid layer (row numbers increase northward from the equator, and column numbers increase eastward from the Prime Meridian). Then, represent the row and column numbers in binary, and merge them bit-wise (i.e., first take the first binary digit of the row number, then the first binary digit of the column number, and so on until all bits are taken). Convert the merged binary number to a decimal integer; this integer is the geospatial grid code for that pixel. Use this geospatial grid code as the row key of a hash table.

[0067] Calculate the total number of minutes elapsed between the peak flux occurrence and the reference time to obtain a continuous time sequence number (an integer). Select a prime number hash function, specifically the prime number 1000003, as the hash base, and use the remainder when the continuous time sequence number is divided by 1000003 as the column key value.

[0068] The row key and column key are concatenated as a string in the form of "row key - column key" to form a composite key. For each pixel location in the synthetic remote sensing feature map, the reflectance values ​​of that pixel location in the five multispectral bands (blue, green, red, red edge, and near-infrared) are extracted to form a reflectance value vector containing five elements; at the same time, the ground flux peak value corresponding to that pixel location is extracted (this value has been obtained in step S2, and the same synthetic remote sensing feature map corresponds to a unique flux peak time, so the flux peak value of each pixel is the same). The reflectance value vector and the flux peak value are used together as a composite value. The composite key and composite value are stored as a key-value pair in a hash bucket. The hash bucket uses open addressing to resolve collisions, and the size of the bucket is set to 1.2 times the total number of pixels in the synthetic remote sensing feature map. The above operation is repeated for all pixel locations until all pixels have been processed, completing the construction of the two-dimensional hash mapping table.

[0069] For each precipitation event, the geographic coordinate range (i.e., the smallest bounding rectangle of the 1 square kilometer area covered by the ground monitoring array) and time window (i.e., the time interval from the start of precipitation to one hour after the flux peak) are obtained. Using this geographic coordinate range and time window as the query key, a search is performed in a two-dimensional hash map table: first, all possible geospatial grid-coded intervals within the geographic coordinate range are calculated; then, the corresponding column key value intervals are calculated based on the flux peak occurrence time within the time window. All key-value pairs that simultaneously satisfy both row and column key value interval conditions are selected. The retrieved key-value pairs are concatenated and sorted according to the flux peak occurrence time (i.e., the original time corresponding to the column key value) to form a data record arranged in chronological order. This data record constitutes the pulse synchronization paired dataset corresponding to a single precipitation event. Data sets from multiple precipitation events together constitute the overall pulse synchronization paired dataset, used for subsequent flux estimation rule correction.

[0070] In S4, based on the pulse synchronization paired dataset, the mapping relationship between remote sensing features and flux in the baseline flux estimation rule is corrected to obtain the updated flux estimation rule, which specifically includes:

[0071] A baseline flux estimation rule is predefined, which states that the greenhouse gas flux at any pixel location is equal to the sum of the reflectance values ​​of that pixel in the five multispectral bands (bands 1 to 5, corresponding to blue, green, red, red-edge, and near-infrared, respectively) multiplied by their respective mapping coefficients. The initial values ​​of the five mapping coefficients are obtained through linear regression of flux and remote sensing data during historical periods of no precipitation. Specifically, all remote sensing images with stable grazing intensity and no precipitation over the past year, along with their synchronized ground flux data, are selected. The least squares method is used to calculate the five coefficient values ​​that minimize the sum of squared residuals between the predicted and measured fluxes, which serve as the baseline mapping coefficients. The pixel reflectance vector (containing the reflectance values ​​of the five bands) and its corresponding flux peak value of each synthetic remote sensing feature map are read pixel by pixel from the pulse-synchronized paired dataset. The five-band reflectance values ​​of the pixel are substituted into the baseline flux estimation rule, and a predicted flux value is calculated by multiplying the five band reflectance values ​​by their respective mapping coefficients and then summing the results. The residual is obtained by subtracting the predicted flux value from the peak flux value corresponding to the pixel.

[0072] Calculate the Mahalanobis distance of the pixel's reflectance vector in the feature space. The feature space is a five-dimensional space composed of reflectance values ​​from five bands. First, extract the five-band reflectance vectors of all pixels from the pulse synchronization pairing dataset, and calculate the covariance matrix of this set of vectors: for each pair of bands, calculate the covariance of the reflectance values ​​of all pixels in these two bands, forming a 5x5 covariance matrix. For the current pixel's five-band reflectance vector, subtract it from the average vector of all pixels (the mean of the five bands) to obtain a difference vector; perform two inner product operations (i.e., first multiply by the inverse matrix on the left and then multiply by the transpose of the difference vector on the right) on this difference vector, and the square root of the result is the Mahalanobis distance. The larger the Mahalanobis distance, the further the pixel's reflectance vector deviates from the overall distribution, and the lower its reliability.

[0073] The reciprocal of the Mahalanobis distance is used as the basis for weight allocation. The residual calculated above is multiplied by (the reciprocal of the Mahalanobis distance of the pixel divided by the sum of the reciprocals of the Mahalanobis distances of all pixels) to obtain a weighted residual. This weighted residual is then allocated to the five mapping coefficients according to the proportion of the reflectance values ​​of the five bands to the total reflectance of all bands: for each band, the reflectance value of that band is divided by the sum of the reflectance values ​​of the five bands, and then multiplied by the weighted residual to obtain the adjustment amount of the mapping coefficient for that band. This adjustment amount is added to the baseline mapping coefficient to obtain the new mapping coefficient after one iteration. The above residual calculation, Mahalanobis distance weight allocation, and mapping coefficient adjustment are performed sequentially on all pixels in the pulse synchronization pairing dataset to complete one iteration. After one iteration, the predicted flux of all pixels is recalculated using the updated mapping coefficients, and the sum of squares of the residuals of all pixels is calculated. This sum of squares is compared with the sum of squares of the residuals from the previous iteration. Convergence is determined when the difference between the two sums of squares is less than one-thousandth of the sum of squares from the previous iteration. The five mapping coefficients obtained upon convergence form the mapping relationships in the updated flux estimation rule. If convergence fails, the updated mapping coefficients are used as the new baseline mapping coefficients, and the above iterative process is repeated until the convergence condition is met.

[0074] In S5, the updated flux estimation rules are applied to whole-grassland remote sensing imagery to calculate greenhouse gas fluxes pixel by pixel, generating a spatial distribution map of greenhouse gas emissions, specifically including:

[0075] Acquired full grassland remote sensing images (including five multispectral bands and thermal infrared imagery) covering the entire monitoring area, and simultaneously acquired a digital elevation model (DEM) for the area. The spatial resolution of the DEM was consistent with the remote sensing images, both being 10 cm. Gradient calculations were performed on the DEM: the slope of each pixel was obtained by dividing the elevation difference between that pixel and its adjacent pixels by the horizontal distance; the aspect of each pixel was obtained by calculating the arctangent of the elevation change rates in the east and north directions, with values ​​ranging from 0 to 360 degrees. Slope values ​​were divided into intervals of 5 degrees (0 to 5 degrees, 5 to 10 degrees… up to 50 degrees and above), and aspect values ​​were divided into intervals of 45 degrees (0 to 45 degrees, 45 to 90 degrees… 315 to 360 degrees). Pixels with the same slope and aspect intervals and spatially adjacent locations on the full grassland remote sensing images were aggregated into a single patch. Isolated patches with an area less than 100 square meters were merged into the adjacent largest patch. The result is multiple homogeneous topographic patches, with pixels within each patch having similar topographic conditions, thus reducing the interference of topography on remote sensing reflectance.

[0076] For each homogeneous topographic patch, the reflectance values ​​of all pixels within the patch are extracted across five multispectral bands. The frequency distribution of reflectance values ​​for each pixel in each band is calculated, and the most frequently occurring reflectance value (the mode) is taken as the representative value for that band. If multiple modes exist for a band, the median is used. The representative values ​​from the five bands are combined into a five-dimensional vector, called the representative feature of the patch. This representative feature reflects the most typical spectral information of ground features within the patch.

[0077] The representative features of the polygons are substituted into the updated flux estimation rule. The updated flux estimation rule, derived from S4, contains five mapping coefficients, denoted as follows: , , , , This corresponds to five multispectral bands (blue, green, red, red-edge, and near-infrared). Let the reflectance values ​​of the five bands representing the characteristics be respectively... , , , , The initial flux The calculation formula is: ;in, The unit is milligrams per square meter per hour; to These are dimensionless mapping coefficients, whose values ​​are obtained in S4 through iterative convergence; to This is the reflectivity value, ranging from 0 to 1.

[0078] Retrieve the three paired samples with the smallest spectral angular distance to the representative feature of the current patch from the pulse-synchronous pairing dataset. The spectral angular distance is calculated as follows: for two five-dimensional vectors A and B, first calculate the dot product of A and B, then calculate the magnitude of A (the square root of the sum of the squares of each component) and the magnitude of B. Divide the dot product by the product of the two magnitudes to obtain the cosine of the angle between the two vectors. The inverse cosine of this cosine (angle) is the spectral angular distance, in degrees. The smaller the spectral angular distance, the more similar the two spectral vectors are. Extract the pixel reflectance vectors (each vector is also a five-dimensional reflectance value) corresponding to all synthetic remote sensing feature maps from the pulse-synchronous pairing dataset. Calculate the spectral angular distance between each vector and the representative feature, and select the three vectors with the smallest distance as paired samples. For each paired sample, extract its corresponding flux peak (denoted as ). , ) and the initial flux calculated in S4 using the updated rules for this sample (denoted as ) The calculation method is the same as (The calculation formula). Calculate the residual for each paired sample. Then calculate the average of these three residuals. Design a weighting coefficient that decreases as the spectral angular distance increases. Let the average spectral angular distance between the representative feature of the current patch and the three paired samples be... (Unit: degrees). Weighting coefficients The calculation formula is: ;in, It is a natural constant; The average spectral angular distance is expressed in degrees; 10 is the attenuation constant, meaning that for every 10-degree increase in spectral angular distance, the weighting coefficient decreases to approximately 0.368 times its original value. This weighting coefficient reduces the correction amount when the spectral angular distance between the paired sample and the representative feature is large, thus avoiding overcorrection.

[0079] Multiply the residual average by the weighting factor. The correction amount is obtained. The correction value is then added to the initial flux to obtain the final representative flux. For each pixel within the current patch, calculate the spectral angular distance between the pixel's five-band reflectance vector and the patch's representative feature (calculation method as above). Let the maximum spectral angular distance of all pixels within the patch be [value missing]. The minimum value is Spectral angular distance for a given pixel Its inverse distance weighting is: if equal If the weight is 1, then the weight is 1; otherwise, the weight is equal to 1. Divide by Multiply this weight by the value representing the flux. This yields the initial flux for that pixel. After performing the above assignment on all pixels, each pixel obtains a greenhouse gas flux value.

[0080] Create a blank raster with the same geographic extent and spatial resolution as the entire grassland remote sensing image, matching the number of rows and columns of pixels in the image. Fill the blank raster with the flux value of each pixel according to its row and column index. After filling, smooth the raster using median filtering: the filter window size is 3 pixels × 3 pixels. For each pixel, take the flux values ​​of itself and its eight neighboring pixels, sort them by value, and take the median value as the final flux for that pixel. After traversing all pixels, a smoothed flux raster is obtained. Add geographic coordinate reference information (consistent with the original remote sensing image) and color mapping (green to yellow for low flux, orange to red for high flux) to this raster, and output it as a GeoTIFF format file, i.e., a spatial distribution map of greenhouse gas emissions. This distribution map can be directly used to identify emission hotspots and guide grazing management decisions.

[0081] The working principle of this invention is as follows: A ground monitoring array is deployed on grazing natural grassland to continuously collect data on greenhouse gas fluxes, precipitation intensity, and soil temperature and humidity. Simultaneously, multi-temporal multispectral and thermal infrared images of the grassland are acquired, forming a ground flux sequence and a remote sensing image sequence. Precipitation events are identified based on precipitation intensity. The flux peak value in the first time period after precipitation is extracted from the ground flux sequence, and the first-temporal image before precipitation and the second-temporal image after precipitation are extracted from the remote sensing image sequence. Through temporal synthesis, a synthetic remote sensing feature map synchronous with the flux peak is generated using the first-temporal image, the second-temporal image, and the high-resolution flux sequence before and after the flux peak. The synthetic remote sensing feature map is associated and stored with the flux peak value to form a pulse-synchronized paired dataset. Based on the pulse-synchronized paired dataset, the mapping relationship between remote sensing features and flux in the baseline flux estimation rule is corrected to obtain an updated flux estimation rule. Finally, the updated flux estimation rule is applied to the entire grassland remote sensing image to calculate greenhouse gas flux pixel-by-pixel, generating a spatial distribution map of greenhouse gas emissions.

[0082] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for detecting greenhouse gas emissions from grazing natural grasslands, characterized in that, Includes the following steps: S1: Deploy ground monitoring arrays on grazing natural grasslands to continuously collect data on greenhouse gas fluxes, precipitation intensity, and soil temperature and humidity. At the same time, acquire multi-temporal multispectral and thermal infrared images of the grasslands to form ground flux sequences and remote sensing image sequences. S2: Identify precipitation events based on precipitation intensity, extract the peak flux in the first time period after precipitation from the ground flux sequence, and extract the first time phase image before precipitation and the second time phase image after precipitation from the remote sensing image sequence; through time-series synthesis, generate a synthetic remote sensing feature map that is in the same phase as the peak flux using the first time phase image, the second time phase image, and the high-resolution flux sequence before and after the peak flux. S3: Associate and store the synthetic remote sensing feature map with the flux peak to form a pulse synchronization paired dataset; S4: Based on the pulse synchronization paired dataset, the mapping relationship between remote sensing features and flux in the baseline flux estimation rule is corrected to obtain the updated flux estimation rule; S5: Apply the updated flux estimation rules to the whole grassland remote sensing image, calculate greenhouse gas flux pixel by pixel, and generate a spatial distribution map of greenhouse gas emissions.

2. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 1, characterized in that, The extraction of peak fluxes from the ground flux sequence during the first time period after precipitation, and the extraction of the first-phase image before precipitation and the second-phase image after precipitation from the remote sensing image sequence, specifically include: The soil moisture change rate after the falling edge of precipitation intensity is extracted from the surface flux sequence. When the change rate turns from positive to negative and the absolute value exceeds the standard deviation of historical background fluctuations, the flux value corresponding to the turning point is taken as the flux peak. The first phase image is the image frame with the smallest spatial variance of surface temperature in the most recent thermal infrared image before the precipitation begins. The second temporal image is the image frame corresponding to the moment when the slope of the vegetation red edge returns to the pre-precipitation level in the first multispectral image after the precipitation ends.

3. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 1, characterized in that, The process of generating the synthetic remote sensing feature map is as follows: The first and second temporal images were subjected to non-downsampled contour wave transforms to decompose the low-frequency subband coefficients and high-frequency directional subband coefficients at multiple scales. The flux ratio between the peak time and the previous time in the high-resolution flux sequence before and after the flux peak is used as the low-frequency weighting factor, and the flux ratio between the peak time and the next time is used as the directional weighting factor for the high-frequency directional subband coefficient. The low-frequency weighting factor is multiplied by the sum of the low-frequency sub-band coefficients of the first and second temporal images, and the directional weighting factor is multiplied by the difference between the high-frequency directional sub-band coefficients of the corresponding directions of the first and second temporal images. The product result is then subjected to non-downsampled contour wave inverse transform to reconstruct the synthetic remote sensing feature map.

4. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 3, characterized in that, The non-subsampled contour wave transformation specifically includes: Calculate the local grayscale variance maps of the first and second temporal images, use the absolute value of the difference between the two local grayscale variance maps as the weight map, and then determine the number of decomposition scale levels of the non-subsampled contour wave transform based on the energy distribution of the weight map. The number of directional sub-bands at each scale level is determined by the number of peaks in the polar coordinate sector after the Fourier transform of the weight map at the corresponding level. Based on the number of decomposition scale levels and the number of directional subbands at each level, non-subsampled Laplacian pyramid decomposition and non-subsampled directional filter bank decomposition are performed on each image to output low-frequency subband coefficients and high-frequency directional subband coefficients at multiple scales.

5. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 1, characterized in that, S3 specifically includes: A two-dimensional hash map is constructed using the spatial pixel locations of the synthesized remote sensing feature maps as row indices and the occurrence time of flux peaks as column indices. The reflectance value vector of each pixel in the synthetic remote sensing feature map and the flux peak value at the corresponding pixel location are stored as key-value pairs in a two-dimensional hash mapping table. Using the geographical coordinates and time window of the precipitation event as the query key, all associated key-value pairs are extracted and concatenated in chronological order to form a pulse synchronization paired dataset.

6. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 5, characterized in that, The construction of the two-dimensional hash mapping table specifically includes: For each spatial pixel location, calculate its geospatial grid code of latitude and longitude coordinates, and use the geospatial grid code as the row key value; For each flux peak occurrence time, it is converted into a continuous time sequence number from the reference time and mapped to a column key value using a prime number hash function; The row key and column key are concatenated to form a composite key. The reflectance vector of the synthetic remote sensing feature map at the corresponding pixel position and the corresponding flux peak are used as the composite value and stored in the hash bucket to complete the construction of the two-dimensional hash mapping table.

7. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 1, characterized in that, S4 specifically includes: Extract the pixel reflectance vector and its corresponding flux peak of each synthetic remote sensing feature map from the pulse synchronization pairing dataset, and calculate the residual between the flux peak and the predicted flux output by the baseline flux estimation rule. Using the Mahalanobis distance between pixel reflectance vectors in feature space as weights, the residuals are inversely distributed to the mapping coefficients of each feature component in the baseline rule according to the distance. This process is iterated until the sum of squared residuals of all paired datasets converges, thus obtaining the updated mapping relationship.

8. The detection method for greenhouse gas emissions from grazing natural grasslands according to claim 1, wherein The generation of the spatial distribution map of greenhouse gas emissions specifically includes: The whole grassland remote sensing image is divided into multiple topographically homogeneous patches according to the topographic slope and aspect. The mode vector of the pixel reflectance vector in each patch is used as the representative feature of the patch. The representative features are substituted into the updated flux estimation rules to calculate the representative flux of the patch. Then, the representative flux is linearly inversely distributed according to the spectral angular distance between the reflectance vector of each pixel and the representative features to obtain the greenhouse gas flux of each pixel. The flux values ​​of all pixels are filled into the blank raster according to the geographic coordinates of the original image, and after median filtering and smoothing, a spatial distribution map of greenhouse gas emissions is output.

9. The method for detecting greenhouse gas emissions from grazing natural grasslands according to claim 8, characterized in that, The calculation of the representative flux of the patch specifically includes: The reflectance values ​​of each band in the representative features are multiplied one by one with the mapping coefficients of the corresponding bands in the updated flux estimation rules, and then summed to obtain the preliminary flux. Then, retrieve the three paired samples with the smallest spectral angular distance to the representative feature from the pulse synchronization paired dataset, and calculate the average residual between the peak flux of the corresponding pixel of the synthesized remote sensing feature map in these three paired samples and the initial flux. The residual average is multiplied by a weighting coefficient that decreases as the spectral angular distance increases, and then superimposed onto the initial flux to obtain the final representative flux.

Images