Methane inversion method and device, electronic equipment and computer medium

Abstract

The present disclosure discloses a methane inversion method and device, electronic equipment and computer medium, which relates to the field of aerospace technology, and the specific implementation scheme is: acquiring multispectral satellite observation data; independently fitting the reflectivity model coefficients of the first wave band and the second wave band respectively, and processing the multispectral satellite observation data to obtain a reflectance distribution function corrected reflectivity; taking each pixel as the center in a local window, estimating the spectral interval scaling coefficient between the first wave band and the second wave band based on a robust loss function and an iterative reweighted least squares method, and obtaining a spectral interval scaling coefficient field; based on the spectral interval scaling coefficient field, calculating a target fractional absorption model of a target time phase and a reference fractional absorption model of a reference time phase respectively; calculating the difference between the two to obtain a differential absorption model; under the maximum a posteriori probability framework, based on a pre-constructed lookup table, performing a nonlinear optimization inversion on the differential absorption model to obtain a methane column concentration enhancement.

Description

Technical Field

[0001] This invention belongs to the field of aerospace technology, particularly the field of remote sensing monitoring technology, and specifically relates to a methane inversion method and apparatus, electronic equipment, and computer-readable storage medium. Background Technology

[0002] Methane (CH4) is the second largest greenhouse gas after carbon dioxide (CO2), with a global warming potential approximately 28 times that of CO2 over the next century. Human activities such as oil and gas extraction, coal mining, and waste disposal are the main sources of rising atmospheric methane concentrations. The United Nations Environment Programme (UNEP) estimates that about 40% of global methane emissions originate from oil and gas systems, with a significant portion of these emissions coming from a relatively small number of but extremely large "super-emission sources." Timely identification and location of these emission sources are crucial for developing effective emission reduction measures.

[0003] Satellite remote sensing technology has become an important tool for monitoring methane emissions due to its wide coverage and sustainable observation capabilities. In scenarios where multispectral satellite data is used to retrieve surface parameters, it is usually necessary to compare and analyze observation data acquired at different times in order to identify subtle signal differences caused by changes in atmospheric composition.

[0004] However, during this process, changes in observation geometry introduce systematic signal biases between different time phases, and the spatial heterogeneity of the Earth's surface also leads to complex spatial variation patterns in spectral relationships. Furthermore, the observational data inevitably contain outliers, and the presence of random noise poses a challenge to obtaining stable and reliable inversion results under low signal-to-noise ratio conditions.

[0005] Therefore, in the aforementioned remote sensing monitoring scenarios, there is a need for a technical solution that can effectively correct signal deviations caused by changes in observation conditions, suppress spatial heterogeneity interference and the influence of anomalous samples, and maintain inversion stability in noisy environments. Summary of the Invention

[0006] The purpose of this disclosure is to solve the technical problems of existing methane monitoring and to provide a methane inversion method and apparatus, electronic equipment and computer-readable storage medium.

[0007] The first aspect of this disclosure provides a methane inversion method, which includes: acquiring multispectral satellite observation data of the target time phase and a reference time phase, wherein the observation data includes at least the apparent reflectance of a first band and the apparent reflectance of a second band, which are sensitive to methane absorption; using the multi-temporal observation data, independently fitting the reflectance model coefficients of the first band and the second band, and using the reflectance model coefficients to uniformly correct the multispectral satellite observation data to a set reference geometric reference, thereby obtaining a reflectance corrected by the reflectance distribution function that eliminates observation geometric differences; and, for the reflectance corrected by the reflectance distribution function, within a local window centered on each pixel, based on stable... The inter-spectral scaling factor between the first and second bands is estimated using a robust loss function and iterative reweighted least squares method, resulting in a pixel-by-pixel spatially varying inter-spectral scaling factor field. Based on this field, the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase are calculated. The difference between the target fractional absorption model and the reference fractional absorption model is then calculated to obtain the differential absorption model. Under the maximum a posteriori probability framework, a nonlinear optimization inversion is performed on the differential absorption model based on a pre-constructed lookup table to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

[0008] A second aspect of this disclosure provides a methane inversion apparatus, comprising: an acquisition unit configured to acquire multispectral satellite observation data of a target time phase and a reference time phase, the observation data including at least the apparent reflectance of a first band and the apparent reflectance of a second band sensitive to methane absorption; a fitting unit configured to independently fit the reflectance model coefficients of the first and second bands using the multi-temporal observation data, and to uniformly correct the multispectral satellite observation data to a set reference geometric reference using the reflectance model coefficients to obtain a corrected reflectance based on the reflectance distribution function that eliminates observation geometric differences; and an estimation unit configured to correct the reflectance based on the reflectance distribution function, centered on each pixel within a local window. Within the system, the inter-spectral scaling coefficients between the first and second bands are estimated based on a robust loss function and iterative reweighted least squares method, resulting in a pixel-by-pixel spatially varying inter-spectral scaling coefficient field. The computation unit is configured to calculate the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase, respectively, based on the inter-spectral scaling coefficient field. The difference between the target fractional absorption model and the reference fractional absorption model is then calculated to obtain the differential absorption model. The inversion unit is configured to perform a nonlinear optimization inversion of the differential absorption model within a maximum a posteriori probability framework, based on a pre-constructed lookup table, to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

[0009] A third aspect of this disclosure provides a computer device including a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the method as described in the first aspect.

[0010] The fourth aspect of this disclosure provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method of the first aspect.

[0011] The methane inversion method and apparatus provided in the embodiments of this disclosure firstly acquire multispectral satellite observation data of the target time phase and a reference time phase. The observation data includes at least the apparent reflectance of a first band and the apparent reflectance of a second band, which are sensitive to methane absorption. Secondly, using the multi-temporal observation data, the reflectance model coefficients of the first and second bands are independently fitted respectively. The reflectance model coefficients are then used to uniformly correct the multispectral satellite observation data to a set reference geometric reference, resulting in a reflectance distribution function corrected reflectance that eliminates observation geometric differences. Thirdly, the reflectance distribution function corrected reflectance is calculated based on a local window centered on each pixel, using a stable... The inter-spectral scaling factor between the first and second bands is estimated using a robust loss function and iterative reweighted least squares method, resulting in a pixel-by-pixel spatially varying inter-spectral scaling factor field. Then, based on this field, the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase are calculated. The difference between the target and reference fractional absorption models is calculated to obtain the differential absorption model. Finally, within the maximum a posteriori probability framework, a nonlinear optimization inversion is performed on the differential absorption model based on a pre-constructed lookup table to obtain the methane column concentration enhancement. The lookup table characterizes the physical mapping relationship between remote sensing observations and methane concentration. Therefore, in the process of multispectral satellite methane inversion, the reflectivity model coefficients of the two bands can be independently fitted separately and the observation data can be uniformly corrected to the reference geometric benchmark, effectively eliminating the systematic signal bias introduced by the difference in observation geometry between different time phases; by estimating the inter-spectral scaling coefficients based on a robust loss function and iterative reweighted least squares method within a local window with each pixel as the center, it can adapt to local surface spectral characteristics, suppress the interference of spatial heterogeneity on spectral relationships, and reduce the impact of anomalous samples on the estimation results; by performing a nonlinear optimization inversion based on a lookup table under the maximum a posteriori probability framework, it can adaptively balance the prior constraints and observation data, thereby maintaining the stability of the inversion solution in a noisy environment, meeting the technical requirements for signal bias correction, heterogeneity interference suppression, and inversion stability in remote sensing monitoring scenarios. Attached Figure Description

[0012] Figure 1 This is a flowchart of an embodiment of the methane inversion method according to the present disclosure; Figure 2 This is a flowchart of another embodiment of the methane inversion method according to the present disclosure; Figure 3 This is a schematic diagram of the structure of one embodiment of the methane inversion apparatus according to the present disclosure; Figure 4 This is a block diagram of an electronic device used to implement the methane inversion method of the embodiments of this disclosure. Detailed Implementation

[0013] Unless otherwise expressly stated, throughout the specification and claims, the term "comprising" or its variations such as "including" or "comprises" shall be understood to include the stated elements or components without excluding other elements or other components.

[0014] The technical solution of the present invention is illustrated below through specific embodiments. It should be understood that the one or more steps mentioned in the present invention do not preclude the existence of other methods and steps before or after the combined steps, or that other methods and steps may be inserted between these explicitly mentioned steps. It should also be understood that these examples are for illustrative purposes only and are not intended to limit the scope of the present invention. Unless otherwise stated, the numbering of each method step is only for the purpose of identifying each method step, and not for limiting the order of each method or limiting the scope of the present invention. Changes or adjustments to their relative relationships, without substantial changes to the technical content, can also be considered as within the scope of the present invention.

[0015] The raw materials and instruments used in the examples are not subject to any specific restrictions on their source; they can be purchased from the market or prepared according to conventional methods known to those skilled in the art.

[0016] In satellite remote sensing monitoring scenarios, multispectral sensors are typically used to acquire surface reflectance information. By comparing and analyzing observation data from different bands and time phases, weak absorption signals caused by changes in atmospheric composition can be identified, and the column concentration enhancement of the target gas can then be retrieved. Taking methane monitoring as an example, the shortwave infrared bands of multispectral satellites have different sensitivities to methane absorption. By comparing the apparent reflectance of the first and second bands, which are sensitive to methane absorption, as well as the difference between the target time phase and the reference time phase, the methane enhancement signal can be extracted. However, in this process, changes in observation geometry, surface spatial heterogeneity, and the presence of random noise can all significantly affect the inversion accuracy.

[0017] Some processing methods directly use the raw apparent reflectance to calculate the ratios or differences between bands or time phases to eliminate the influence of surface background. While this method is computationally simple, it fails to adequately consider the systematic bias in reflectance caused by observational geometric differences between different time phases, and it does not adaptively correct for spatial variations in surface spectral relationships. When the observation scene contains high heterogeneity and the surface state changes between time phases, the direct difference results are prone to retaining artifacts related to surface changes. These artifacts may have spatial morphology and spectral characteristics highly similar to the actual plume, leading to misidentification.

[0018] In other processing methods, a global constant is used to describe the spectral scaling relationship between different bands, and multi-temporal statistical averaging is used to suppress random noise. This approach can achieve certain results under homogeneous surface conditions, but in heterogeneous and complex scenarios, the global constant cannot adapt to the local spectral feature differences of different land cover types, leading to systematic errors at the boundaries. In addition, the sample selection for statistical averaging lacks optimization for noise features, has insufficient robustness to anomalous samples, and has limited noise suppression efficiency.

[0019] In view of this, the embodiments of this specification provide a methane inversion method, which establishes a cascaded processing chain of geometric correction and spectral adaptive correction, and introduces physical forward model constraints and adaptive regularization under the maximum a posteriori probability framework to achieve comprehensive suppression of observation geometric bias, surface heterogeneity interference and random noise. Figure 1 A flowchart 100 illustrating an embodiment of a methane inversion method is shown, which includes the following steps: Step 101: Obtain multispectral satellite observation data for the target time phase and the reference time phase.

[0020] In this embodiment, the observation data includes at least the apparent reflectance of a first band and the apparent reflectance of a second band, which are sensitive to methane absorption. The first and second bands can be shortwave infrared bands, wherein the center wavelength of the first band is located around 1610 nm, and the center wavelength of the second band is located around 2190 nm.

[0021] In this embodiment, "time phase" is a technical term in the field of remote sensing, referring to the specific time or date when a satellite passes over the target area for observation. The target time phase is the observation time of the methane emission to be detected, and the reference time phase is the clean observation time assuming no methane plume.

[0022] In this embodiment, apparent reflectance is atmospheric top reflectance, which refers to the reflectance converted from the radiation directly received by the satellite sensor. It includes the combined results of surface reflection information, atmospheric scattering (path radiation), and gas absorption (such as methane absorption), and has not yet undergone atmospheric and geometric correction.

[0023] In this embodiment, the observation data includes at least the apparent reflectance of the first band and the apparent reflectance of the second band. The first band is sensitive to methane absorption, and the second band is also sensitive to methane absorption, but with a different sensitivity than the first band. Furthermore, the observation data may also include observational geometric parameters such as the solar zenith angle, observation zenith angle, and relative azimuth angle for each pixel, used for subsequent reflectance model coefficient estimation and geometric correction. After acquiring the observation data, preprocessing can be performed, such as cloud masking, geometric registration, and radiometric calibration, to ensure that the target time phase and the reference time phase are comparable in spatial location and radiometric measure.

[0024] Step 102: Using multi-temporal observation data, independently fit the reflectivity model coefficients of the first band and the second band respectively, and use the reflectivity model coefficients to uniformly correct the multispectral satellite observation data to the set reference geometric benchmark, so as to obtain the corrected reflectivity of the reflectivity distribution function that eliminates the observation geometric differences.

[0025] In this embodiment, because the solar angle (where it shines) and satellite perspective (from which it looks) differ each time the satellite observes, the same ground feature may appear to have different brightness. The reflectance distribution function corrected reflectance refers to the pure surface reflectance normalized to a unified angular reference after eliminating the geometric differences in observation due to the solar angle and satellite perspective.

[0026] In this embodiment, reflectance model coefficient estimation and geometric correction are performed independently for the first and second bands, ensuring that while preserving the differences in their respective surface directional reflectance characteristics, the two bands are reclassified to the same geometric reference, thereby eliminating systematic reflectance deviations caused by changes in observed geometry between different time phases. Because the coefficients are estimated and corrected independently for the two bands, both the physical differences between the bands are preserved, and the two corrected bands are ensured to be comparable on the same geometric reference.

[0027] In this embodiment, the multi-temporal observation data can be multiple historical observation data of the target area within a set time window. Step 102 includes: assigning time decay weights to each historical observation data according to the time interval between the acquisition time of each historical observation data and the target temporal phase, wherein the smaller the time interval, the larger the assigned time decay weight; using the multi-temporal observation data with assigned time decay weights, fitting reflectance model coefficients pixel by pixel using the weighted least squares method. Using the reflectance model coefficients, calculating the first reflectance distribution function value of the pixel under the current actual observation geometry and the second reflectance distribution function value under the reference geometric reference; using the ratio of the second reflectance distribution function value to the first reflectance distribution function value, performing multiplicative correction on the apparent reflectance of the first band and the second band respectively, to obtain the reflectance distribution function corrected reflectance that eliminates the observation geometry difference.

[0028] Step 103: Correct the reflectance of the reflectance distribution function. With each pixel as the center, estimate the inter-spectral scaling coefficient between the first band and the second band within the local window based on the robust loss function and the iterative reweighted least squares method, and obtain the inter-spectral scaling coefficient field that varies spatially for each pixel.

[0029] In this embodiment, the interspectral scaling factor refers to the proportional relationship between the first band and the second band that are sensitive to methane. Since the surface is not completely uniform, the ratio of these two bands is different for different ground features. The interspectral scaling factor field indicates that the interspectral scaling factor is not a globally unique constant, but a two-dimensional matrix or image that varies per pixel space.

[0030] In this embodiment, surface spectral inhomogeneity is further eliminated based on the geometrically corrected reflectance. A local window is constructed centered on each pixel, and all valid pixels within the window constitute the sample set for that location, used to estimate the inter-spectral scaling factor at that pixel. By estimating through a pixel-by-pixel sliding window, the spectral scaling factor of each pixel adapts to the local surface spectral characteristics, upgrading the global constant description to a spatially adaptive field, thereby suppressing systematic errors at the boundaries of different land cover types in heterogeneous scenes.

[0031] In one embodiment of this disclosure, within a local window, the relationship between the first band and the second band is represented by a linear model including an intercept; an iterative reweighted least squares method is used to simultaneously estimate the interspectral scaling factor and the intercept term; wherein the intercept term is used to absorb the effects of atmospheric path radiation and surface offset, and the intercept term is eliminated in the two temporal differences and does not contribute to the differential absorption model.

[0032] Optionally, before performing the iterative reweighted least squares method, the methane inversion method further includes: identifying effective pixels within a local window and assigning them effective pixel weights, which are used to remove invalid observation data; using the effective pixels assigned effective pixel weights, calculating the initial interspectral scaling factor through the weighted least squares method, and using the initial interspectral scaling factor as the iteration starting point of the iterative reweighted least squares method.

[0033] Step 104: Based on the interspectral scaling coefficient field, calculate the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase respectively; calculate the difference between the target fractional absorption model and the reference fractional absorption model to obtain the differential absorption model.

[0034] In this embodiment, the fractional absorption model refers to the ratio / difference between the theoretical value of the second band predicted using the first band and the actual observed value of the second band in a single observation. The differential absorption model refers to subtracting the fractional absorption model of the target time phase from the fractional absorption model of the reference time phase (taking the difference). The differential absorption model is the core mathematical structure for completely eliminating surface heterogeneity artifacts.

[0035] This disclosure supports two model forms: the no-intercept term model. Linear scaling through the origin; model with intercept term Additional intercept term included It can absorb the effects of atmospheric path radiation and surface offset. Among them Indicates the reflectivity of the first band, B11. This represents the reflectivity of the second band, B12. The model with the intercept term is the default; the fractional absorption model in two forms is also given below.

[0036] In this embodiment, the fractional absorption model is defined as shown in equations (1) and (2): No intercept term model : (1) Model with intercept term Intercept term Does not change the sign of the methane signal: (2) In this embodiment, the inter-spectral scaling coefficient field estimated in step 103 is used to calculate fractional absorption models for both the target and reference time phases. The fractional absorption model can be expressed as the difference between the corrected first-band reflectance and the inter-spectral scaling second-band reflectance, divided by the first-band reflectance. This model design allows the inter-spectral scaling coefficient to act as a spectral normalization factor, normalizing the local surface spectral relationship to the reference spectral baseline. After calculating the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase, the difference between the two is obtained to obtain the differential absorption model. This differential operation further eliminates the surface background contribution between time phases, retaining only the differential signal caused by changes in atmospheric composition. For pixels where the inter-spectral scaling coefficient is invalid, a preset constant can be used as a backoff value to ensure the continuity of the calculation process. Furthermore, when selecting the reference time phase, observations with a season similar to the target time phase can be prioritized to reduce false positives caused by changes in surface conditions.

[0037] Step 105: Under the maximum a posteriori probability framework, based on the pre-built lookup table, perform a nonlinear optimization inversion on the differential absorption model to obtain the methane column concentration enhancement.

[0038] In this embodiment, a lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

[0039] In this embodiment, step 105 includes: in each iteration of the nonlinear optimization, using the currently assumed methane column concentration enhancement as an input parameter; using the input parameter to interpolate the lookup table to obtain the physical forward simulation value corresponding to the currently assumed concentration; calculating the observation residual between the differential absorption model and the physical forward simulation value, and updating the assumed methane column concentration enhancement based on the observation residual in conjunction with prior constraints, until the set convergence condition is met.

[0040] Optionally, the above methane inversion method further includes: setting a prior probability distribution for the methane column concentration enhancement to be inverted within the maximum a posteriori probability framework; wherein the prior mean of the prior probability distribution is set to zero to characterize the physical background state of the target area without point source emissions; during the nonlinear optimization inversion process, the prior mean set to zero is used as a regularization constraint to suppress non-physical extreme inversion values ​​at low signal-to-noise ratio pixels, thereby avoiding false positives of methane in the inversion results.

[0041] In this embodiment, the differential absorption model is inverted into the methane column concentration enhancement. A pre-constructed lookup table simulates the physical mapping relationship between the differential absorption model and the methane column concentration enhancement under different observation conditions using a radiative transfer model, thereby transforming the complex physical forward computation into an efficient lookup table interpolation operation. Within the maximum a posteriori probability framework, a cost function is constructed that includes observation terms and prior terms. The observation term measures the difference between the observed values ​​of the differential absorption model and the physical forward simulation values ​​obtained through lookup table interpolation, and is normalized using the observation noise standard deviation. The prior term measures the difference between the methane column concentration enhancement to be inverted and the prior mean, and is normalized using the prior standard deviation. By balancing the observation and prior terms, the observed data can be fully trusted under high signal-to-noise ratio conditions, while prior constraints can be used to avoid solution divergence under low signal-to-noise ratio conditions.

[0042] The methane inversion method provided in this disclosure firstly acquires multispectral satellite observation data for the target and reference time phases. The observation data includes at least the apparent reflectance of a first band and a second band, both sensitive to methane absorption. Secondly, using the multi-temporal observation data, reflectance model coefficients for the first and second bands are independently fitted, and these coefficients are used to uniformly correct the multispectral satellite observation data to a set reference geometric benchmark, resulting in a reflectance distribution function-corrected reflectance that eliminates observational geometric differences. Thirdly, the reflectance distribution function-corrected reflectance is applied to each pixel within a local window, based on robust loss... The inter-spectral scaling factor between the first and second bands is estimated using the loss function and iterative reweighted least squares method, resulting in a pixel-by-pixel spatially varying inter-spectral scaling factor field. Then, based on the inter-spectral scaling factor field, the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase are calculated respectively. The difference between the target fractional absorption model and the reference fractional absorption model is obtained to obtain the differential absorption model. Finally, under the maximum a posteriori probability framework, a nonlinear optimization inversion is performed on the differential absorption model based on a pre-constructed lookup table to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration. Therefore, in the process of multispectral satellite methane inversion, the reflectivity model coefficients of the two bands can be independently fitted separately and the observation data can be uniformly corrected to the reference geometric benchmark, effectively eliminating the systematic signal bias introduced by the difference in observation geometry between different time phases; by estimating the inter-spectral scaling coefficients based on a robust loss function and iterative reweighted least squares method within a local window with each pixel as the center, it can adapt to local surface spectral characteristics, suppress the interference of spatial heterogeneity on spectral relationships, and reduce the impact of anomalous samples on the estimation results; by performing a nonlinear optimization inversion based on a lookup table under the maximum a posteriori probability framework, it can adaptively balance the prior constraints and observation data, thereby maintaining the stability of the inversion solution in a noisy environment, meeting the technical requirements for signal bias correction, heterogeneity interference suppression, and inversion stability in remote sensing monitoring scenarios.

[0043] Optionally, when calculating the differential absorption model, if there is a change in the surface state between the target time phase and the reference time phase, the above methane inversion method further includes: independently estimating the inter-spectral scaling factor of the reference time phase; and performing a trade-off or smoothing regularization process on the inter-spectral scaling factor of the target time phase and the inter-spectral scaling factor of the reference phase to suppress false positive methane signals caused by abrupt changes in the surface state.

[0044] In some optional implementations of this disclosure, the above-mentioned use of multi-temporal observation data to independently fit the reflectivity model coefficients of the first and second bands, and using the reflectivity model coefficients to uniformly correct the multispectral satellite observation data to a set reference geometric reference, thereby obtaining the reflectivity correction function that eliminates observation geometric differences, includes: using multi-temporal observation data, estimating the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient of the first and second bands pixel by pixel through least squares fitting; based on the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient, using the arithmetic mean of the scene angles of the target temporal phase and the reference temporal phase as the reference geometric reference for reduction, thereby obtaining the reflectivity correction function that eliminates observation geometric differences.

[0045] In this optional implementation, the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient are the three core parameters of the Ross-Li kernel-driven bidirectional reflectance distribution function. Isotropic scattering assumes that the Earth's surface reflects light uniformly in all directions, like a Lambertian object. Volume scattering refers to light that has been scattered multiple times after entering a medium such as vegetation canopy before being emitted again. Geometric scattering is the optical occlusion effect caused by surface roughness (such as the shadows of trees and buildings).

[0046] In this optional implementation, the volume scattering kernel corresponding to the volume scattering coefficient and the geometric scattering kernel corresponding to the geometric scattering coefficient are calculated based on the solar zenith angle, the observation zenith angle, and the relative azimuth angle during satellite observation.

[0047] In this optional implementation, the Ross-Li nuclear-driven reflectivity model (hereinafter referred to as the BRDF model) is used to correct the reflectivity of the first band B11 and the second band B12 of the target time phase and the reference time phase, respectively. Specifically, it includes the following sub-steps: Note: The Ross-Li kernel-driven BRDF model is a very important semi-empirical, linear kernel-driven model in the field of remote sensing. It is mainly used to describe and simulate the bidirectional reflectance characteristics of the Earth's surface under different solar and / or observation geometry conditions.

[0048] use The BRDF model coefficients for each pixel were estimated using least-squares fitting based on observation data from different time phases. The Ross-Li kernel-driven BRDF model is represented as shown in equation (3): (3) In equation (3), Indicates apparent reflectance, Represents the isotropic scattering coefficient. Represents the volume scattering coefficient. Represents the geometric scattering coefficient. Relative azimuth angle. Volume scattering nucleus. and geometric scattering kernel Observational geometry calculations are performed using equations (4) and (5) respectively: (4) (5) In equation (4), The zenith angle of the sun. To observe the zenith angle, Let be the phase angle, satisfying In equation (5), , These are the projected angles of the sun and the observed zenith angle, respectively. This is the projection phase angle.

[0049] The coefficients are estimated using the least squares solution as shown in equation (6): (6) In equation (6), the matrix Observation vector Let N be the reflectance vector of effective pixels. Two sets of BRDF model coefficients are obtained by independently estimating the first band B11 and the second band B12 respectively. , , and , , .

[0050] For a given observation geometry ( , , First, calculate the BRDF value under this geometry. The reflectance is corrected by obtaining the reflectance distribution function through equation (7). : (7) In equation (7), The BRDF value under the reference geometry is obtained by taking the average scene angles of the target and reference phases. The arithmetic mean is used as a unified reference geometry to ensure that the two scenes are reduced to the same benchmark.

[0051] In this optional implementation, at least five different temporal observation data points can be used to estimate the bidirectional reflectance distribution function model coefficients for each pixel in each band through least-squares fitting. This model can employ a semi-empirical linear kernel-driven model, expressing apparent reflectance as a linear combination of isotropic scattering, volume scattering, and geometric scattering components, where the volume scattering kernel and geometric scattering kernel are calculated from the observed geometric parameters. Through least-squares solving, the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient for each pixel in each band can be obtained. After estimating the model coefficients, the arithmetic mean of the scene angles for the target and reference temporal phases is calculated as a unified reference geometric benchmark. For each pixel in each temporal phase, the current bidirectional reflectance distribution function value is calculated using its estimated model coefficients and the current observed geometry. Simultaneously, the bidirectional reflectance distribution function value under the reference geometric benchmark is calculated. Multiplying the original apparent reflectance by the ratio of the function value under the reference geometry to the function value under the current geometry yields the corrected reflectance distribution function. After this correction, the difference in reflectance of the same ground feature at different time phases mainly comes from actual changes in the ground surface and changes in atmospheric composition, rather than differences in observed geometry.

[0052] In some optional implementations of this disclosure, the above-mentioned method of correcting reflectance for the reflectance distribution function and estimating the inter-spectral scaling coefficient between the first and second bands within a local window centered on each pixel, based on a robust loss function and iterative reweighted least squares, includes: correcting reflectance based on the reflectance distribution function and calculating the residual of the current iteration; calculating a robust scaling estimate of the residual using the median absolute deviation; calculating the Hubel weight of each sample within the local window by combining a set Hubel threshold and the robust scaling estimate, wherein, for abnormal samples with residuals greater than the set threshold, their Hubel weights are reduced; and updating the weighted least squares using the Hubel weights to obtain new inter-spectral scaling coefficients, until the iterative convergence condition is met or the maximum number of iterations is reached.

[0053] This disclosure proposes a Local Adaptive Spectral Scaling (LASS) method, which estimates the local spectral scaling coefficients through a pixel-by-pixel sliding window to adaptively address the surface feature spectral relationship between the first band B11 and the second band B12. LASS supports two model forms: a model without an intercept term. Linear scaling through the origin; model with intercept term Additional intercept term included It can absorb the effects of atmospheric path radiation and surface offset. The model with the intercept term is the default; the formulas for each sub-step are given below for both forms. For each cell location Based on the position of that pixel Centered ( (If odd, default 51) Valid pixels within the neighborhood window are used as samples to estimate the inter-spectral scaling factors. Initial estimation uses ordinary weighted least squares: During initialization, the interspectral scaling coefficients of the model without intercept terms are shown in equation (8): (8) The interspectral scaling factor for the model with intercept term is given by equation (9): (9) In equations (8) and (9), Scaling factor between terms without intercept initial value, Refers to the pixel coordinates of value, and These are the scaling factors for the model with intercept term. and The initial values, where Refers to the pixel coordinates of value, Refers to the pixel coordinates of value, For Centered window, The effective cell weight (initialized to 1, and the invalid cell weight to 0). and The reflectance correction results for the first band B11 and the second band B12 in the reflectance correction function are shown respectively. To correct the reflectance of the first band B11 reflectance in the reflectance correction result of the reflectance distribution function, the pixel value in the i-th row and j-th column is used. The pixel value in the i-th row and j-th column of the second band B12 reflectance correction result for the reflectance distribution function correction.

[0054] Huber-IRLS is a robust regression algorithm that combines the Huber loss function with iteratively reweighted least squares. It uses Huber-IRLS for robust estimation of outliers that may exist within the window (such as incompletely removed cloud remnants, water body boundaries, methane plume edges, etc.) to improve the robustness of the inter-spectral scaling factor estimation. The iterative process is as follows: In the k-th iteration, the current residual is calculated first: Residuals of a model without intercept term As shown in equation (10): (10) Residuals of a model with intercept term As shown in equation (11): (11) Next, robust scaling estimates of the residuals are calculated using the Median Absolute Deviation (MAD) scaling estimator as shown in Equation (12): (12) Set Hubel threshold Calculate Hubel weights As shown in equation (13): (13) Weighted least squares update using Hubel weights: In the k-th iteration, the interspectral scaling coefficients of the model without intercept term As shown in equation (14): (14) Interspectral scaling factors of models with intercept terms As shown in equation (15): (15) In equations (14) and (15), The scaling factor between the spectrum of the term without intercept in the k-th iteration. The value, and The scaling factor of the model with intercept term in the k-th iteration. and The value, This represents the weight value in the i-th row and j-th column of the effective cell weights, with each summation term having a range of [value missing]. Repeat the above process until convergence. When the maximum number of iterations is reached (default 20), the Hubel weighting function switches from quadratic growth to linear growth at large residuals, reducing the impact of outliers on the estimation of interspectral scaling factors, while maintaining the estimation efficiency when the residuals are small.

[0055] In this optional implementation, in each iteration, the residuals of each sample within the window are first calculated based on the current inter-spectral scaling coefficients. Then, a robust scaling estimate of the residuals is calculated using the median absolute deviation, and a Hubel threshold is set accordingly. For samples whose absolute residual values ​​do not exceed the threshold, unit weights are assigned to maintain Gaussian optimality; for samples whose absolute residual values ​​exceed the threshold, the weights are decayed by the reciprocal of their absolute residual values, thereby reducing the impact of outliers on the estimation results. Weighted least squares are then performed using the updated weights to obtain new inter-spectral scaling coefficients. The above process is repeated until the difference between the coefficients obtained from two adjacent iterations is less than the convergence tolerance or the preset maximum number of iterations is reached. Since the Hubel weight function switches from quadratic loss to linear loss at large residuals, it is equivalent to automatically identifying outliers and adaptively reducing their weights in the estimation, eliminating the need for a preset hard threshold for outliers. Even under harsh conditions with a high proportion of outliers, the accuracy and robustness of the estimation can still be maintained.

[0056] Optionally, the interspectral scaling factor between the first and second bands can also be estimated using a linear model with an intercept term. The intercept term is used to absorb the effects of atmospheric path radiation and surface offset, and the intercept term is eliminated in the two time-phase differences and does not contribute to the differential absorption model.

[0057] In some embodiments of this disclosure, the above-described methane inversion method further includes: during the weighted least squares update process, when performing weighted summation calculations for samples within a local window, obtaining the weighted summation result of samples within the local window based on a pre-constructed integral map of the entire image, so that the complexity of the weighted summation calculation is independent of the size of the local window.

[0058] In this embodiment, an integral map can be pre-constructed for the array to be summed. This integral map contains the cumulative sum information of the array, allowing the sum of any rectangular sub-region to be obtained in constant time via a four-corner lookup table. By applying the integral map mechanism to the numerator and denominator summation terms in weighted least squares updates, the computational complexity of the weighted summation of all windows in each iteration can be reduced from being related to the window area to being linearly related to the image size. This significantly improves computational efficiency while maintaining the accuracy of large window estimation, making pixel-by-pixel large window sliding estimation feasible in engineering.

[0059] In some embodiments of this disclosure, the above-described methane inversion method further includes: determining whether the proportion of effective pixels within a local window is lower than a preset threshold; if it is lower than the preset threshold, using the effective estimation region as a seed, an iterative filling strategy based on morphological dilation is adopted to expand outwards successively, and the mean of the outer window is used to initialize the inter-spectral scaling coefficient of the invalid region during each expansion, until all pixels obtain effective inter-spectral scaling coefficients.

[0060] When the proportion of effective pixels within a window is lower than a preset threshold (such as in densely cloud-covered areas), directly estimating the inter-spectral scaling factor value is unstable. An iterative filling strategy based on morphological dilation is adopted: using the effective estimation region as a seed, it expands outwards successively, initializing the inter-spectral scaling factor of the invalid region with the mean of the large window each time, iterating until all pixels have effective values.

[0061] In this embodiment, the preset threshold can be set to approximately 30% depending on the application scenario. When the proportion of effective pixels within the window is lower than this threshold due to dense cloud cover or missing data, directly estimating the inter-spectral scaling factor may be unstable. In this case, using neighboring effective estimation regions as seeds, morphological dilation is performed outwards sequentially. Each time dilation occurs, the inter-spectral scaling factor of the newly included region is initialized with the mean of the outer effective window. This process is iterated until all pixels obtain effective values. This filling strategy can provide effective inter-spectral scaling factors for all pixels in areas with missing data, ensuring the spatial continuity of subsequent inversion steps and avoiding gaps in the inversion results due to local missing data.

[0062] In some alternative implementations of this disclosure, the above-mentioned nonlinear optimization inversion of the differential absorption model based on a pre-constructed lookup table, within the maximum a posteriori probability framework, yields the methane column concentration enhancement. This includes: constructing a cost function containing observation terms and prior terms, wherein the observation terms are constructed based on the difference between the differential absorption model and the physical forward simulation values, as well as the standard deviation of the observation noise; the physical forward simulation values ​​are obtained by interpolating the lookup table; the prior terms are constructed based on the difference between the methane column concentration enhancement to be inverted and the prior mean, as well as the prior standard deviation; and solving the cost function using the Levenberg-Marquardt algorithm, switching between the steepest descent method and the Gauss-Newton method using an adaptive damping factor to obtain the methane column concentration enhancement that minimizes the cost function.

[0063] In this optional implementation, the cost function is specifically expressed as the sum of the ratio of the square of the observation residual to the variance of the observation noise and the ratio of the square of the prior residual to the variance of the prior noise; where the observation residual is the difference between the differential absorption model and the physical forward simulation value, and the prior residual is the difference between the methane column concentration enhancement to be inverted and the prior mean.

[0064] In this optional implementation, physical forward simulation refers to using an atmospheric radiative transfer model (such as MODTRAN), inputting an assumed methane concentration, and using physical formulas to calculate the optical signal that the satellite should receive.

[0065] In this optional implementation, because the forward physics model calculation is extremely slow, in engineering practice, the satellite signals corresponding to various methane concentrations are usually calculated in advance and stored in a lookup table. That is, the lookup table sets up a table that corresponds various methane concentrations with the corresponding satellite signals. During inversion, the dictionary interpolation value of the lookup table is directly looked up, which greatly improves the calculation speed.

[0066] In this optional implementation, the differential absorption model Inversion into methane column concentration increment (Unit: ppm). The Levenberg-Marquardt (LM) algorithm is adopted, that is, under the maximum a posteriori probability framework, the LM algorithm is used to solve such a MAP estimation problem. The LM algorithm is an optimization algorithm specifically used to solve nonlinear least squares problems. The cost function is shown in Equation (16): (16) In equation (16), Obtained by interpolation using a lookup table pre-built using a radiative transfer model; Differential absorption model The standard deviation of the observation noise; This is a hypothetical value (without methane enhancement). The prior standard deviation is used. This is achieved through an adaptive damping factor. Switching between the Gauss-Newton method and the steepest descent method ensures global convergence.

[0067] In this optional implementation, the Levenberg-Marquardt algorithm introduces an adaptive damping factor to automatically adjust the search direction during iteration based on the degree of local linearization: when the residual is large or the Jacobian matrix is ​​ill-conditioned, it favors the steepest descent direction to ensure global convergence; when it is close to the optimal solution and the residual is small, it favors the Gauss-Newton direction to accelerate convergence. This adaptive mechanism effectively avoids gradient descent oscillations and Gauss-Newton divergence problems, ensuring the numerical stability of the inversion process under complex terrain and noisy conditions.

[0068] In some optional implementations of this disclosure, the strategy for determining the aforementioned prior standard deviation includes an empirical Bayesian strategy or a data-driven adaptive strategy: The empirical Bayesian strategy includes scanning candidate prior standard deviations on a background image without methane plumes and selecting the candidate value that minimizes the background noise as the prior standard deviation; The data-driven adaptive strategy includes dynamically adjusting the prior standard deviation according to the magnitude of the observation noise standard deviation. When the observation noise standard deviation increases, the prior standard deviation is reduced to enhance the prior regularization constraint, thereby avoiding divergence in the inversion solution under low signal-to-noise ratio conditions.

[0069] In this optional implementation, the regularization intensity can be automatically adapted to changes in observation conditions through data-driven adaptive adjustment of the prior standard deviation. This maintains low variance in methane-free regions to avoid false alarms and high sensitivity in real plume regions to avoid missed alarms.

[0070] Figure 2The illustrated embodiment is in Figure 1 Based on the steps shown, the sub-processes within each step and their data flow relationships are further illustrated. Figure 2 In the illustrated embodiment, after acquiring multispectral satellite observation data, the multi-temporal data (target temporal data, reference temporal data, and adjacent temporal data) from the multispectral satellite observation data are used for pixel-by-pixel least squares fitting to estimate the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient for each band. A unified reference geometric benchmark is constructed based on the arithmetic mean of the scene angles of the target and reference temporal phases, completing the reflection distribution function correction. The corrected reflectance enters the local adaptive interspectral scaling factor estimation stage, which may include initial weighted least squares estimation, Hubell iterative reweighted robust iterative update, constant complexity weighted summation calculation based on integral image, and morphological dilation iterative filling triggered when the effective pixel ratio is insufficient. After iterative convergence, a complete interspectral scaling factor field is obtained, based on which the fractional absorption models of the target and reference temporal phases are calculated and the difference is calculated to output the differential absorption model. Subsequently, physical forward simulation values ​​are obtained by interpolation based on a pre-built lookup table. A cost function containing observation and prior terms is constructed by combining the observation noise standard deviation and the adaptively adjusted prior standard deviation. The Levenberg-Marquardt algorithm is used to optimize the solution by switching between the Gauss-Newton method and the steepest descent method with an adaptive damping factor, and the methane column concentration enhancement is output.

[0071] Figure 2 The embodiment shown forms a complete closed loop from data input to result output by embedding sub-processing mechanisms within each major step, which can comprehensively suppress geometric deviations, heterogeneous interferences and random noise under complex surface conditions.

[0072] Further reference Figure 3 As an implementation of the methods shown in the above figures, this disclosure provides an embodiment of a methane inversion apparatus, which is similar to... Figure 1 Corresponding to the method embodiments shown, this device can be specifically applied to various electronic devices.

[0073] like Figure 3As shown, the methane inversion device 300 provided in this embodiment includes: an acquisition unit 301, a fitting unit 302, an estimation unit 303, a calculation unit 304, and an inversion unit 305. The acquisition unit 301 can be configured to acquire multispectral satellite observation data of the target time phase and a reference time phase. The observation data includes at least the apparent reflectance of a first band and the apparent reflectance of a second band, which are sensitive to methane absorption. The fitting unit 302 can be configured to independently fit the reflectance model coefficients of the first and second bands using the multi-temporal observation data, and then use the reflectance model coefficients to uniformly correct the multispectral satellite observation data to a set reference geometric benchmark, obtaining a reflectance correction reflectance that eliminates observation geometric differences. The estimation unit 303 can be configured to correct the reflectance using the reflectance distribution function, and within a local window centered on each pixel, estimate the interspectral scaling coefficients between the first and second bands based on a robust loss function and iterative reweighted least squares method, obtaining a pixel-by-pixel spatially varying interspectral scaling coefficient field. The aforementioned calculation unit 304 can be configured to calculate the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase based on the inter-spectral scaling coefficient field; and to obtain the differential absorption model by subtracting the target fractional absorption model and the reference fractional absorption model. The aforementioned inversion unit 305 can be configured to perform a nonlinear optimization inversion on the differential absorption model based on a pre-constructed lookup table within the maximum a posteriori probability framework to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

[0074] In this embodiment, the specific processing and technical effects of the acquisition unit 301, fitting unit 302, estimation unit 303, calculation unit 304, and inversion unit 305 in the methane inversion device 300 can be found in the following references: Figure 1 The relevant descriptions of steps 101, 102, 103, 104, and 105 in the corresponding embodiments will not be repeated here.

[0075] In one embodiment of this disclosure, the fitting unit 302 is configured to: use multi-temporal observation data to estimate the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient of the first band and the second band respectively by least squares fitting; based on the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient, use the arithmetic mean of the scene angles of the target temporal phase and the reference temporal phase as a reference geometric benchmark for reduction to obtain the corrected reflectance of the reflection distribution function that eliminates the observation geometric difference.

[0076] In one embodiment of this disclosure, the estimation unit 303 is configured to: correct reflectance based on the reflectance distribution function and calculate the residual of the current iteration; calculate a robust scaling estimate of the residual using the median absolute deviation; calculate the Hubel weight of each sample within the local window by combining a set Hubel threshold and the robust scaling estimate, wherein for outlier samples whose residuals are greater than the set threshold, their Hubel weights are reduced; and perform weighted least squares update using the Hubel weights to obtain new interspectral scaling coefficients until the iteration convergence condition is met or the maximum number of iterations is reached.

[0077] In some embodiments of this disclosure, the methane inversion apparatus 300 further includes a weighting unit (not shown in the figure), which is configured to: during the weighted least squares update process, when performing weighted summation calculation on samples within a local window, obtain the weighted summation result of samples within the local window based on a pre-constructed integral map of the full image, so that the complexity of the weighted summation calculation is independent of the size of the local window.

[0078] In some embodiments of this disclosure, the methane inversion device 300 further includes an expansion unit (not shown in the figure), which is configured to: determine whether the proportion of effective pixels within a local window is lower than a preset threshold; if it is lower than the preset threshold, use the effective estimation region as a seed, adopt an iterative filling strategy based on morphological expansion to expand outwards successively, and initialize the inter-spectral scaling coefficient of the invalid region with the mean of the outer window during each expansion, until all pixels obtain effective inter-spectral scaling coefficients.

[0079] In one embodiment of this disclosure, the inversion unit 305 is configured to: construct a cost function comprising observation terms and prior terms, wherein the observation terms are constructed based on the difference between the differential absorption model and the physical forward simulation value and the standard deviation of the observation noise, and the physical forward simulation value is obtained by interpolation of a lookup table; the prior terms are constructed based on the difference between the methane column concentration enhancement to be inverted and the prior mean and the prior standard deviation; the cost function is solved using the Levenberg-Marquardt algorithm, and the methane column concentration enhancement that minimizes the cost function is obtained by switching between the steepest descent method and the Gauss-Newton method using an adaptive damping factor.

[0080] In one embodiment of this disclosure, the strategy for determining the prior standard deviation includes an empirical Bayesian strategy or a data-driven adaptive strategy: the empirical Bayesian strategy includes scanning candidate prior standard deviations on a background image without methane plumes and selecting the candidate value that minimizes the background noise as the prior standard deviation; the data-driven adaptive strategy includes dynamically adjusting the prior standard deviation according to the magnitude of the observation noise standard deviation, and reducing the prior standard deviation to enhance the prior regularization constraint when the observation noise standard deviation increases, thereby avoiding divergence in the inversion solution under low signal-to-noise ratio conditions.

[0081] According to embodiments of this disclosure, this disclosure also provides an electronic device, a readable storage medium, and a computer program product.

[0082] Figure 4 A schematic block diagram of an example electronic device 400 that can be used to implement embodiments of the present disclosure is shown. The electronic device is intended to represent various forms of digital computers, such as laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device may also represent various forms of mobile devices, such as personal digital processors, cellular phones, smartphones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their patterns are merely examples and are not intended to limit the implementation of the present disclosure described and / or claimed herein.

[0083] like Figure 4 As shown, the electronic device 400 includes a computing unit 401, which can perform various appropriate actions and processes according to a computer program stored in a read-only memory (ROM) 402 or a computer program loaded from a storage unit 408 into a random access memory (RAM) 403. The RAM 403 may also store various programs and data required for the operation of the electronic device 400. The computing unit 401, ROM 402, and RAM 403 are interconnected via a bus 404. An input / output (I / O) interface 405 is also connected to the bus 404.

[0084] Multiple components in electronic device 400 are connected to I / O interface 405, including: input unit 406, such as keyboard, mouse, etc.; output unit 407, such as various types of displays, speakers, etc.; storage unit 408, such as disk, optical disk, etc.; and communication unit 409, such as network card, modem, wireless transceiver, etc. Communication unit 409 allows electronic device 400 to exchange information / data with other devices through computer networks such as the Internet and / or various telecommunications networks.

[0085] The computing unit 401 can be a variety of general-purpose and / or special-purpose processing components with processing and computing capabilities. Some examples of the computing unit 401 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various special-purpose artificial intelligence (AI) computing chips, various computing units running machine learning model algorithms, a digital signal processor (DSP), and any suitable processor, controller, microcontroller, etc. The computing unit 401 performs the various methods and processes described above, such as the methane inversion method. For example, in some embodiments, the methane inversion method may be implemented as a computer software program tangibly contained in a machine-readable medium, such as storage unit 408. In some embodiments, part or all of the computer program may be loaded and / or installed on the electronic device 400 via ROM 402 and / or communication unit 409. When the computer program is loaded into RAM 403 and executed by the computing unit 401, one or more steps of the methane inversion method described above may be performed. Alternatively, in other embodiments, the computing unit 401 may be configured as the methane inversion method by any other suitable means (e.g., by means of firmware).

[0086] Various embodiments of the systems and techniques described above herein can be implemented in digital electronic circuit systems, integrated circuit systems, field-programmable gate arrays (FPGAs), application-specific integrated circuits (ASICs), application-specific standard products (ASSPs), systems-on-a-chip (SoCs), payload-programmable logic devices (CPLDs), computer hardware, firmware, software, and / or combinations thereof. These various embodiments may include implementations in one or more computer programs that can be executed and / or interpreted on a programmable system including at least one programmable processor, which may be a dedicated or general-purpose programmable processor, capable of receiving data and instructions from a storage system, at least one input device, and at least one output device, and transmitting data and instructions to the storage system, the at least one input device, and the at least one output device.

[0087] The program code used to implement the methods of this disclosure may be written in any combination of one or more programming languages. This program code may be provided to the processor or controller of a general-purpose computer, special-purpose computer, or other programmable methane inversion apparatus, such that when executed by the processor or controller, the program code causes the patterns / operations specified in the flowcharts and / or block diagrams to be implemented. The program code may be executed entirely on a machine, partially on a machine, as a standalone software package partially on a machine and partially on a remote machine, or entirely on a remote machine or server.

[0088] In the context of this disclosure, a machine-readable medium can be a tangible medium that may contain or store a program for use by or in conjunction with an instruction execution system, apparatus, or device. A machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. A machine-readable medium can be, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. More specific examples of machine-readable storage media include electrical connections based on one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.

[0089] To provide interaction with a user, the systems and techniques described herein can be implemented on a computer having: a display device for displaying information to the user (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor); and a keyboard and pointing device (e.g., a mouse or trackball) through which the user provides input to the computer. Other types of devices can also be used to provide interaction with the user; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form (including sound input, voice input, or tactile input).

[0090] The systems and technologies described herein can be implemented in computing systems that include backend components (e.g., as a data server), or computing systems that include middleware components (e.g., an application server), or computing systems that include frontend components (e.g., a user computer with a graphical user interface or web browser through which a user can interact with implementations of the systems and technologies described herein), or any combination of such backend, middleware, or frontend components. The components of the system can be interconnected via digital data communication of any form or medium (e.g., a communication network). Examples of communication networks include local area networks (LANs), wide area networks (WANs), and the Internet.

[0091] It should be understood that the various forms of processes shown above can be used to rearrange, add, or delete steps. For example, the steps described in this disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this disclosure can be achieved, and this is not limited herein.

[0092] The foregoing description of specific exemplary embodiments of the invention is for illustrative and explanatory purposes. These descriptions are not intended to limit the invention to the precise forms disclosed, and it will be apparent that many changes and variations can be made in accordance with the foregoing teachings. The exemplary embodiments were chosen and described in order to explain the specific principles of the invention and its practical application, thereby enabling those skilled in the art to implement and utilize various different exemplary embodiments of the invention, as well as various different choices and variations. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A methane inversion method, characterized in that, The method includes: Acquire multispectral satellite observation data of the target time phase and the reference time phase, wherein the observation data includes at least the apparent reflectance of the first band and the apparent reflectance of the second band, which are sensitive to methane absorption. Using multi-temporal observation data, the reflectance model coefficients of the first band and the second band are independently fitted respectively. Then, using the reflectance model coefficients, the multispectral satellite observation data is uniformly corrected to the set reference geometric benchmark to obtain the reflectance corrected by the reflectance distribution function that eliminates observation geometric differences. The reflectance is corrected for the reflection distribution function. Within a local window centered on each pixel, the inter-spectral scaling coefficient between the first band and the second band is estimated based on a robust loss function and iterative reweighted least squares method, resulting in a pixel-by-pixel spatially varying inter-spectral scaling coefficient field. Based on the inter-spectral scaling coefficient field, the target fractional absorption model for the target time phase and the reference fractional absorption model for the reference time phase are calculated respectively; the difference between the target fractional absorption model and the reference fractional absorption model is obtained to obtain the differential absorption model. Within the maximum a posteriori probability framework, a nonlinear optimization inversion is performed on the differential absorption model based on a pre-constructed lookup table to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

2. The method according to claim 1, characterized in that, The process involves independently fitting reflectance model coefficients for the first and second bands using multi-temporal observation data, and then using these reflectance model coefficients to uniformly correct the multispectral satellite observation data to a set reference geometric benchmark, thereby obtaining a reflectance correction function that eliminates observational geometric differences. This includes: Using multi-temporal observation data, the isotropic scattering coefficient, volume scattering coefficient, and geometric scattering coefficient of the first and second bands are estimated pixel by pixel through least squares fitting. Based on the isotropic scattering coefficient, the volume scattering coefficient, and the geometric scattering coefficient, the arithmetic mean of the scene angles of the target time phase and the reference time phase is used as the reference geometric benchmark for reduction to obtain the reflection distribution function corrected reflectivity that eliminates observation geometric differences.

3. The method according to claim 1, characterized in that, The step of correcting the reflectance of the reflectance distribution function, centered on each pixel within a local window, and estimating the inter-spectral scaling factor between the first band and the second band based on a robust loss function and iterative reweighted least squares method, includes: The reflectivity is corrected based on the reflection distribution function, and the residual of the current iteration is calculated. A robust scaling estimate of the residuals is calculated using the median absolute deviation. By combining the set Hubel threshold and the robust scaling estimator, the Hubel weight of each sample within the local window is calculated. For outlier samples with residuals greater than the set threshold, their Hubel weights are reduced. The weighted least squares update is performed using the Hubel weights to obtain new interspectral scaling coefficients until the iterative convergence condition is met or the maximum number of iterations is reached.

4. The method according to claim 3, characterized in that, The method further includes: During the weighted least squares update process, when performing weighted summation calculation on samples within the local window, the weighted summation result of samples within the local window is obtained based on the pre-constructed integral map of the entire image, so that the complexity of the weighted summation calculation is independent of the size of the local window.

5. The method according to claim 1, characterized in that, The method further includes: Determine whether the proportion of effective pixels within the local window is lower than a preset threshold; If the value is below the preset threshold, the effective estimated region is used as the seed, and an iterative filling strategy based on morphological dilation is adopted to expand outwards one by one. During each expansion, the mean of the outer window is used to initialize the inter-spectral scaling factor of the invalid region until all pixels obtain an effective inter-spectral scaling factor.

6. The method according to claim 1, characterized in that, Within the maximum a posteriori probability framework, based on a pre-built lookup table, a nonlinear optimization inversion is performed on the differential absorption model to obtain the methane column concentration enhancement, including: A cost function is constructed that includes observation terms and prior terms. The observation terms are constructed based on the difference between the differential absorption model and the physical forward simulation values, as well as the standard deviation of the observation noise. The physical forward simulation values ​​are obtained by interpolating the lookup table. The prior terms are constructed based on the difference between the methane column concentration enhancement to be inverted and the prior mean, as well as the prior standard deviation. The cost function is solved using the Levenberg-Marquardt algorithm. By switching between the steepest descent method and the Gauss-Newton method using an adaptive damping factor, the methane column concentration enhancement that minimizes the cost function is obtained.

7. The method according to claim 6, characterized in that, The strategy for determining the prior standard deviation includes either an empirical Bayesian strategy or a data-driven adaptive strategy. The empirical Bayesian strategy includes: scanning candidate prior standard deviations on a background image without methane plumes, and selecting the candidate value that minimizes the background noise as the prior standard deviation. The data-driven adaptive strategy includes: dynamically adjusting the prior standard deviation according to the magnitude of the observation noise standard deviation; when the observation noise standard deviation increases, decreasing the prior standard deviation to enhance the prior regularization constraint, thereby avoiding divergence of the inversion solution under low signal-to-noise ratio conditions.

8. A methane inversion apparatus, characterized in that, The device includes: The acquisition unit is configured to acquire multispectral satellite observation data of the target time phase and the reference time phase, wherein the observation data includes at least the apparent reflectance of a first band and the apparent reflectance of a second band that are sensitive to methane absorption. The fitting unit is configured to independently fit the reflectance model coefficients of the first band and the second band using multi-temporal observation data, and use the reflectance model coefficients to uniformly correct the multispectral satellite observation data to a set reference geometric benchmark, thereby obtaining a reflectance correction reflectance that eliminates observation geometric differences; The estimation unit is configured to correct the reflectance for the reflection distribution function and estimate the inter-spectral scaling coefficient between the first band and the second band within a local window centered on each pixel, based on a robust loss function and iterative reweighted least squares method, to obtain a pixel-by-pixel spatially varying inter-spectral scaling coefficient field. The calculation unit is configured to calculate the target fractional absorption model of the target time phase and the reference fractional absorption model of the reference time phase based on the interspectral scaling coefficient field; and to obtain the differential absorption model by subtracting the target fractional absorption model and the reference fractional absorption model. The inversion unit is configured to perform a nonlinear optimization inversion of the differential absorption model based on a pre-built lookup table within a maximum a posteriori probability framework to obtain the methane column concentration enhancement. The lookup table is used to characterize the physical mapping relationship between remote sensing observations and methane concentration.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 7.

Images