Cloud computing method for regional daily evapotranspiration based on multi-source remote sensing data fusion

Abstract

The application discloses a regional daily evapotranspiration cloud computing method based on multi-source remote sensing data fusion. The method comprises the following steps: obtaining multi-source remote sensing images, a digital elevation model and site meteorological parameters of a target region; performing standardization preprocessing, pixel-level statistical synthesis and cutting on a remote sensing cloud computing platform to generate a target region fusion image; performing multi-channel feature extraction based on the fusion image to obtain pixel-by-pixel ground surface feature parameters; performing topographic relief dynamics three-dimensional spatialization correction on the site meteorological parameters by using the digital elevation model to obtain pixel-scale meteorological driving parameters; inputting each parameter into a ground surface energy balance system model, and solving instantaneous latent heat flux and evapotranspiration fraction based on dry and wet boundary constraints; and completing time scale expansion based on the relative stability assumption during the day and outputting daily evapotranspiration results. The application can solve the problems of data source revisit limitation and meteorological distortion in complex terrain areas, and realize efficient, high-precision and automatic inversion of regional evapotranspiration.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing inversion, surface energy balance calculation and geographic information processing, and in particular to a regional diurnal scale evapotranspiration cloud computing method based on multi-source remote sensing data fusion. Background Technology

[0002] Evapotranspiration, as a crucial variable connecting the surface water cycle, energy balance, and ecological processes, is a key input indicator for assessing water resource evolution trends, optimizing agricultural irrigation scheduling, and conducting in-depth analysis of eco-hydrological effects. Modern technology utilizes remote sensing to conduct macroscopic monitoring of regional evapotranspiration, overcoming the spatial discontinuity of traditional ground-based monitoring stations and obtaining large-scale, spatiotemporally continuous information on surface water consumption.

[0003] Current regional evapotranspiration inversion methods largely rely on a single satellite remote sensing data source to drive the surface energy balance model. In practical operational applications, this type of method has limitations. For example, a single sensor is limited by a fixed revisit period, and under variable weather conditions with frequent cloud cover or shadows, it is prone to generating large areas of missing effective observation data, resulting in discontinuous inversion time series. Summary of the Invention

[0004] The purpose of this invention is to provide a regional diurnal evapotranspiration cloud computing method based on multi-source remote sensing data fusion, in order to solve the aforementioned problems in the prior art.

[0005] Technical solution: A regional diurnal evapotranspiration cloud computing method based on multi-source remote sensing data fusion, comprising:

[0006] Acquire multi-source remote sensing images, digital elevation models, diurnal radiation data, and station meteorological parameters within the target area;

[0007] Multi-source remote sensing images are preprocessed using a remote sensing cloud computing platform to obtain a standardized image set;

[0008] A pre-defined statistical synthesis strategy is used to perform pixel-level statistical synthesis on a standardized image set and then crop it to obtain a fused image of the target area. Based on this, multi-channel feature extraction is performed to obtain pixel-by-pixel surface feature parameters.

[0009] Based on the digital elevation model, the meteorological parameters of the stations are subjected to three-dimensional spatial correction of topographic undulation dynamics to obtain the meteorological driving parameters at the pixel scale.

[0010] The pixel-by-pixel surface feature parameters, digital elevation model, and pixel-scale meteorological driving parameters are input into the surface energy balance system model to solve for the instantaneous latent heat flux and evapotranspiration fraction.

[0011] Based on the assumption that the evapotranspiration fraction is relatively stable during the day, the instantaneous latent heat flux is combined with the available energy at the day-scale calculated based on day-scale radiation data to extend the time scale and output the regional day-scale evapotranspiration results.

[0012] Beneficial effects include solving the problems of data source revisit limitations and meteorological distortion in complex terrain areas, and achieving high-efficiency and high-precision automated inversion of regional evapotranspiration. Attached Figure Description

[0013] Figure 1 This is a flowchart of the regional daily-scale evapotranspiration cloud computing method based on multi-source remote sensing data fusion according to the present invention.

[0014] Figure 2 This is a flowchart of the process of obtaining a standardized image set according to the present invention.

[0015] Figure 3 This is a flowchart of the present invention for obtaining pixel-by-pixel surface feature parameters.

[0016] Figure 4 This is a flowchart of the present invention for obtaining pixel-scale meteorological driving parameters.

[0017] Figure 5 This is a flowchart of the present invention for solving the instantaneous latent heat flux and evapotranspiration fraction of the target region. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0019] As mentioned in the background section, current regional evapotranspiration inversion relies on a single satellite remote sensing data source to drive the surface energy balance model, which has technical limitations. Through analysis and research, the applicant also discovered other shortcomings. For example, the processing flow for regional evapotranspiration inversion involving different surface characteristic parameters and station meteorological driving data is independent, resulting in fragmented local computations and difficulty in coordinated data scheduling. Furthermore, this approach often focuses on instantaneous flux calculations at the moment of transit, lacking the ability to downscale and reconstruct meteorological elements in accordance with the terrain in areas with complex topographic relief, making it difficult to stably obtain high-frequency, physically consistent diurnal inversion data.

[0020] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0021] Example 1

[0022] This embodiment provides a regional diurnal evapotranspiration cloud computing method based on multi-source remote sensing data fusion, which includes the following steps:

[0023] Step 101: Acquire multi-source remote sensing images, digital elevation models, daily-scale radiation data, and station meteorological parameters within the target area.

[0024] In this embodiment, the target area refers to the preset geographic spatial range for which evapotranspiration inversion is to be performed, and a continuous observation period between preset start and end dates needs to be set as the target time period. The acquired multi-source remote sensing images have complementarity in spatial and temporal observation frequencies and can cover Landsat image datasets and Sentinel-2 image datasets. The digital elevation model is used to provide spatially distributed three-dimensional geographic information of surface undulations. Station meteorological parameters are provided by real-time collection or pre-storage from meteorological observation stations distributed within the region, including air temperature, relative humidity, air pressure, wind speed, and downwave radiation basic meteorological elements measured at discrete time points.

[0025] Correspondingly, daily-scale radiation parameter data for the target time period can also be obtained. This data can be obtained by time integration based on hourly observations using a surface radiation balance model, or it can be extracted from publicly available reanalysis meteorological datasets.

[0026] Step 102: Preprocess the multi-source remote sensing images based on the remote sensing cloud computing platform to obtain a standardized image set.

[0027] The remote sensing cloud computing platform possesses massive distributed computing and storage capabilities and can be deployed using mainstream distributed cloud architectures. During preprocessing tasks, it leverages the platform's server-side lazy evaluation mechanism to construct a full-stack automated computation graph. This mechanism allocates computing power in the cloud to perform preprocessing operations such as cleaning, masking for cloud removal, radiometric calibration, and reflectivity scaling, avoiding the download of massive amounts of raster data generated during intermediate processes to local terminals. This embodiment outputs a standardized image set with a unified matrix structure by implementing band alignment and null-filling strategies across multiple satellite platforms, for subsequent heterogeneous data fusion.

[0028] Step 103: Perform pixel-level statistical synthesis on the standardized image set using a preset statistical synthesis strategy and then perform cropping to obtain the fused image of the target area.

[0029] Due to limitations in satellite transit cycles and cloud cover interference during the target time period, remote sensing images from a single transit often suffer from significant missing effective pixels. Performing temporal-dimensional pixel-level statistical synthesis on a standardized image set in the cloud can remove residual cloud-contaminated pixels and outliers, reconstructing a spatially continuous and temporally representative surface reflectance base map. Correspondingly, based on the vector boundary file of the target area, spatial boundary clipping is performed on the synthesized base map to output a spatially aligned fused image of the target area.

[0030] Step 104: Perform multi-channel feature extraction based on the fused image of the target area to obtain pixel-by-pixel surface feature parameters.

[0031] Solving the surface energy balance system relies on multi-dimensional surface physical parameters. In this embodiment, the optical reflectance and thermal infrared radiation information of different spectral channels in the fused image of the target area are analyzed to recursively solve for the normalized vegetation index, surface albedo, vegetation cover, surface emissivity, and surface temperature. These parameters are extracted and constructed into a two-dimensional raster matrix set, which is then encapsulated as pixel-by-pixel surface feature parameters and input into the downstream computing engine.

[0032] Step 105: Based on the digital elevation model, perform three-dimensional spatialization correction of the meteorological parameters of the station on topographic relief dynamics to obtain the pixel-scale meteorological driving parameters.

[0033] Traditional evapotranspiration inversion typically applies data from a single weather station to the entire image area, which can easily lead to thermodynamic and aerodynamic distortions in environments with complex terrain.

[0034] Therefore, in this embodiment, the elevation difference information of each pixel in the digital elevation model is extracted, and the meteorological elements are spatially downscaled and reconstructed in combination with atmospheric laws. Elements such as temperature, air pressure, and air density are transformed into a three-dimensional spatial distribution field strongly correlated with terrain, generating pixel-scale meteorological driving parameters, which can provide high-precision physical boundary conditions for the energy model.

[0035] Step 106: Input the pixel-by-pixel surface feature parameters, digital elevation model and pixel-scale meteorological driving parameters into the surface energy balance system model, and solve the instantaneous latent heat flux and evapotranspiration fraction of the target area based on the dry and wet boundary constraints.

[0036] Specifically, the surface energy balance system model uses pixel-scale meteorological driving parameters and pixel-by-pixel surface characteristic parameters as joint inputs to calculate components such as net surface radiation and soil heat flux. When extrapolating sensible heat flux, limiting sensible heat flux reference frames are constructed under dry and wet boundary conditions, respectively, and the actual sensible heat flux is dynamically constrained within the constraint intervals constructed by these two limiting reference frames. This embodiment solves for relative evaporation through a dual boundary constraint mechanism, and, combined with the energy balance relationship, outputs the instantaneous latent heat flux and corresponding evapotranspiration fraction at the moment of satellite transit.

[0037] Furthermore, it also includes the introduction of surface roughness parameters, specifically including:

[0038] Step 1061: Extract the normalized vegetation index from the pixel-by-pixel surface feature parameters.

[0039] Surface roughness characterizes the physical shearing and impeding effect of the surface on the underlying airflow and is used in boundary layer dynamics model calculations. In this embodiment, a set of pixel-by-pixel surface feature parameter matrices is extracted, and a normalized vegetation index matrix reflecting the morphological density of the vegetation canopy is extracted.

[0040] Step 1062: Calculate the leaf area index and momentum roughness length of the target area based on the normalized vegetation index.

[0041] In this embodiment, the leaf area index is used to characterize the cross-sectional area of ​​contact exchange between vegetation and the atmosphere, and the momentum roughness length is related to the mechanical drag characteristics of the ground surface. The normalized vegetation index is used to drive a preset empirical logarithmic model to calculate the momentum roughness length at the pixel scale, and the leaf area index distribution matrix of each pixel is derived based on the transformation parameters.

[0042] Step 1063: Calculate the thermal roughness length based on the momentum roughness length and the preset thermal impedance parameters.

[0043] In this embodiment, the thermal roughness length is used to describe the additional interfacial impedance of laminar heat transfer. Specifically, preset heat transfer impedance parameters are read, and the momentum roughness length is converted into a thermal roughness length matrix through an exponential operation mechanism to complete the surface boundary parameter set of the aerodynamic impedance model.

[0044] Step 1064: The surface roughness parameters, including leaf area index, momentum roughness length, and thermal roughness length, along with the pixel-scale meteorological driving parameters, are synchronously input into the surface energy balance system model.

[0045] In this embodiment, the generated pixel-scale air temperature, air pressure, air density, and water vapor pressure deficits can be used to eliminate the spurious thermodynamic homogeneity caused by complex terrain. The generated leaf area index and roughness parameter matrix describe the microscopic surface drag properties. This physical matrix is ​​used as a joint driving condition and input into the surface energy balance system model for subsequent spatial iterative solutions of radiation balance and surface heat flux components.

[0046] In one alternative implementation, the heat transfer impedance parameter is used to characterize the difference in momentum transfer and heat transfer efficiency, and its specific parameter value can be determined based on the surface roughness element characteristics through existing empirical or semi-empirical models.

[0047] Step 107: Based on the assumption that the evapotranspiration fraction is relatively stable during the day, the instantaneous latent heat flux is combined with the available energy at the day-scale calculated based on the day-scale radiation data to extend the time scale and output the regional day-scale evapotranspiration results.

[0048] Since the instantaneous latent heat release rate is insufficient to support the needs of water resource assessment and agricultural irrigation scheduling, this embodiment assumes that the evapotranspiration fraction in nature is relatively stable during normal daytime periods. Based on this physical characteristic, the instantaneous evapotranspiration ratio derived from the satellite's transit time is extrapolated to the entire day. This stable ratio coefficient is then integrated with the entire day's data to obtain the available energy on a daily scale, which is then converted from instantaneous physical flux to cumulative daily consumption. The output is packaged into multi-channel raster data and connected to a cloud application terminal or geographic information system platform.

[0049] Example 2

[0050] Based on the above embodiment one, this embodiment further explains the process of standardization preprocessing and fusion cropping of multi-source remote sensing images.

[0051] In one possible implementation, the multi-source remote sensing imagery includes Landsat imagery and Sentinel-2 imagery; the multi-source remote sensing imagery is preprocessed using a remote sensing cloud computing platform to obtain a standardized image set, including the following steps:

[0052] Step 201: Perform bitwise operations on the pixel quality assessment bands of the Landsat image, identify cloud and shadow pixels and perform mask removal, perform radiometric calibration on the optical bands based on preset scaling factors and offsets, and perform brightness temperature conversion on the thermal infrared bands to obtain calibrated land images.

[0053] Specifically, the system loads a Landsat surface reflectance image dataset from a cloud drive or distributed image database and extracts a pixel quality assessment band. This band records the quality characteristic information of each pixel in the form of a binary bitmask. Correspondingly, it reads the preset cloud data bit index and shadow data bit index and performs a bitwise AND operation between the pixel quality assessment band and the index. When the operation result indicates that the current pixel is covered by clouds or is in terrain shadow, it is determined to be an invalid pixel, and all spectral channel values ​​of that pixel are assigned null values ​​or masked for removal.

[0054] To reduce the effects of atmospheric absorption and scattering, in this embodiment, a preset scaling factor and offset corresponding to each optical band are obtained, and radiometric calibration is performed on the Landsat optical bands selected by the mask to convert their original digital quantization values ​​into real surface reflectance values.

[0055] Furthermore, the radiance image of the thermal infrared spectral band is loaded, and nonlinear calculations are performed using preset radiative correction constants and brightness temperature conversion models to convert the radiance of the thermal infrared channel into the corresponding thermal infrared brightness temperature field, outputting a calibrated land image that includes the surface reflectance spectral band and the brightness temperature band.

[0056] Step 202: Identify cloud, shadow, and cirrus cloud pixels based on the scene classification band of Sentinel-2 image and perform mask removal. Scale the reflectivity of the optical band based on a preset scaling factor to obtain a cloud-free Sentinel image.

[0057] To simultaneously clean Sentinel-2 imagery with higher spatial resolution and revisit frequency, this embodiment retrieves a set of Sentinel-2 imagery within a target time period and target area, extracting scene classification bands for each image. Scene classification bands refer to a pixel attribute identifier matrix generated by a classification algorithm. Pixels identified as clouds, terrain shadows, and upper-level cirrus clouds are identified, and their values ​​are assigned null values ​​in all spectral channels using a matrix mask operator. For the noise-removed blue, green, red, near-infrared, and multi-channel short-wave infrared optical bands, a preset scaling factor is obtained, and the pixel values ​​are scaled proportionally to ensure their dimensions and numerical ranges match the optical reflectance of the calibrated land imagery.

[0058] Unlike existing technologies that use only a single satellite data source, this embodiment proposes the following brightness temperature occupancy band construction strategy to address the issue of inconsistent band structures between different satellite sensors, specifically including:

[0059] Step 203: Construct a brightness temperature occupancy band with the same band structure as the cloud-free sentinel image of the non-thermal infrared band, and uniformly assign null values ​​or mask values ​​to the pixel values ​​of the brightness temperature occupancy band to obtain the calibrated sentinel image.

[0060] Since the Sentinel-2 satellite does not carry a thermal infrared sensor, its raw data structure lacks the brightness temperature spectrum band required for retrieving surface temperature. To achieve cross-platform alignment of multi-source satellites in the cloud data matrix, this embodiment employs a structure alignment mechanism during the preprocessing stage.

[0061] Specifically, for the processed Sentinel-2 imagery, new matrix channels are dynamically created as brightness temperature occupant bands. The number of rows and columns of the brightness temperature occupant bands is consistent with the current Sentinel imagery, serving as placeholders for the band structure. To reduce interference from occupant channels lacking thermal infrared information in subsequent surface temperature calculations, the pixel values ​​of each geographic coordinate point within the brightness temperature occupant band are uniformly assigned specific null value identifiers, specifically non-numeric symbols or preset no-data codes. After completing the topology stitching of this channel, a calibrated Sentinel imagery containing optical bands and virtual brightness temperature occupant channels is output, ensuring that the multi-source remote sensing images have the same matrix structure in the channel dimension.

[0062] Step 204: Unify the band naming of the calibration land image and the calibration sentinel image, resample the spatial resolution of the calibration sentinel image to be consistent with that of the calibration land image, and obtain the standardized image set after performing pixel grid alignment.

[0063] To reduce the differences in channel naming conventions between different satellite sensors, corresponding spectral bands in the calibrated land imagery and the calibrated sentinel imagery are renamed and uniformly mapped to standard abbreviations, including blue light band, green light band, red light band, near-infrared band, first shortwave infrared band, second shortwave infrared band, and brightness temperature band. Furthermore, since the spatial resolution of the sentinel-2 imagery is higher than that of the land satellite imagery, a preset resampling operator is invoked to uniformly reduce the spatial resolution of each band in the calibrated sentinel imagery and resample it to a standard resolution consistent with that of the calibrated land imagery in order to achieve pixel-level temporal synthesis. In this embodiment, the standard resolution can be adjusted according to the spatial resolution requirements of the actual application.

[0064] Furthermore, using the geographic projection coordinate system and pixel tile grid of the calibrated land image as a reference, spatial position correction and projection transformation are performed on the resampled calibrated sentinel image to achieve grid alignment of the two data sources in geographic space and pixel grid, and then merge and package them to generate a regular standardized image set.

[0065] In some optional implementations, the data acquisition end can also access other multispectral satellite data sources with visible, near-infrared, and short-wave infrared bands, such as Gaofen series satellite imagery. When using this data source to replace remote sensing imagery, the corresponding radiometric calibration factor, offset, and band mapping rules need to be adjusted according to the sensor characteristics of the newly accessed satellite. For replacement satellites that do not carry thermal infrared bands, the brightness temperature occupancy band construction process still needs to be performed, and the remaining processes such as spatial alignment and data cropping remain consistent. That is, the cloud computing framework proposed in this embodiment has good compatibility with multi-source heterogeneous data and functional scalability.

[0066] Example 3

[0067] Building upon Embodiment 2 above, this embodiment further provides an adaptive pixel-level statistical synthesis scheme to address different cloud coverage conditions. In one possible implementation, a preset statistical synthesis strategy is used to perform pixel-level statistical synthesis on a standardized image set and then cropping it to obtain the pixel-level statistical synthesis in the fused image of the target region, including the following steps:

[0068] Step 301: Merge the Landsat image set and the Sentinel-2 image set in the standardized image set.

[0069] Specifically, a set of Landsat images and a set of Sentinel-2 images, both processed with radiometric calibration, reflectivity scaling, and band structure alignment, are acquired. Along the time dimension, the two independent image sets are sequentially stitched together to construct a time-series image data cube containing multi-source observation frequencies. This increases the effective observation frequency of the target area within a preset target time period, which is used for subsequent removal of transient meteorological disturbances and extraction of stable surface features.

[0070] Step 302: Using any one of the median composite method, maximum vegetation index composite method, weighted average method, or cloud probability weighted method, perform temporal pixel-level statistical composite on the merged image set, automatically remove null values ​​or mask values, and generate time-series fused images.

[0071] Specifically, the proportion of cloud-contaminated pixels in the merged image set within the target time period is statistically analyzed. This involves reading the mask label matrix generated during the preprocessing stage from the merged time-series image data cube. The total number of pixels marked as null values ​​or cloud shadow mask values ​​within the target area is obtained, and the ratio of this value to the total number of pixels in the target area is calculated to obtain the average cloud contamination coverage proportion within the target time period, thus assessing the usability of the remote sensing data for the current time period.

[0072] Accordingly, when the coverage ratio is lower than a preset cloud cover threshold, median composite method is used for statistical synthesis. That is, the preset cloud cover threshold is an empirical value; when the statistically obtained coverage ratio is lower than this threshold, it indicates that the meteorological conditions within the current target time period are relatively good, effectively representing relatively sufficient observation data. Along the time dimension, all effective pixel values ​​at the same geographic coordinate point are sorted, and the reflectance, intensity, or brightness temperature values ​​located in the middle position are extracted and used as the final composite value for that coordinate point. In this embodiment, the median composite method can smooth out local abnormal fluctuations caused by atmospheric fluctuations or minor sensor errors, outputting a time-series fused image with high spatiotemporal continuity and stability. Simultaneously, it automatically filters out residual sporadic null values ​​or mask markers using the sorting distribution pattern.

[0073] When the coverage ratio is higher than or equal to the cloud cover threshold, the system automatically switches to either the maximum vegetation index synthesis method or the cloud probability weighted method for statistical synthesis to ensure the continuity of the spatiotemporal sequence of land surface parameters. In other words, when the coverage ratio reaches or exceeds the preset cloud cover threshold, it indicates that there is persistent cloud cover during the target time period, resulting in insufficient effective pixels. Continuing to use the median synthesis method could easily lead to the synthesis results deviating from the true land surface condition. The judgment condition is met, and the synthesis operator is switched.

[0074] In an alternative approach, the maximum vegetation index synthesis method is employed. This involves calculating the normalized vegetation index (NDI) in the pixel time series, extracting the observation time at which the NDI value reaches its maximum, and using all spectral band values ​​corresponding to that time as the synthesized value for that coordinate point. In other words, when the measured area is cloud-polluted, the reduced spectral reflectance of the vegetation index is utilized to capture a clearer surface observation slice by locking onto the maximum value.

[0075] In an alternative approach, a cloud probability weighting method is used to perform statistical synthesis. Specifically, based on a pre-configured cloud detection probability layer, a weight coefficient that is non-linearly negatively correlated with the cloud contamination probability is assigned to each pixel in each time slice. The spectral values ​​of the pixels at each observation time are multiplied by their corresponding weight coefficients and then summed to generate a comprehensively weighted time-series fused image. This improves the utilization rate of limited clear-sky observation data under adverse weather conditions and maintains the spatiotemporal continuity of the input boundaries for subsequent inversion processes.

[0076] Furthermore, in one possible implementation, a pre-defined statistical synthesis strategy is used to perform pixel-level statistical synthesis on a standardized image set and then cropping is performed to obtain the cropped image of the target region. Specifically, this includes:

[0077] Step 303: Load the pre-stored target region vector boundary file; use the target region vector boundary file to perform boundary cropping on the time-segment fused image to obtain a target region fused image with regular spatial boundaries.

[0078] Specifically, the pre-stored target area vector boundary file is loaded, and the spatial geometric topology information of the polygons in the vector boundaries is extracted to generate a binary spatial mask matrix with the same resolution as the time-series fused image. Using the spatial mask matrix, matrix bitwise multiplication is performed on the time-series fused image, retaining the valid pixel values ​​inside the polygon boundaries and clearing or assigning null values ​​to the pixel values ​​outside the boundaries, thereby obtaining a well-formed target area fused image on the geographic spatial range and data boundaries, which is used to drive the subsequent batch inversion of surface physical parameters and energy model.

[0079] Example 4

[0080] Based on Example 1, this example further explains the specific algorithm implementation for multi-channel feature extraction. In one possible implementation, multi-channel feature extraction is performed based on the fused image of the target area to obtain pixel-by-pixel surface feature parameters, including the following steps:

[0081] Step 401: Extract the near-infrared and red light bands from the fused image of the target area, and calculate the normalized vegetation index based on the ratio between the two.

[0082] Specifically, the spectral matrix of the fused image of the target area is read, and the near-infrared and red band reflectances are separated. To characterize the growth status and cover intensity of surface vegetation, the difference between the two band reflectances is calculated and used as the numerator. The sum of the two is calculated and used as the denominator. The ratio of the numerator to the denominator is then used to achieve spatial normalization. This operation reduces some of the radiation errors caused by topographic relief and changes in solar altitude angle, resulting in a pixel-by-pixel normalized vegetation index matrix.

[0083] Step 402: Extract blue light, red light, near-infrared and multi-channel shortwave infrared bands in parallel, combine the reflectance of each band and remove the system bias term, that is, extract the blue light band, red light band, near-infrared band and multi-channel shortwave infrared band from the fused image of the target area, use the preset empirical weighted model to linearly combine the reflectance of each band and remove the system bias term to obtain the surface albedo.

[0084] In this embodiment, surface albedo is used to describe the overall reflectivity of the Earth's surface to shortwave solar radiation. Blue, red, near-infrared, first shortwave infrared, and second shortwave infrared band data are extracted simultaneously from the fused image. Based on the surface cover type of the target area, a set of transformation weighting coefficients related to the spectral response function of each channel are pre-configured. The reflectivity matrix of each band is multiplied by the corresponding preset weighting coefficient, and matrix accumulation is performed. Correspondingly, a constant-form systematic bias term is subtracted to correct for the residual effects of atmospheric path radiation, thereby outputting a surface albedo field representing broadband reflectivity characteristics.

[0085] Step 403: Spatial normalization of the normalized vegetation index is performed, that is, using the preset soil end-member vegetation index and vegetation end-member vegetation index, the normalized vegetation index is spatially normalized and the 0 to 1 interval is truncated to obtain the vegetation coverage.

[0086] In this embodiment, vegetation cover reflects the proportion of the vertically projected area of ​​the vegetation canopy within a pixel. Specifically, soil end-member vegetation index and vegetation end-member vegetation index are obtained based on the geographical environment of the study area. For each pixel, the difference is obtained by subtracting the soil end-member vegetation index from its normalized vegetation index, and then dividing this difference by the difference between the vegetation end-member and soil end-member vegetation indices. Considering that water bodies, abnormally bright features, or dense forests may cause the calculation results to overflow the theoretical boundary, it is determined whether the normalized values ​​are within the theoretical range. Values ​​below a preset lower limit are reset to the lower limit value, and values ​​above a preset upper limit are reset to the upper limit value, thus completing the truncation limit of the 0 to 1 interval and outputting a vegetation cover that conforms to the physical reality.

[0087] Step 404: Based on vegetation cover, perform weighted summation using a proportional coefficient. That is, use a preset background value of bare soil specific emissivity and a total vegetation cover proportional coefficient to perform weighted summation to obtain the surface specific emissivity.

[0088] In this embodiment, surface emissivity characterizes the efficiency of heat energy radiated outward from the target surface. Based on vegetation cover, it is considered as the area weight of the surface vegetation component and the bare soil component. The vegetation cover is multiplied by a preset total vegetation cover ratio coefficient to obtain the emissivity increment of the vegetation portion. This increment is added to a preset background value of bare soil emissivity to achieve linear weighted reconstruction of the radiation characteristics of different land features within the mixed pixel, thereby obtaining a surface emissivity matrix that is positively correlated with the density of surface vegetation growth.

[0089] Step 405: Extract the effective brightness temperature band from the fused image of the target area and perform nonlinear conversion on it, that is, combine the surface emissivity and the preset single-window radiation transfer simplification constant to perform nonlinear conversion to obtain the surface temperature.

[0090] In this embodiment, effective brightness temperature data observed by thermal infrared sensors is extracted from the fused image matrix. The brightness temperature detected by the satellite in orbit fails to account for differences in surface emissivity, easily leading to radiation errors. Therefore, the system invokes a single-window radiative transfer mechanism for temperature inversion. Specifically, the natural logarithm of the surface emissivity is calculated, multiplied by the effective brightness temperature band matrix, and a simplified constant derived from Planck's law. This product is incremented by 1 and used as the denominator, with the original effective brightness temperature matrix as the numerator for division. This reduces the thermal radiation deviation caused by uneven surface emissivity, restoring the true pixel-by-pixel surface temperature field. The corresponding calculation formula is:

[0091] T_s =BT / (1+(λ×BT / ρ _rad )×ln(ε _s ));

[0092] ρ _rad =h×c / σB;

[0093] Where BT represents the effective brightness temperature band data, λ is the center wavelength of the thermal infrared sensor, h is Planck's constant, c is the speed of light, σB is Boltzmann's constant, and ε _s ρ is the surface emissivity. _rad is Planck's second radiation constant.

[0094] Step 406: Integrate the five parameters, in other words, integrate the normalized vegetation index, surface albedo, vegetation cover, surface emissivity, and surface temperature as pixel-by-pixel surface feature parameters.

[0095] After completing the recursive inversion of all surface optical and thermodynamic parameters, the independently generated two-dimensional raster matrices are stacked and stitched together along the channel dimension. This process integrates the dispersed physical quantities into a multidimensional feature tensor, which serves as a unified data interface for downstream processing, used to solve for pixel-by-pixel parameters of surface energy balance.

[0096] In one optional implementation, the channel weighting coefficients and system bias terms of the empirical weighted model can be determined based on the spectral response function of the target sensor using conventional narrow-band to wide-band conversion empirical relationships. The soil end-member vegetation index and vegetation end-member vegetation index can be determined through statistical analysis based on the land cover type of the study area. The background value of bare soil emissivity and the full vegetation cover ratio coefficient can be determined based on the thermal infrared spectral characteristics of the target land surface type using conventional empirical models of surface emissivity. The simplified constant of single-window radiative transfer can be derived from the ratio of the second radiation constant of Planck's radiation law to the center wavelength of the thermal infrared band, and directly calculated based on the center wavelength of the thermal infrared band of the satellite used.

[0097] Example 5

[0098] Unlike traditional methods that assume spatial uniformity of meteorological elements, this embodiment physically reconstructs the three-dimensional spatial correction process of station meteorological parameters and the derivation of dynamic roughness parameters. In one possible implementation, the station meteorological parameters are subjected to three-dimensional spatial correction based on topographic relief dynamics using a digital elevation model to obtain pixel-scale meteorological driving parameters, including the following steps:

[0099] Step 501: Obtain the measured air temperature, measured air pressure, and measured relative humidity from the meteorological parameters of the station.

[0100] Atmospheric environmental attributes at discrete time points are extracted from data streams from cloud-based meteorological databases or automatic weather stations. Measured air temperature, air pressure, and relative humidity recorded by the measurement stations during that time period are extracted and used as initial reference values ​​for subsequent topographic dynamics downscaling and three-dimensional spatial extrapolation.

[0101] Furthermore, based on the digital elevation model, pixel-by-pixel spatial downscaling and dynamic correction are performed on the measured air temperature, measured air pressure, and measured relative humidity through surface elevation undulations. Specifically, spatial downscaling and correction are performed on air temperature and air pressure based on the digital elevation model, and pixel-scale water vapor pressure deficit is calculated based on the corrected relative humidity.

[0102] Step 502: Calculate the relative elevation difference between the elevation of the corresponding pixel in the digital elevation model and the altitude of the preset weather station, and combine the preset standard temperature vertical lapse rate to correct the measured air temperature pixel by pixel to obtain the air temperature at the pixel scale.

[0103] Specifically, the elevation values ​​of each independent pixel in the digital elevation model are extracted, and the altitude of the benchmark meteorological station providing the measured meteorological data is read. The absolute physical difference between the pixel elevation and the meteorological station altitude is calculated. In the lower troposphere, air temperature exhibits an approximately linear decay characteristic with increasing altitude. The temperature compensation term is obtained by multiplying the elevation difference by the standard temperature lapse rate. This temperature compensation term is subtracted from the measured air temperature to output the pixel-scale air temperature corresponding to the geographic coordinates. The corresponding calculation formula is as follows:

[0104] T _a (x,y)=T _a_stn -Г×(DEM(x,y)-z _stn );

[0105] Among them, T _a (x,y) represents the pixel-scale air temperature at pixel (x,y), T _a_stn Here, Γ represents the measured air temperature at the weather station, Г represents the standard temperature vertical lapse rate, DEM(x,y) represents the elevation at pixel (x,y), and z represents the vertical elevation at the station. _stn The altitude of the preset weather station is used. In this embodiment, the parameter value of the standard temperature vertical lapse rate can be adaptively adjusted according to the actual atmospheric conditions of the target area.

[0106] Step 503: Based on the relationship between air pressure and altitude, and combined with gravitational acceleration, dry air gas constant, and the ratio of air temperature before and after correction, the measured air pressure is corrected for altitude correlation to obtain the air pressure at the pixel scale.

[0107] In real-world environments, surface air pressure is not only affected by altitude, but also exhibits a nonlinear coupling relationship with the vertical temperature profile of the air column.

[0108] Calculate the ratio of the spatially corrected pixel-scale air temperature to the measured air temperature at the station. Using gravitational acceleration as the numerator and the product of the dry air gas constant and the vertical lapse rate of standard air temperature as the denominator, calculate the ratio and use it as the dynamic power exponent. Perform this power exponentiation on the air temperature ratio, then multiply it by the measured air pressure at the station to output the pixel-scale air pressure with three-dimensional spatial consistency. The corresponding calculation formula is as follows:

[0109] P(x,y)=P _stn ×(T _a (x,y) / T _a_stn ) (g / (Rd×Г)) ;

[0110] Where P(x,y) is the pixel-scale air pressure at pixel (x,y), P _stn For the actual air pressure measured at the weather station, T _a (x,y) represents the corrected pixel-scale air temperature, T _a_stn denoted as the measured air temperature at the weather station, g is the gravitational acceleration, and Rd is the gas constant for dry air.

[0111] Step 504: Substitute the pixel-scale air temperature and pixel-scale air pressure into the ideal gas law to deduce the physical spatial distribution and obtain the pixel-scale air density.

[0112] After obtaining the continuous temperature and pressure fields, they are input into the ideal gas law to deduce the local air density. Specifically, the dry air gas constant is multiplied by the pixel-scale air temperature as the equivalent energy denominator, and the pixel-scale air pressure is used as the numerator for division. The resulting spatial distribution field of air density covering the target area is output, and its calculation formula is as follows:

[0113] ρ _a (x,y)=P(x,y) / (Rd×T _a (x,y));

[0114] Where, ρ _a (x,y) represents the pixel-scale air density at pixel (x,y), and P(x,y) represents the pixel-scale air pressure.

[0115] Step 505: Calculate the saturated water vapor pressure of the corresponding pixel based on the air temperature at the pixel scale, and calculate the actual water vapor deficit by combining the measured relative humidity, thus obtaining the water vapor pressure deficit at the pixel scale.

[0116] The vapor pressure deficit is a key hydrometeorological impedance element driving evapotranspiration. In this embodiment, the system uses pixel-scale air temperature as the driving variable and employs the Clausius-Clapeyron empirical relationship to calculate the saturated vapor pressure of the pixel under the current temperature boundary. Specifically, assuming that the relative humidity in the target area is spatially approximately uniform, the measured relative humidity of the station is normalized by dividing it by a constant of 100. The relative dryness is calculated by subtracting this normalized value from the constant 1. This dryness is then multiplied by the saturated vapor pressure to output the vapor pressure deficit field describing the atmospheric water absorption potential. The calculation formula is as follows:

[0117] VPD(x,y)=e _s (x,y)×(1-RH _stn / 100);

[0118] Where VPD(x,y) is the pixel-scale vapor pressure deficit at pixel (x,y), e _s (x,y) represents the saturated vapor pressure (RH) calculated based on the pixel-scale air temperature. _stn The formula represents the relative humidity measured at the weather station. In this formula, 100 is the conversion factor for percentages to decimals.

[0119] In this embodiment, considering the limited spatial range of the target area, typically a single irrigation district or watershed, it is assumed that the spatial variation of relative humidity in the horizontal direction can be ignored, i.e., the measured relative humidity at the station is applied to all pixels, and the saturated vapor pressure is calculated using the Tetens empirical formula.

[0120] Step 506: Generate pixel-scale meteorological driving parameters, including pixel-scale air temperature, pixel-scale air pressure, pixel-scale air density, and pixel-scale water vapor pressure deficit.

[0121] Specifically, the dispersed atmospheric state variables are repackaged and encapsulated to output pixel-scale meteorological driving parameters, and the parameter matrix is ​​spatially aligned with the remote sensing image grid.

[0122] In this embodiment, before performing the sensible heat flux extrapolation at the wet-dry boundary, the system iteratively solves for the aerodynamic impedance based on measured wind speed, momentum roughness length, and thermal roughness length using the logarithmic wind profile equation and the Moning-Obukhov similarity theory. This iterative solution process is a standard calculation method in boundary layer meteorology in this field and can be routinely implemented based on pixel-scale meteorological driving parameters and roughness parameters.

[0123] Example 6

[0124] This embodiment further explains the internal solution logic of the Earth's surface energy balance system model.

[0125] In one possible implementation, before extrapolating the sensible heat flux of the target area under dry and wet boundary conditions, the surface energy balance system model further includes solving for the soil heat flux, specifically including the following steps:

[0126] Step 601: Based on the surface albedo, surface emissivity, surface temperature and pixel-scale meteorological driving parameters in the pixel-by-pixel surface characteristic parameters, the net radiation of the target area is calculated using the radiative transfer equation.

[0127] In this embodiment, net radiation is used to characterize the energy balance difference between the solar shortwave radiation absorbed by the Earth's surface and the longwave radiation emitted outward from the Earth's surface. Specifically, surface albedo, surface emissivity, and corrected pixel-scale air temperature and surface temperature are obtained, along with measured downwave radiation and measured wind speed from the station's meteorological parameters. Based on pixel-scale air temperature and pixel-scale vapor pressure deficit, the atmospheric emissivity of the target area is estimated using a preset empirical formula for atmospheric emissivity. Net radiation is calculated using a preset energy balance equation. In this calculation process, the downscaled reconstructed pixel-scale air temperature is used instead of the conventional boundary layer temperature to reduce the thermodynamic errors caused by large-scale meteorological driving. The corresponding calculation formula is as follows:

[0128] R _n =(1-α)×Rswd+σ×ε _s ×(εair×(T _a (x,y) / 100) 4 -(T _s / 100) 4 );

[0129] Among them, R _n The net radiation of the target area is given by α, the surface albedo is given by Rswd, the downwave shortwave radiation is given by σ, the Stefan-Boltzmann constant is given by temperature scaling, εair is given by atmospheric emissivity, and T is given by T. _a (x,y) represents the air temperature at the pixel scale, T _s εair represents the surface temperature, superscripted with 4 (thermodynamic temperature to the fourth power). 100 in this formula is a temperature scaling factor. When measured data is available, the measured Rswd is used; when measured data is lacking, Rswd is calculated based on the solar zenith angle model. εair is calculated based on the actual water vapor pressure at the pixel scale using an empirical formula, where the actual water vapor pressure is determined jointly by the saturated water vapor pressure and the measured relative humidity.

[0130] Step 602: Based on the vegetation cover of the target area, the net radiation is linearly weighted and allocated by combining the preset soil heat flux ratio coefficients of full vegetation cover and bare soil to obtain the soil heat flux.

[0131] Soil heat flux is used to describe the share of energy conducted into the soil interior. A single conversion coefficient is insufficient to accommodate the energy distribution mechanisms of complex land surfaces. In this embodiment, the system extracts the vegetation cover parameter and uses it as the spatial distribution weight for land surface types. The preset ratio coefficients of total vegetation cover and bare soil heat flux are read. The bare soil exposure ratio is obtained by subtracting the vegetation cover from a constant 1. This bare soil exposure ratio is multiplied by the difference between the two soil heat flux ratios, and then the product is added to the total vegetation cover ratio coefficient. Finally, this product is multiplied by the net radiation to output the soil heat flux distribution matrix. The corresponding calculation formula is as follows:

[0132] G _0 =R _n ×(γ _c +(1-cover)×(γ _s -γ _c ));

[0133] Among them, G _0 For soil heat flux, R _n For net radiation, γ _c γ is the soil heat flux ratio under full vegetation cover conditions, where cover is the vegetation cover degree, and γ is the soil heat flux ratio. _s This represents the soil heat flux ratio coefficient under bare soil conditions.

[0134] In one alternative implementation, the specific values ​​of the heat flux ratio coefficient for soil with full vegetation cover and the heat flux ratio coefficient for bare soil can be adaptively adjusted according to the actual surface cover type and climatic conditions.

[0135] In one possible implementation, solving for the instantaneous latent heat flux and evapotranspiration fraction of the target region specifically includes:

[0136] Step 603: Based on the surface energy balance system model, calculate the dry boundary sensible heat flux, the wet boundary sensible heat flux, and the actual sensible heat flux of the target area under dry boundary conditions.

[0137] Specifically, based on the aerodynamic impedance network in the surface energy balance system model, the actual sensible heat flux of the target area is solved by using the surface temperature in the pixel-by-pixel surface characteristic parameters and the air temperature in the pixel-scale meteorological driving parameters.

[0138] After obtaining the fundamental energy components, the surface energy balance system model incorporates the dry-wet boundary theoretical framework. Specifically, a dry boundary condition is defined to represent the boundary condition where soil moisture is depleted, latent heat flux approaches zero, and all available energy is converted into sensible heat flux; a wet boundary condition is defined to represent another boundary condition where water supply is sufficient and sensible heat flux is limited by the potential evapotranspiration rate. Based on aerodynamic impedance and temperature field information, the sensible heat fluxes at the dry and wet boundaries under the two theoretical limiting boundaries are derived respectively.

[0139] For example, by assuming the latent heat flux is 0, all available energy is converted into sensible heat to obtain the dry boundary sensible heat flux, i.e., the dry boundary sensible heat flux equals net radiation minus soil heat flux. Assuming zero surface drag, the latent heat flux corresponding to potential evapotranspiration is calculated using the Penman-Montis formula. Then, the difference between the soil heat flux and the potential latent heat flux is obtained by subtracting the difference from the net radiation, thus yielding the wet boundary sensible heat flux. The Penman-Montis formula is the standard method in this field for estimating potential evapotranspiration and can be conventionally calculated based on pixel-scale meteorological driving parameters and aerodynamic impedance parameters.

[0140] Furthermore, under dry boundary conditions, assuming complete depletion of surface moisture and evapotranspiration approaching zero, the latent heat flux LE _dry The value is 0, based on the surface energy balance equation, the dry boundary sensible heat flux H _dry The calculation formula is:

[0141] H _dry =R _n -G _0 ;

[0142] Among them, R _n For net radiation, G _0 This refers to soil heat flux.

[0143] Under moist boundary conditions, assuming sufficient moisture supply and evapotranspiration reaching its potential maximum rate, the formula for calculating the sensible heat flux at the moist boundary is:

[0144] H _wet =[(R _n -G _0 )×γ _psy -ρa×cp×VPD / rah] / (Δ _s +γ _psy );

[0145] Where, γ _psy ρa is the hygroscopic constant, cp is the specific heat capacity at constant pressure, VPD is the vapor pressure deficit at the pixel scale, rah is the aerodynamic impedance, and Δ is the hygroscopic constant. _s The slope of the saturated vapor pressure-temperature curve is given. Aerodynamic impedance is calculated iteratively based on measured wind speed, momentum roughness length, and thermal roughness length using a pre-defined logarithmic wind profile model and the Moning-Obukhov similarity theory. This constructs the physical envelope space of actual surface heat flux fluctuations.

[0146] Step 604: Dynamically limit the actual sensible heat flux within the constraint range between the dry boundary sensible heat flux and the wet boundary sensible heat flux.

[0147] When dealing with complex underlying surfaces, iterative solutions are prone to generating divergent and outlier solutions. To address this, the system utilizes the sensible heat flux at both the dry and wet boundaries to jointly establish a dynamic numerical constraint range. When the actual sensible heat flux exceeds the upper limit of the dry boundary or falls below the lower limit of the wet boundary, the computation module executes a threshold truncation mechanism, assigning the out-of-limit pixel values ​​to critical boundary values ​​and pulling them back into the constraint range that conforms to the laws of thermodynamics, thereby maintaining the numerical stability of the model for large-scale region inversion.

[0148] Step 605: Calculate relative evaporation based on actual sensible heat flux, dry boundary sensible heat flux, and wet boundary sensible heat flux.

[0149] Relative evaporation, as an intermediate evaluation index of the degree of evapotranspiration potential release, is obtained by processing the distribution location of sensible heat flux. Accordingly, the calculation module extracts the energy difference between the dry boundary sensible heat flux and the actual sensible heat flux, divides it by the difference between the dry boundary sensible heat flux and the wet boundary sensible heat flux, and outputs the relative evaporation coefficient, which represents the degree of actual evapotranspiration approaching the potential maximum evapotranspiration.

[0150] Step 606: Based on relative evaporation and the ratio of sensible heat flux at the wet boundary to that at the dry boundary, the evapotranspiration fraction is calculated, and the instantaneous latent heat flux is obtained based on the energy balance relationship.

[0151] In this embodiment, the evapotranspiration fraction is defined as the proportion of latent heat flux to available surface energy. The ratio of sensible heat flux at the wet boundary to that at the dry boundary is extracted.

[0152] Specifically, the ratio of the sensible heat flux at the wet boundary to the available energy at the surface is subtracted from the constant 1 to obtain the reference value for the evapotranspiration fraction under wet boundary conditions. The relative evaporation coefficient is then multiplied pixel-by-pixel by the reference value for the evapotranspiration fraction under wet boundary conditions to output the final evapotranspiration fraction. The corresponding calculation formula is as follows:

[0153] EF = EFr × (1 - H) _wet / (R _n -G _0 ));

[0154] Where EF is the evapotranspiration fraction, EFr is the relative evaporation, and H _wet For the sensible heat flux at the wet boundary, R _n For net radiation, G _0 This refers to soil heat flux.

[0155] Based on this, and combining the relative evaporation coefficient, the evapotranspiration fraction of each pixel is calculated using the model's energy segmentation operator. The available surface energy field is obtained by subtracting the soil heat flux from the net radiation. Multiplying the evapotranspiration fraction by this available energy field, and then calculating the instantaneous latent heat flux of the target area at the moment of satellite transit according to the surface energy balance equation, this physical quantity characterizes the latent heat energy flow absorbed during water vaporization and can be used for hydrological evapotranspiration conversion.

[0156] Example 7

[0157] This embodiment explains the mechanism by which instantaneous latent heat flux extends to regional diurnal evapotranspiration results.

[0158] In one possible implementation, based on the assumption of relatively stable diurnal evapotranspiration fraction, the instantaneous latent heat flux is combined with diurnal available energy for temporal scaling to output regional diurnal evapotranspiration results, including the following steps:

[0159] Step 701: Using the instantaneous latent heat flux in combination with the preset latent heat of vaporization constant and water density constant, physical flux conversion is performed to obtain the instantaneous evaporation flux at the moment of satellite transit.

[0160] The instantaneous latent heat flux output by the surface energy balance system model belongs to the heat flow field, and its dimension is W / m². 2 To convert this into evapotranspiration in hydrology, the system performs a dimensional conversion between energy and matter. Specifically, it retrieves preset latent heat of vaporization and water density constants. The product of the latent heat of vaporization and water density constants is used as the conversion denominator, and the instantaneous latent heat flux is used as the numerator for division. This operation eliminates the energy unit, deriving and converting the output to represent the instantaneous evapotranspiration flux at the satellite's transit time, representing the total rate of water evaporation and transpiration per unit time. When using the International System of Units (SI), the output dimension of this division operation is m / s, which is then converted to millimeters per hour. The corresponding calculation formula is as follows:

[0161] E=LE / (L _e ×ρ _w )×3600×1000;

[0162] Where E is the instantaneous evapotranspiration flux, LE is the instantaneous latent heat flux, and L... _e ρ is the preset latent heat of vaporization constant. _w 3600 is the preset water density constant, 3600 is the conversion factor from seconds to hours, and 1000 is the conversion factor from meters to millimeters.

[0163] Step 702: Obtain the daily net radiation and daily soil heat flux of the target area based on the radiative transfer model or reanalysis meteorological dataset, and subtract the two to obtain the daily available energy.

[0164] To provide an all-weather energy-driven boundary for timescale extrapolation, the system needs to construct a diurnal energy budget field. Specifically, it obtains the all-day integrated diurnal net radiation matrix provided by a surface radiation balance model or a third-party reanalysis meteorological dataset, and then obtains the diurnal soil heat flux matrix conducted downwards for the corresponding time period. Subtracting the diurnal soil heat flux from the diurnal net radiation yields a difference matrix. This difference matrix represents the total remaining heat energy available for water evapotranspiration in the target area throughout the day, and is packaged and defined as the diurnal available energy output.

[0165] Step 703: Based on the assumption that the evapotranspiration fraction is relatively stable during the day, the evapotranspiration fraction is multiplied pixel by pixel with the daily available energy and converted by combining the latent heat of vaporization constant and the water density constant to obtain the total regional daily evapotranspiration result for the whole day.

[0166] In this embodiment, the assumption of relatively stable diurnal evapotranspiration fraction means that during the daytime period without extreme weather changes, such as from sunrise to sunset, the proportion of latent heat flux to available surface energy remains approximately constant, and its diurnal variation usually does not exceed a preset range of the average value for that period. Since satellite overpass observations can only capture a fleeting technical slice, and water evapotranspiration is affected by solar radiation and wind speed, exhibiting nonlinear fluctuations throughout the day, this embodiment utilizes the assumption of relatively stable diurnal evapotranspiration fraction to overcome this timescale barrier. This physical assumption is based on the observation mechanism that the proportion of latent heat flux to available energy remains approximately constant during the daytime. The instantaneous evapotranspiration fraction is considered as a proportional coefficient for the entire day, and it is multiplied pixel-by-pixel with the diurnal available energy matrix. Based on this, the timescale extension is performed on the heat integral result after multiplication by combining the latent heat of vaporization and the time integral conversion constant, resulting in an extended calculation that outputs a regional diurnal evapotranspiration result representing the total cumulative evapotranspiration throughout the day. The data format is a single-channel two-dimensional millimeter-level grid field.

[0167] Example 8

[0168] To illustrate the engineering usability and commercial value of the above embodiments, this embodiment, based on the same inventive concept as the above method embodiments, provides an application implementation scheme based on multi-channel raster data. In one possible implementation, it further includes a step of data business-oriented export, specifically including:

[0169] Step 801: The normalized vegetation index, surface albedo, vegetation cover, surface emissivity, surface temperature, and regional diurnal evapotranspiration results in the pixel-by-pixel surface feature parameters are converted into multi-channel raster data.

[0170] After a series of algorithms are processed in the cloud, a multi-dimensional physical spatial distribution matrix of the target area can be obtained. Based on this, the system calls the data encapsulation module to combine and overlay the six sets of two-dimensional data structures—normalized vegetation index, surface albedo, vegetation cover, surface emissivity, surface temperature, and regional diurnal evapotranspiration—with a unified geographic reference system and coordinate projection mapping relationship.

[0171] Accordingly, a file format encoding operation is performed to integrate the discrete intermediate features with the final physical output and encode them into a multi-channel raster data file. That is, the multi-channel raster data is serialized and stored using the standard geotagged image file format, preserving floating-point numerical values ​​and spatial topological attributes.

[0172] Step 802: Output the multi-channel raster data to a cloud drive or geographic information system platform to provide indicator parameter support for agricultural irrigation monitoring, water resource assessment, or eco-hydrological analysis.

[0173] To transform cloud-based data output into engineering and business services for end users, this embodiment establishes an application programming interface (API) connection with external storage or a business platform. Encapsulated multi-channel raster data is pushed via network transmission protocols and persistently stored on a cloud hard drive or loaded into the database of a spatial geographic information system platform. Based on this output mechanism, water conservancy and agricultural management departments can extract the required hierarchical data through a visual client.

[0174] Furthermore, after outputting multi-channel raster data to the geographic information system platform, the process also includes:

[0175] Step 803: Based on the output regional daily evapotranspiration results and the pre-stored theoretical water requirement data of the target crop during a specific growth period, assess the water stress status of crops in the target area.

[0176] For commercial applications of agricultural water resource management, this embodiment uses a preset target crop as the monitoring object. After obtaining the regional daily evapotranspiration result matrix, an external agricultural knowledge base is simultaneously searched to obtain reference data on the theoretical water requirement of the target crop at its current preset growth stage.

[0177] For example, taking the rice growing season as an example, the standard water consumption quota for rice within this time window is extracted, and the difference is calculated or proportionally converted with the pixel-level regional daily evapotranspiration result. When the actual calculated total daily evapotranspiration is lower than the warning threshold set by the theoretical water demand, it is determined that the crops at that geographic coordinate have not obtained sufficient water to support transpiration, and the spatial water stress status of crops in the target area is assessed.

[0178] Step 804: Generate and output the corresponding regional agricultural irrigation scheduling strategy based on the water stress state.

[0179] After obtaining the pixel-level distribution matrix of water stress states, spatial clustering analysis is performed on the distributed stress states based on the pre-configured irrigation district water conservancy network topology and water allocation rules. For irrigation units with high water stress values ​​and contiguous areas, the required water compensation is calculated, generating a regional agricultural irrigation scheduling strategy that includes target irrigation area nodes and compensation water amounts. This scheduling strategy is output to the automated irrigation control terminal or management node in the form of control commands or report data, completing the closed loop from satellite remote sensing multi-source data acquisition to quantitative management of field water resources.

[0180] Based on the above embodiments, this application also provides a regional daily-scale evapotranspiration cloud computing system based on multi-source remote sensing data fusion. The system includes a processor coupled to a memory. The memory stores instructions, which are executed by the processor to implement the methods described in the above method embodiments.

[0181] This application also provides an electronic device, including: one or more processors and a memory, wherein the one or more processors are coupled to the memory, the memory is used to store computer program code, the computer program code including computer instructions, and when the one or more processors execute the computer instructions, the electronic device executes the methods executed in the above-described method embodiments.

[0182] This application also provides a chip, including: a processor and a communication interface, wherein the processor calls a computer program or instruction through the communication interface to execute the methods executed in the above-described method embodiments.

[0183] This application also provides a computer program product, which, when run on a computer, executes the methods described in the above-described method embodiments.

[0184] This application also provides a computer-readable storage medium storing a computer program. When the computer program is run on a computer, the computer executes the methods described in the above-described method embodiments.

[0185] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. A regional diurnal evapotranspiration cloud computing method based on multi-source remote sensing data fusion, characterized in that, include: Acquire multi-source remote sensing images, digital elevation models, diurnal radiation data, and station meteorological parameters within the target area; Multi-source remote sensing images are preprocessed using a remote sensing cloud computing platform to obtain a standardized image set; A pre-defined statistical synthesis strategy is used to perform pixel-level statistical synthesis on a standardized image set and then crop it to obtain a fused image of the target area. Based on this, multi-channel feature extraction is performed to obtain pixel-by-pixel surface feature parameters. Based on the digital elevation model, the meteorological parameters of the stations are subjected to three-dimensional spatial correction of topographic undulation dynamics to obtain the meteorological driving parameters at the pixel scale. The pixel-by-pixel surface feature parameters, digital elevation model, and pixel-scale meteorological driving parameters are input into the surface energy balance system model to solve for the instantaneous latent heat flux and evapotranspiration fraction. Based on the assumption that the evapotranspiration fraction is relatively stable during the day, the instantaneous latent heat flux is combined with the day-scale available energy calculated based on day-scale radiation data to extend the time scale and output the day-scale evapotranspiration results for the region. The standardized image set obtained includes: Multi-source remote sensing imagery includes Landsat imagery and Sentinel-2 imagery; Based on the pixel quality assessment band of Landsat imagery, bitwise operations are performed to identify cloud and shadow pixels and perform mask removal. Radiometric calibration is performed on the optical band and brightness temperature conversion is performed on the thermal infrared band to obtain calibrated land imagery. Based on the scene classification band of Sentinel-2 imagery, cloud, shadow, and cirrus cloud pixels are identified and masked removal is performed. Reflectivity scaling is applied to the optical band to obtain cloud-free Sentinel imagery. For the declouded sentinel image in the non-thermal infrared band, a brightness temperature occupancy band with the same band structure is constructed, and its pixel value is uniformly assigned to null value or mask value to obtain the calibrated sentinel image. The standardized land image and the standardized sentinel image are named in the same band. The spatial resolution of the standardized sentinel image is uniformly resampled to be consistent with that of the standardized land image. After performing pixel grid alignment, a standardized image set is obtained. The pixel-scale meteorological driving parameters are obtained, including: Obtain the measured air temperature, measured air pressure, and measured relative humidity from the meteorological parameters of the station; Based on the digital elevation model, the measured air temperature, measured air pressure and measured relative humidity are subjected to pixel-by-pixel spatial downscaling and dynamic correction by the surface elevation undulation. Pixel-scale meteorological driving parameters are generated, including pixel-scale air temperature, pixel-scale air pressure, pixel-scale air density, and pixel-scale water vapor pressure deficit. Based on the digital elevation model, pixel-by-pixel spatial downscaling and dynamic correction are performed on measured air temperature, measured air pressure, and measured relative humidity by analyzing surface elevation undulations, including: The relative elevation difference between the corresponding pixel elevation of the digital elevation model and the preset meteorological station elevation is calculated, and the measured air temperature is corrected pixel by pixel in combination with the preset standard temperature vertical lapse rate to obtain the air temperature at the pixel scale. Based on the relationship between air pressure and altitude, combined with gravitational acceleration, dry air gas constant, and the ratio of air temperature before and after correction, the measured air pressure is corrected for altitude correlation to obtain the air pressure at the pixel scale. By substituting the pixel-scale air temperature and pixel-scale air pressure into the ideal gas law, the physical spatial distribution is extrapolated to obtain the pixel-scale air density. The saturated water vapor pressure of the corresponding pixel is calculated based on the air temperature at the pixel scale, and the actual water vapor deficit is calculated by combining the measured relative humidity, thus obtaining the water vapor pressure deficit at the pixel scale.

2. The method according to claim 1, characterized in that, Pixel-level statistical synthesis of the fused image of the target region is obtained, including: The Landsat image set in the standardized image set will be merged with the Sentinel-2 image set; Using any one of the median composite method, maximum vegetation index composite method, weighted average method, or cloud probability weighted method, temporal pixel-level statistical composite is performed on the merged image set, automatically removing null values ​​or mask values ​​to generate time-series fused images.

3. The method according to claim 2, characterized in that, The cropping process obtained from the fused image of the target region includes: Load the pre-stored target region vector boundary file; By using the target region vector boundary file to perform boundary cropping on the time-series fused image, a target region fused image with regular spatial boundaries is obtained.

4. The method according to claim 1, characterized in that, The pixel-by-pixel surface feature parameters are obtained, including: Extract the near-infrared and red bands from the fused image of the target area and calculate the normalized vegetation index; The reflectance of the land surface is obtained by extracting blue light, red light, near-infrared light and multi-channel short-wave infrared bands in parallel, combining the reflectance of each band and removing the system bias term. The vegetation cover was obtained by spatially normalizing the normalized vegetation index. Based on vegetation cover, the surface emissivity is obtained by weighting and summing using a proportionality coefficient. The effective brightness temperature bands in the fused image of the target area are extracted and nonlinearly converted to obtain the surface temperature. The five parameters are integrated and output as pixel-by-pixel surface feature parameters.

5. The method according to claim 1, characterized in that, It also includes the introduction of surface roughness parameters, specifically including: Extracting normalized vegetation index from pixel-by-pixel surface feature parameters; Calculate the leaf area index and momentum roughness length of the target area based on the normalized vegetation index; The thermal roughness length is derived based on the momentum roughness length and the preset heat transfer impedance parameters. The surface roughness parameters, including leaf area index, momentum roughness length, and thermal roughness length, along with pixel-scale meteorological driving parameters, are synchronously input into the surface energy balance system model.

6. The method according to claim 1, characterized in that, Solving for the instantaneous latent heat flux and evapotranspiration fraction in the target region includes: Based on the surface energy balance system model, the sensible heat flux at the dry boundary under dry boundary conditions, the sensible heat flux at the wet boundary under wet boundary conditions, and the actual sensible heat flux in the target area are extrapolated and calculated. The actual sensible heat flux is dynamically constrained between the dry boundary sensible heat flux and the wet boundary sensible heat flux; Relative evaporation is calculated based on actual sensible heat flux, dry boundary sensible heat flux, and wet boundary sensible heat flux; Based on relative evaporation and the ratio of sensible heat flux at the wet boundary to that at the dry boundary, the evapotranspiration fraction is calculated, and the instantaneous latent heat flux is obtained based on the energy balance relationship. Before deriving the sensible heat flux in the target area under dry and wet boundary conditions, the surface energy balance system model also includes the step of solving for the soil heat flux, specifically: Based on the surface albedo, surface emissivity, surface temperature, and pixel-scale meteorological driving parameters in the pixel-by-pixel surface characteristic parameters, the net radiation of the target area is calculated using the radiative transfer equation. Based on the vegetation cover of the target area, the net radiation is linearly weighted and allocated by combining the preset soil heat flux ratio coefficients of full vegetation cover and bare soil to obtain the soil heat flux.

7. The method according to claim 1, characterized in that, Output the diurnal evapotranspiration results for the region, including: The instantaneous evapotranspiration flux at the moment of satellite passage is obtained by combining the instantaneous latent heat flux with the latent heat of vaporization constant and the water density constant for physical flux conversion. Obtain the daily net radiation and daily soil heat flux of the target area based on the radiative transfer model or reanalysis meteorological dataset, and subtract the two to obtain the daily available energy. Based on the assumption that the evapotranspiration fraction is relatively stable during the day, the evapotranspiration fraction is multiplied pixel by pixel with the daily available energy and converted using the latent heat of vaporization constant and the water density constant to obtain the total regional daily evapotranspiration result for the whole day.

8. The method according to claim 1, characterized in that, It also includes the steps of exporting data for business purposes, specifically including: The normalized vegetation index, surface albedo, vegetation cover, surface emissivity, surface temperature, and regional diurnal evapotranspiration results in the pixel-by-pixel surface feature parameters are converted into multi-channel raster data. Multi-channel raster data is output to cloud drives or geographic information system platforms to provide indicator parameter support for agricultural irrigation monitoring, water resource assessment, or eco-hydrological analysis.

Images