Aerodynamic resistance and canopy resistance compensation and equivalent quantification evaluation method

Abstract

The present application provides an aerodynamic resistance and canopy resistance compensation and equivalent quantitative evaluation method, relates to the micro-meteorology technical field, and the present application acquires meteorological observation data, canopy height and observed latent heat flux to construct a reference parameter point; a two-dimensional grid and a response surface are constructed in a preset resistance search interval; the global equivalent area is screened and the compensation index is calculated by calculating the latent heat flux error; the local vulnerability index is obtained by recalculating the latent heat flux by applying disturbance to the reference parameter point; finally, the quantitative evaluation result is obtained by jointly analyzing the compensation index and the local vulnerability index. The present application effectively removes the fitting interference of parameter coupling, quantifies the risk boundary of the model, and improves the objectivity and reliability of the model simulation under complex meteorological conditions.

Description

Technical Field

[0001] This invention relates to the field of micrometeorology, specifically to a method for compensating for aerodynamic drag and canopy drag, and for quantitative evaluation of equivalence. Background Technology

[0002] In micrometeorological and ecohydrological studies of surface-atmosphere water and heat exchange, aerodynamic drag and canopy drag are two core physical parameters characterizing energy distribution and water transport mechanisms. Aerodynamic drag represents the physical resistance of atmospheric turbulence mixing to the transport process, while canopy drag reflects the comprehensive control of vegetation physiological characteristics on water exchange. With the deep integration of meteorological data IoT sensing technology and land surface process numerical simulation technology, accurately analyzing and characterizing the dynamic interaction patterns of these two drag parameters under different meteorological conditions has become a key technical bottleneck for improving the reliability of simulations in complex environments.

[0003] Regarding the acquisition and application of the aforementioned core drag parameters, existing micrometeorological parameterization techniques typically treat them as independent solution objects or simple unidirectional input variables. For example, in conventional drag parameter analysis schemes, aerodynamic drag is solidified separately using empirical formulas based on empirical observation data, and then substituted as a constant into the physical equations to unidirectionally deduce canopy drag; or deterministic optimization algorithms are used to forcibly search for a set of drag values ​​in the parameter space that minimizes the error of the output model. These existing techniques mainly focus on obtaining static parameter combinations that superficially make the equations hold true.

[0004] However, the aforementioned existing technologies have significant technical shortcomings at the level of analyzing the underlying physical mechanisms. In actual natural environments, there is a strong nonlinear coupling and parameter compensation effect between aerodynamic drag and canopy drag. Existing technologies only pursue the static extraction of single parameters or the numerical fitting of the overall equation, completely ignoring the equivalent interference that may arise from the physical evolution of different combinations of drag parameters. This approach, which ignores the internal coupling and compensation mechanism of parameters, cannot identify the false rationality caused by the mutual cancellation of parameter errors; once the external meteorological conditions deviate from the norm, the parameter vulnerability masked by the fitting spuriousness will be rapidly amplified. Therefore, existing technical solutions lack quantitative means for the equivalent range between aerodynamic drag and canopy drag, and cannot reveal the model risk boundary brought about by parameter compensation from the underlying logic.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: The method for compensating for aerodynamic drag and canopy drag, and for quantitative evaluation of equivalence, includes the following steps: Step 1: Obtain meteorological observation data for the target area, including canopy height and observed latent heat flux; Step 2: Calculate reference aerodynamic drag based on the meteorological observation data and the canopy height. Using the reference aerodynamic drag as a fixed parameter and the observed latent heat flux as the target, calculate the reference canopy drag by inversion using the Penman equation. Then, construct reference parameter points based on the reference aerodynamic drag and the reference canopy drag. Step 3: Construct a two-dimensional grid within the preset aerodynamic drag search range and canopy drag search range; calculate the predicted latent heat flux based on the aerodynamic drag, canopy drag and meteorological observation data of the two-dimensional grid; and construct a two-dimensional response surface through the two-dimensional grid and the predicted latent heat flux. Step 4: In the two-dimensional response surface, calculate the latent heat flux error based on the predicted latent heat flux and the observed latent heat flux, select the set of grid points whose latent heat flux error meets the preset relative error threshold as the global equivalent region, calculate the area of ​​the global equivalent region, and calculate the compensation index based on the total area of ​​the search interval of the two-dimensional grid and the area of ​​the global equivalent region. Step 5: Using the reference parameter point as the center, apply a preset proportion of perturbation to the aerodynamic drag and the canopy drag respectively, recalculate the latent heat flux, and calculate the local vulnerability index based on the latent heat flux, the aerodynamic drag and the canopy drag; Step 6: Perform a joint analysis based on the compensation index and the local vulnerability index to obtain a quantitative assessment result.

[0008] Furthermore, the meteorological observation data includes: net radiation, soil heat flux, air temperature, atmospheric pressure, relative humidity, air density, and wind speed.

[0009] Furthermore, the calculation of the reference aerodynamic drag includes: The zero-plane displacement height and surface roughness length are determined based on the canopy height. Atmospheric stability is calculated based on the air temperature and wind speed, and a stability correction function is calculated based on the atmospheric stability and the Monin-Obukhov similarity theory. The reference aerodynamic drag is calculated based on the logarithmic wind speed profile by combining the wind speed, the zero plane displacement height, the surface roughness length, and the stability correction function. The inversion of the reference canopy drag includes: calculating the saturated vapor pressure difference and the slope of the saturated vapor pressure curve based on the air temperature and the relative humidity; and calculating the humidity constant based on the atmospheric pressure. Substitute the saturated vapor pressure difference, the slope of the saturated vapor pressure curve, the net radiation, the soil heat flux, the humidity constant, the air density, and the reference aerodynamic drag into the Penman equation. Use the reference aerodynamic drag as a fixed parameter, the observed latent heat flux as the target value, and the reference canopy drag as the parameter to be inverted. Substitute the parameter to be inverted into the Penman equation to calculate the corresponding latent heat flux. When the relative error between the latent heat flux and the observed latent heat flux is less than or equal to a preset convergence threshold, the corresponding canopy resistance is determined as the reference canopy resistance.

[0010] Furthermore, the construction of the two-dimensional mesh includes: The preset aerodynamic drag search range and canopy drag search range are discretized according to the preset step size to form aerodynamic drag nodes and canopy drag nodes; A two-dimensional regular mesh is constructed using the aerodynamic drag nodes as the first coordinate axis and the canopy drag nodes as the second coordinate axis. The aerodynamic drag and canopy drag corresponding to each grid node in the two-dimensional grid are used as input parameters. Combined with the meteorological observation data, they are substituted into the Penman equation to calculate the corresponding predicted latent heat flux. A two-dimensional response surface is constructed based on the aerodynamic drag, canopy drag, and predicted latent heat flux corresponding to each grid node in the two-dimensional grid.

[0011] Furthermore, the determination of the global equivalent region includes: Calculate the absolute value of the difference between the predicted latent heat flux and the observed latent heat flux at each grid point; Calculate the ratio of the absolute value to the observed latent heat flux to obtain the latent heat flux error for each grid point; The set of grid points whose latent heat flux error is less than or equal to a preset relative error threshold is divided into a global equivalent region.

[0012] Furthermore, the area corresponding to a single grid cell is determined based on the grid step size of the two-dimensional grid; The number of grid cells contained in the global equivalent region is counted, and the area of ​​the global equivalent region is calculated based on the number of grid cells and the area corresponding to a single grid cell. Contour extraction is performed on the boundary of the global equivalent region to obtain the geometric feature parameters of the boundary of the global equivalent region. The geometric feature parameters include at least one of boundary curvature, aspect ratio, and principal axis direction. The ratio of the global equivalent region area to the total area of ​​the search interval is corrected based on the geometric feature parameters to obtain the compensation index.

[0013] Furthermore, the calculation of the local vulnerability index includes: Using the logarithmic central difference method, positive and negative perturbations of a preset proportion are applied to the aerodynamic drag and canopy drag of the reference parameter point, respectively, and the corresponding latent heat flux is calculated by substituting them back into the Penman equation. Based on the central difference approximation, the first sensitivity index of the latent heat flux to the aerodynamic drag and the second sensitivity index of the latent heat flux to the canopy drag are calculated respectively. Compare the first sensitivity index with the second sensitivity index, and extract the maximum absolute value of the two as the local vulnerability index.

[0014] Furthermore, a meteorological state space is constructed using the net radiation and the saturated water vapor pressure difference as coordinate axes; The meteorological state space is divided into grids according to a preset interval, and the compensation index and the local vulnerability index are mapped as attribute features to the corresponding net radiation and saturated water vapor pressure difference coordinate points in the meteorological state space.

[0015] Furthermore, based on the compensation index and local vulnerability index within each meteorological state grid, regression analysis is performed on each meteorological state grid to determine the parameter constraint capability and model risk boundary corresponding to each meteorological state grid. Based on the parameter constraints and the model risk boundary, the risk level of the corresponding meteorological state is determined, and the compensation index, local vulnerability index, equivalent zone morphology, centerline slope and risk level corresponding to each meteorological state are output as quantitative assessment results.

[0016] Compared with the prior art, the beneficial effects of the present invention are: This invention utilizes the Penman equation to calculate reference aerodynamic drag and reference canopy drag based on meteorological observation data, canopy height, and observed latent heat flux in the target area, constructing reference parameter points. Then, a two-dimensional grid is built within preset aerodynamic drag and canopy drag search intervals, forming a two-dimensional response surface by calculating the predicted latent heat flux. Subsequently, the latent heat flux error is calculated using the predicted and observed latent heat flux within the two-dimensional response surface. A set of grid points satisfying a preset relative error threshold is selected as the global equivalent region, and a compensation index is calculated based on the total area of ​​the two-dimensional grid search interval and the area of ​​the global equivalent region. Existing technologies often employ single-target fitting, neglecting the nonlinear coupling between aerodynamic drag and canopy drag, easily falling into the equivalence illusion caused by the mutual cancellation of parameter errors. This technical solution transforms the difficult-to-observe parameter compensation characteristics into quantifiable global equivalent region area and compensation index, breaking the limitations of traditional methods that seek a single static solution. It accurately eliminates fitting interference caused by parameter coupling, ensuring the true physical meaning of parameter analysis and the reliability of the calculation process. This invention also applies a preset proportion of perturbation to aerodynamic drag and canopy drag, centered on a reference parameter point, to recalculate latent heat flux. Based on the latent heat flux and corresponding drag values, a local vulnerability index is calculated. Finally, a comprehensive quantitative assessment result is obtained by jointly analyzing this local vulnerability index and a compensation index. Traditional analytical methods can only obtain static parameters and cannot detect the resilience of parameter combinations to small fluctuations, leading to serious biases in prediction models under complex weather conditions. This technical solution, by actively applying a preset proportion of perturbation, keenly captures and quantifies the risk level of physical parameters under small deviations, revealing the physical defects hidden beneath the surface-fitted values. Simultaneously, by jointly analyzing the local vulnerability index reflecting resilience and the compensation index reflecting equivalence, it overcomes the limitations of traditional methods relying solely on single residual evaluation, defining the risk boundaries and application limitations of the model in complex evolution processes, and achieving an objective and in-depth comprehensive quantitative assessment. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall method flow of the present invention.

[0018] Figure 2 This is a schematic diagram of the two-dimensional response surface of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0020] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0021] Example: Please see Figures 1-2 The present invention provides a technical solution: The method for compensating for aerodynamic drag and canopy drag, and for quantitative evaluation of equivalence, includes the following steps: Step 1: Obtain meteorological observation data, canopy height, and observed latent heat flux for the target area.

[0022] In this embodiment, the meteorological observation data includes: net radiation, soil heat flux, air temperature, atmospheric pressure, relative humidity, air density, and wind speed.

[0023] Meteorological observation data and latent heat flux of the target area are collected in real time and synchronously using meteorological instruments installed in actual farmland or ecological observation stations. In practical applications, the acquisition frequency of instantaneous data is usually set to a relatively high frequency, such as 10Hz, and then time block averaging is performed according to a preset time step, preferably 30 minutes, to output a continuous half-hour scale meteorological dataset. For canopy height, the value needs to be set according to the actual phenological growth stage of vegetation in the target area. For example, during the heading and jointing stage of a crop, the measured and preset input canopy height is 1.5 meters.

[0024] Among them, net radiation and soil heat flux represent the net available energy driving the phase change evaporation of surface water; the combined calculation of air temperature and relative humidity represents the current dynamic gap for atmospheric absorption of evaporating water vapor; atmospheric pressure and air density represent the basic background parameters used to determine the thermodynamic constants at the current altitude; wind speed represents the direct driving force of the dominant aerodynamic turbulent mixing intensity; canopy height represents the geometric morphological benchmark that determines surface roughness and zero-plane displacement; and observed latent heat flux represents the energy consumed by actual water evaporation under the current environment, serving as the physical reference true value when inverting parameters of the entire algorithm model.

[0025] In the micrometeorological physics scenario, surface evapotranspiration is essentially a material exchange process constrained by both the radiation energy balance equation and the aerodynamic transport equation. Simultaneously acquiring all the aforementioned meteorological elements aims to completely and comprehensively reconstruct the true physical boundary conditions of the vegetation-atmosphere interaction surface at this moment, from both the energy supply and aerodynamic perspectives. Using a preset time averaging step of 30 minutes precisely encompasses the large-scale turbulent eddies containing the main energy transfer components from a fluid dynamics perspective, while effectively filtering out random measurement noise from high-frequency instruments and accurately capturing the intraday evolution characteristics of typical meteorological elements such as afternoon dry heat and nighttime temperature inversion. The meteorological data ensures that the Penman equation derivation remains accurate, providing a data foundation for the subsequent accurate delineation of equivalent zones.

[0026] Step 2: Calculate reference aerodynamic drag based on the meteorological observation data and the canopy height. Using the reference aerodynamic drag as a fixed parameter and the observed latent heat flux as the target, calculate the reference canopy drag by inversion using the Penman equation. Then, construct reference parameter points based on the reference aerodynamic drag and the reference canopy drag.

[0027] In this embodiment, the calculation of the reference aerodynamic drag includes: The zero-plane displacement height and surface roughness length are determined based on the canopy height.

[0028] Canopy height data was extracted, and the zero-plane displacement height and surface roughness length of the vegetation were directly calculated using empirical proportional relationships. Specifically, the zero-plane displacement height was set to two-thirds of the canopy height, and the surface roughness length included momentum surface roughness and thermal surface roughness. The momentum surface roughness was set to 0.123 times the canopy height, and the thermal surface roughness was set to 0.1 times the momentum surface roughness.

[0029] Among them, canopy height represents the average physical height of vegetation in the target area; zero plane displacement height represents the equivalent reference surface displacement of the wind speed profile as it is lifted upward due to the obstruction of the vegetation canopy; and surface roughness length represents the characteristic scale of the drag and friction effect of small surface undulations on near-surface airflow.

[0030] The aforementioned proportional conversion logic has a clear physical meaning in micrometeorology. The geometry of the Earth's surface directly determines the frictional resistance of airflow. Directly reading the actual observed canopy height for conversion can reflect the dynamic growth process of farmland or ecological underlying surfaces, such as the growth of crops from the seedling stage to the jointing stage, in the aerodynamic resistance, avoiding the serious deviation in wind resistance calculations in the later stages of the growing season caused by using fixed constants.

[0031] Atmospheric stability is calculated based on the air temperature and wind speed, and a stability correction function is calculated based on the atmospheric stability and the Monin-Obukhov similarity theory.

[0032] The formula for calculating the stability correction function is: In the formula, For Richardson's number, For gravitational acceleration, the value in this embodiment is taken as [value missing]. , The observation altitude of meteorological instruments. Zero planar displacement height To observe the temperature gradient difference between atmospheric temperature at altitude and surface aerodynamic temperature, The temperature is the temperature under the Kelvin scale. For wind speed, For atmospheric stability.

[0033] Based on this, and using the Monin-Obukhov similarity theory, and according to the calculated sign of the dimensionless atmospheric stability, empirical integral relationships are used to calculate the stability correction functions for momentum and heat respectively: when At that time, the atmosphere exhibits an unstable stratification state: when The atmosphere exhibits a stable or neutral stratification state: In the formula, For intermediate transition parameters, This is the momentum stability correction function. This is the thermal stability correction function.

[0034] Instead of using energy flux to infer stability, this method uses air temperature and wind speed (i.e., the overall Richardson number) to derive atmospheric stability forward, thus eliminating the circular reasoning logic in parameter calculations. Since latent heat flux is the core target variable for subsequent inversion and equivalence assessment, using heat flux data to calculate stratification stability prematurely would lead to physical contamination between aerodynamic drag and the final latent heat target. Determining stability through the pure trade-off between ambient temperature gradient and wind speed allows for a completely independent and objective characterization of the true thermodynamic deformation of the lower atmosphere. When thermal convection is strong, the nonlinear logarithmic integral function provides exponential compensation for transmission efficiency; while when nighttime inversions suppress turbulence, the linear decay function accurately simulates the obstruction and closure of transmission channels, thus providing the cleanest and most independent physical constraint boundary for subsequently eliminating model fitting artifacts caused by meteorological state deviations.

[0035] The reference aerodynamic drag is calculated based on the logarithmic wind speed profile by combining the wind speed, the zero-plane displacement height, the surface roughness length, and the stability correction function.

[0036] Specifically, the formula for calculating aerodynamic drag is as follows: In the formula, For reference to aerodynamic drag, The observation altitude of meteorological instruments. Zero planar displacement height The momentum surface roughness length, The length of the thermal surface roughness. This is the momentum stability correction function. This is a correction function for thermal stability. is the Kalman constant, which is preset to 0.41 in this embodiment. This refers to wind speed.

[0037] Traditional logarithmic wind speed profiles hold true only under absolutely neutral atmospheric conditions, assuming a perfect logarithmic decrease in wind shear drag. However, real-world field water and heat exchange processes are always accompanied by strong thermodynamic fluctuations. By explicitly introducing a subtraction correction factor into the numerator of the formula, this calculation logic mathematically transforms the originally rigid logarithmic decay curve into an adaptive profile that dynamically bends with atmospheric warm and cold stratification.

[0038] By utilizing observational parameters and standardized surface morphology parameters, a drag anchor point was constructed that closely approximates the physical boundary of real turbulence. Only by completely severing and eliminating the computational bias inherent in aerodynamic drag due to abrupt environmental changes can this anchor point be used as a fixed reference to invert canopy drag without forcing the canopy drag to passively subtract this error. This fundamentally prevents cross-contamination of errors between the two drag parameters, thus laying an indisputable physical foundation for subsequently constructing an accurate two-dimensional response surface and an equivalent space for truly quantifying parameters with different names within a standardized parameter plane.

[0039] The inversion of the reference canopy drag includes: calculating the saturated vapor pressure difference and the slope of the saturated vapor pressure curve based on the air temperature and the relative humidity; and calculating the humidity constant based on the atmospheric pressure.

[0040] Based on the collected air temperature and relative humidity, the slope of the saturated vapor pressure difference versus saturated vapor pressure curve was extracted using the Tetens empirical formula, and the humidity constant was calculated in conjunction with atmospheric pressure. The calculation formula is as follows: In the formula, The saturated vapor pressure, The temperature is the temperature under the Kelvin scale. For saturated water vapor pressure difference, Relative humidity, The slope of the water vapor pressure curve. The humidity constant is The specific heat capacity of air at constant pressure is preset to [value] in this embodiment. ; Atmospheric pressure The ratio of the molecular weight of water vapor to that of dry air is preset to a constant of 0.622; The latent heat of water vaporization is preset to be .

[0041] Substitute the saturated vapor pressure difference, the slope of the saturated vapor pressure curve, the net radiation, the soil heat flux, the humidity constant, the air density, and the reference aerodynamic drag into the Penman equation. Using the reference aerodynamic drag as a fixed parameter and the observed latent heat flux as the target value, and the reference canopy drag as the parameter to be inverted, substitute the parameter to be inverted into the Penman equation to calculate the corresponding latent heat flux.

[0042] Then, a constrained inversion of the reference canopy drag is performed. The intermediate variables and the reference aerodynamic drag are substituted into the Penman equation with latent heat flux as the dependent variable for iterative solution: In the formula, For latent heat flux, The slope of the water vapor pressure curve. Net radiation, For soil heat flux, air density, The specific heat capacity of air at constant pressure. For saturated water vapor pressure difference, For reference to aerodynamic drag, The humidity constant is These are the parameters to be inverted that are continuously updated within the interval.

[0043] When the relative error between the latent heat flux and the observed latent heat flux is less than or equal to a preset convergence threshold, the corresponding canopy resistance is determined as the reference canopy resistance.

[0044] Using a one-dimensional numerical search algorithm, within a pre-defined physically reasonable range of canopy resistance, for example... The internal dynamics generate and update the resistance value to be inverted, and calculate the relative residual error after each update: In the formula, To observe latent heat flux, For latent heat flux, This represents the relative residual error.

[0045] When the relative error is determined to be less than or equal to the preset convergence threshold, the preset convergence threshold is set to 0.001, the iteration is terminated immediately, and the resistance value to be inverted at this time, which is in the denominator of the equation, is forcibly output and confirmed as the unique reference canopy resistance.

[0046] The saturated vapor pressure difference quantifies the absolute dynamic hunger of external dry air for absorbing moisture, while the slope of the curve and the humidity constant constitute the thermodynamic balance of whether net surface radiation energy is used to heat the air or evaporate moisture. The precise stripping of these environmental parameters provides an absolutely clean background for calculating vegetation stomatal behavior.

[0047] The Penman model is essentially a highly coupled nonlinear energy equation. Without constraints, directly seeking the parameter combination that makes the equation valid can easily lead to the illusion of error swallowing and substitution between the two drag parameters. This control algorithm, by locking one end and approximating the true value, uses a reference aerodynamic drag derived from pure meteorological fluid dynamics as a static benchmark and utilizes the latent heat flux observed in the field as the sole, uncompromising target. This numerical inversion method forcibly eliminates the interference of multiple solution spaces, forcing the canopy drag to reveal its most realistic physiological drag value under the current specific radiation and humidity conditions. This not only gives the reference parameter point strong real-world physical significance but also provides an absolute coordinate origin without any offset for subsequent two-dimensional mesh deployment.

[0048] Step 3: Construct a two-dimensional grid within the preset aerodynamic drag search range and canopy drag search range. Calculate the predicted latent heat flux based on the aerodynamic drag, canopy drag, and meteorological observation data of the two-dimensional grid. Construct a two-dimensional response surface using the two-dimensional grid and the predicted latent heat flux.

[0049] In this embodiment, the construction of the two-dimensional mesh includes: The preset aerodynamic drag search range and canopy drag search range are discretized according to the preset step size to form aerodynamic drag nodes and canopy drag nodes; A two-dimensional regular mesh is constructed by using the aerodynamic drag nodes as the first coordinate axis and the canopy drag nodes as the second coordinate axis.

[0050] The preset aerodynamic drag search range is set to According to the preset step size One-dimensional discretization is performed to extract the aerodynamic drag node sequence; the preset canopy drag search interval is set as follows. One-dimensional discretization is performed according to a preset step size to extract the canopy drag node sequence. Then, using the aerodynamic drag nodes as the first coordinate axis and the canopy drag nodes as the second coordinate axis, a two-dimensional regular grid matrix with a fixed resolution is generated in memory through Cartesian product operations of all permutations. Each coordinate intersection point in this matrix represents a specific set of coordinates. Parameter combinations.

[0051] At the underlying mechanism of micrometeorological parameterized inversion, traditional sensitivity analysis is often limited to single-parameter perturbations, i.e., fixing one resistance and fine-tuning another to observe output changes. This dimensionality reduction completely severs the real nonlinear coupling relationship between the two resistances. By setting a clear search interval and discretizing to construct a two-dimensional regular grid, we are essentially laying out an absolutely flat base map within a unified parameter plane for subsequent investigation of multiple solutions to the physical equations. This makes the evolution of parameters no longer a scattered, single-line trial and error, but a gridded spatial traversal under a global perspective.

[0052] The aerodynamic drag and canopy drag corresponding to each grid node in the two-dimensional grid are used as input parameters. Combined with the meteorological observation data, they are substituted into the Penman equation to calculate the corresponding predicted latent heat flux.

[0053] The data of each pair of coordinate nodes in the two-dimensional grid matrix, namely aerodynamic drag and canopy drag, are substituted into the Penman equation simultaneously with meteorological observation data. The predicted latent heat flux at the coordinates of that grid point is explicitly solved and output. The formula for calculating the predicted latent heat flux is as follows: In the formula, For the first The aerodynamic drag node and the first Predicted latent heat flux at each canopy resistance node The slope of the water vapor pressure curve. Net radiation, For soil heat flux, air density, The specific heat capacity of air at constant pressure. For saturated water vapor pressure difference, The humidity constant is For the first One aerodynamic drag, For the first Individual canopy resistance, and For coordinate index.

[0054] A two-dimensional response surface is constructed based on the aerodynamic drag, canopy drag, and predicted latent heat flux corresponding to each grid node in the two-dimensional grid.

[0055] After all grid points have been traversed and calculated, the horizontal axis nodes, vertical axis nodes, and corresponding predicted latent heat flux values ​​are mapped to the three-dimensional geometric space, thus constructing a continuous surface with the parameter plane as the base and the latent heat flux value as the height, i.e., a two-dimensional response surface.

[0056] At the level of physical equations and control theory, the denominator of the Penman equation simultaneously includes the sum and ratio of aerodynamic drag and canopy drag. This inherent mathematical structure determines its strong compromising characteristic for both types of drag. Explicitly representing the predicted latent heat flux as a two-dimensional response surface completely abandons the local perspective that focuses only on a single residual value. This operation concretizes the originally obscure and abstract algebraic calculation of parameters into a visually observable topographic relief map. In this three-dimensional topographic map, if the predicted latent heat flux corresponding to different drag combinations are almost at the same level, this clearly and rigorously exposes the strength of mutual compensation between physical parameters from both geometric and numerical perspectives. The successful construction of this response surface is the core prerequisite for achieving computable, visualized, and quantifiable output for problems with multiple solutions and non-uniqueness.

[0057] Step 4: In the two-dimensional response surface, calculate the latent heat flux error based on the predicted latent heat flux and the observed latent heat flux, select the set of grid points whose latent heat flux error meets the preset relative error threshold as the global equivalent region, calculate the area of ​​the global equivalent region, and calculate the compensation index based on the total area of ​​the search interval of the two-dimensional grid and the area of ​​the global equivalent region.

[0058] In this embodiment, the determination of the global equivalent region includes: Calculate the absolute value of the difference between the predicted latent heat flux and the observed latent heat flux at each grid point; The ratio of the absolute value to the observed latent heat flux is calculated to obtain the latent heat flux error corresponding to each grid point.

[0059] For each grid node in the two-dimensional response surface, extract the predicted latent heat flux corresponding to that node. Subtract the observed latent heat flux from the predicted latent heat flux and take the absolute value to obtain the absolute deviation of that node. Then, divide the absolute deviation by the value of the observed latent heat flux to calculate a dimensionless percentage value, which is determined as the latent heat flux error of the corresponding grid point.

[0060] The absolute error value alone loses its physical meaning for intertemporal comparisons due to the varying base values ​​of latent heat flux under different seasons and weather patterns. For example, latent heat flux is extremely high on sunny days in summer and extremely low on cloudy days in winter. By performing dimensionless processing, the prediction biases of the equations under all complex meteorological conditions can be uniformly mapped to a standard relative error dimension. This conventional but core processing eliminates the interference caused by fluctuations in the environmental energy base, ensuring that the assessment has an objective and unique comparative benchmark under different weather patterns such as sunny, cloudy, and overcast days.

[0061] The set of grid points whose latent heat flux error is less than or equal to a preset relative error threshold is divided into a global equivalent region.

[0062] In this embodiment, considering the typical measurement uncertainty of real field observation instruments, the preset relative error threshold is preferably set to 0.05. All parameter coordinate points in the two-dimensional response surface that meet this tolerance condition are extracted and clustered to form a specific spatial set, defined as the global equivalent region. The core set selection logic is as follows: In the formula, For the global equivalent region, For the first One aerodynamic drag, For the first Individual canopy resistance, For the first The aerodynamic drag node and the first Predicted latent heat flux at each canopy resistance node To observe latent heat flux, The preset relative error threshold is preferably 0.05.

[0063] Because aerodynamic drag and canopy drag have a strong nonlinear compensation effect on the energy distribution in the evaporation equation, they inevitably lead to a large number of drastically different combinations of drag parameters. The latent heat flux results calculated from these combinations are extremely similar, making the models appear accurate from external observations. The essence of constructing this set using the above formulas is not to find a single optimal solution, but rather to do the opposite: by setting a physically acceptable reasonable error boundary—a 5% tolerance—to directly delineate a multi-solution space on a vast parameter plain. This clearly indicates that under the current specific meteorological conditions, the model's parameters are highly susceptible to mutual compromise and masking of equivalent zones. This not only transforms the abstract physical phenomenon of heterogeneous parameters with the same effect into a tangible and measurable mathematical entity, but also completely breaks through the technical blind spot of conventional evaluations that only focus on the highest fitting accuracy. It lays a rigorous data structure foundation for further quantitative evaluation of the severity of this compensation effect, i.e., the area of ​​this equivalent zone.

[0064] In this embodiment, the area corresponding to a single grid cell is determined based on the grid step size of the two-dimensional grid.

[0065] Extract the preset grid step size, which is set to [value] in this embodiment. Multiplying the two directly yields the physical area of ​​a single grid cell, i.e. .

[0066] At the physical mechanism level, the parameter space of the real environment is in an infinitely continuous state, making it difficult to measure directly. Rasterizing it with a uniform minimum fineness and assigning it definite discrete area primitives is to rigorously transform the subsequent multi-solution fuzzy statistics of parameters into geometric measurement operations with clear dimensions, providing a clear and proportional physical base map for the error boundary of the equivalent region.

[0067] The number of grid cells contained in the global equivalent region is counted, and the area of ​​the global equivalent region is calculated based on the number of grid cells and the area corresponding to each individual grid cell.

[0068] Traverse the two-dimensional response surface, count the total number of grid points falling within the equivalent region set, and then calculate the absolute spatial area by discrete summation. The formula for calculating the global equivalent region area is: In the formula, The area of ​​the global equivalent region. For the number of grid cells, This represents the grid step size.

[0069] Essentially, it performs area summation on a two-dimensional parametric plane. It abandons computationally expensive partial differential operations and cleverly transforms the phenomenon of homonyms in nonlinear equations, which is difficult to analyze directly, into pixel-level surface accumulation of discrete grids, thus determining the absolute spatial scale where parameters overlap.

[0070] The boundary of the global equivalent region is contour extracted to obtain the geometric feature parameters of the boundary of the global equivalent region. The geometric feature parameters include at least one of boundary curvature, aspect ratio, and principal axis direction.

[0071] The global equivalent region in the two-dimensional response surface is treated as a digital image matrix, and the outer closed contour of this region is extracted using a morphological edge tracing algorithm. In this embodiment, the Suzuki85 algorithm is preferably used to calculate the principal axis length and secondary axis length of the circumscribed geometry of the closed contour, and the ratio of the principal axis length to the secondary axis length is defined as the aspect ratio, which is used as the extracted geometric feature parameter.

[0072] The formula for calculating geometric characteristic parameters is: In the formula, Geometric feature parameters, Main axis length, This is the length of the secondary axis.

[0073] Knowing only the area value cannot reveal the directional patterns of error generation. The geometry of the global equivalent region directly reflects the physical properties of the compensation direction between parameters. For example, a region that tends to be circular indicates that the contributions of the two drag parameters to the model deviation are independent and equal; while a narrow and elongated region with a very large aspect ratio strongly suggests that there is an extremely sensitive linear substitution relationship between the two parameters in a specific slope direction.

[0074] The ratio of the global equivalent region area to the total area of ​​the search interval is corrected based on the geometric feature parameters to obtain the compensation index.

[0075] The formula for calculating the compensation index is: In the formula, As a compensation indicator, The area of ​​the global equivalent region. The total area of ​​the search interval. These are geometric characteristic parameters.

[0076] By dividing by the total area of ​​the search domain, the compensation phenomenon is strictly mapped to a standardized range of 0 to 1, eliminating the magnitude differences between different measurement stations and meteorological baselines, and achieving absolute comparability across stations. Simultaneously, a geometric correction coefficient is introduced, ensuring that the final evaluation index encompasses both the range of parameter compensation occurrence and the drastic tendency resulting from it. This index accurately characterizes the identifiability of the Penman equation's internal parameters under current climatic conditions, providing highly penetrating numerical evidence for assessing the engineering reliability of the inversion results.

[0077] Step 5: Using the reference parameter point as the center, apply a preset proportion of perturbation to the aerodynamic drag and canopy drag respectively, recalculate the latent heat flux, and calculate the local vulnerability index based on the latent heat flux, the aerodynamic drag and the canopy drag.

[0078] In this embodiment, the calculation of the local vulnerability index includes: Using the logarithmic central difference method, positive and negative perturbations of a preset proportion are applied to the aerodynamic drag and canopy drag of the reference parameter point, respectively, and the corresponding latent heat flux is calculated by substituting them back into the Penman equation. Based on the central difference approximation, the first sensitivity index of the latent heat flux to the aerodynamic drag and the second sensitivity index of the latent heat flux to the canopy drag are calculated respectively.

[0079] Extract reference parameter points, namely reference aerodynamic drag and reference canopy drag, and set the perturbation step size in logarithmic space. In this embodiment, a preset step size is used. This is equivalent to applying a relative proportional fluctuation of approximately 1%. Keeping one resistance parameter fixed at a reference value, and applying positive and negative exponential scaling to another resistance parameter, the resulting four new parameter combinations are substituted into the Penman equation to obtain the corresponding latent heat flux. The first-order partial derivative is then approximated using the central difference method. The formula for calculating the sensitivity index is: In the formula, As the primary sensitivity indicator, As the second sensitivity indicator, The logarithm is the base of the natural constant. To predict latent heat flux, For reference to aerodynamic drag, For reference to canopy resistance, and These are the preset perturbation step sizes. The resulting positive magnification factor and negative reduction factor This is the preset perturbation step size.

[0080] Introducing logarithmic space and exponential perturbations essentially transforms absolute deviation into a constant percentage rate of change test, ensuring that the sensitivity of the two physical dimensions is measured under a unified and fair scale. Furthermore, the use of central difference instead of unidirectional difference, probing only in the positive or negative direction, is because the Penman equation is a fractional function with inherent nonlinear characteristics. Unidirectional probing is easily misled by local extrema or curve skewness, while central difference simultaneously pulls symmetrically to both sides, effectively offsetting the bias caused by nonlinear curvature using second-order truncation accuracy, thus accurately and realistically capturing the absolute normal slope of the reference point in the current state.

[0081] Compare the first sensitivity index with the second sensitivity index, and extract the maximum absolute value of the two as the local vulnerability index.

[0082] The absolute values ​​of the first and second sensitivity indices are positiveed respectively. Then, the two absolute values ​​are compared, and the highest value is directly assigned to the local vulnerability index.

[0083] When meteorological conditions shift or sensor measurements introduce minor errors, the collapse of prediction results is often not determined by the most stable parameter, but by the most uncontrolled and drastically magnified weakness. Instead of averaging or weighting, directly identifying the direction with the largest absolute value allows for a highly sensitive and conservative identification of the worst-case local sensitivity state the model is in under the current, given meteorological conditions. This provides the most direct and safe warning threshold to prevent sudden amplification of prediction bias.

[0084] Step 6: Perform a joint analysis based on the compensation index and the local vulnerability index to obtain a quantitative assessment result.

[0085] In this embodiment, the meteorological state space is constructed using the net radiation and the saturated water vapor pressure difference as coordinate axes; The meteorological state space is divided into grids according to a preset interval, and the compensation index and the local vulnerability index are mapped as attribute features to the corresponding net radiation and saturated water vapor pressure difference coordinate points in the meteorological state space.

[0086] Specifically, the time series of net radiation and saturated water vapor pressure difference are extracted from meteorological observation data. Net radiation is set as... The axis is preset with its physical data coverage range, such as and according to the preset interval Perform grid-based segmentation; simultaneously, set the saturated water vapor pressure difference to... The axis is preset with its physical data coverage range, such as and according to the preset interval The data is segmented into grids. The intersection of the two coordinate axes creates a meteorological state space matrix composed of a discretized two-dimensional bin array.

[0087] Net radiation represents the upper limit of total available energy that can be mobilized by a phase transition in the environment, while the saturated vapor pressure difference directly quantifies the dynamic thirst of the atmosphere in accommodating additional moisture. Orthogonally mapping these two decisive factors allows for a rigorous and complete delineation of typical weather patterns in the real world. For example, high net radiation coupled with a high vapor pressure difference represents intense sunny and hot stress; low net radiation coupled with a low vapor pressure difference directly corresponds to a period of cloudy, rainy, cold, and mild weather. Binning and gridding essentially transforms the continuously evolving and chaotic natural meteorological flow over time into a discrete, controllable physical testing environment container.

[0088] For any meteorological observation sample at any time step, the actual net radiation and actual saturated vapor pressure difference at that moment are read. Using a rounding-down indexing algorithm, the specific grid coordinates of the sample within the current meteorological state space matrix are determined. Subsequently, compensation indicators and local vulnerability indicators are treated as a set of attribute features and directly written into these specific grid coordinates. As massive amounts of historical time-series data spanning months or even years are continuously input and traversed, each meteorological state grid will gradually accumulate and fill with multiple sets of parameter compensation and vulnerability assessment sample data belonging to that specific weather pattern.

[0089] By inputting compensation and vulnerability indices into this meteorological coordinate system, we can not only determine whether there are errors in the model's inversion parameters, but also accurately determine whether the model is more likely to lose parameter uniqueness under extreme drought stress or strong radiation. This lays a highly interpretive physical diagnostic foundation for clarifying and outlining the specific weather conditions under which the model will fail.

[0090] In this embodiment, regression analysis is performed on each meteorological state grid based on the compensation index and local vulnerability index within each meteorological state grid to determine the parameter constraint capability and model risk boundary corresponding to each meteorological state grid.

[0091] For any discrete grid in the meteorological state space, representing a specific combination of net radiation and saturated vapor pressure difference, all historical time-step sample points falling within that grid are extracted. The compensation index corresponding to each sample point is used as the independent variable, and its corresponding local vulnerability index as the dependent variable, constructing a two-dimensional sample scatter set. Subsequently, a linear regression analysis based on the least squares method is performed on this sample scatter set to obtain the slope, intercept, and residual standard deviation of the regression equation. Finally, based on the statistical benchmark parameters generated by the regression analysis, the parameter constraint capability and model risk boundary corresponding to this meteorological state grid are quantitatively calculated. The formulas for calculating the parameter constraint capability and model risk boundary are: In the formula, The first in the meteorological state grid A local vulnerability indicator, The slope The first in the meteorological state grid One compensation indicator, The intercept is... For the first The regression residuals of each sample point For parameter constraint capability, This is the compensation index for all sample points in the meteorological state grid. For the model risk boundary, The 95th percentile of the compensation index for all sample points within the meteorological state grid; This represents the standard deviation of the residuals.

[0092] The compensation index, in its physical geometry, measures the horizontal area of ​​the global multi-solution space, i.e., how many distinct combinations of spurious parameters can be pieced together to produce the same observed flux; while the local vulnerability index measures the vertical collapse slope of the parameter combination in the worst-case direction, i.e., the rate at which these parameters lose accuracy in the face of minor meteorological fluctuations. Under a specific weather pattern, these two do not exist in isolation, but rather have a deep physical coupling relationship.

[0093] Based on the parameter constraints derived from the regression equation, a negative exponential function is cleverly used to transform the combined product of the slope and the average compensated area into an intuitive score. When the regression slope is extremely steep and the average equivalent area is extremely large, this score rapidly approaches 0. In a micrometeorological sense, this directly indicates that the current radiation or drought environment has completely lost its physical constraint on the resistance parameters of the underlying layer, allowing the parameters to wander freely in a vast error space.

[0094] The calculated model risk boundary is determined using the statistical confidence upper bound theory, which involves adding 1.96 times the residual standard deviation to the regression prediction extreme value to forcibly define a worst-case envelope in mathematical space. This not only covers the average deviation trend under this weather pattern but also fully absorbs the unexplained residuals caused by high-frequency turbulent random fluctuations in nature. Therefore, it provides an absolutely conservative quantitative red line with an extremely high engineering safety margin for subsequent determination of whether the model is distorted.

[0095] Based on the parameter constraints and the model risk boundary, the risk level of the corresponding meteorological state is determined, and the compensation index, local vulnerability index, equivalent zone morphology, centerline slope and risk level corresponding to each meteorological state are output as quantitative assessment results.

[0096] The parameter constraint capability and model risk boundary are introduced into the piecewise threshold function for state delimitation and risk classification. The judgment logic operation is as follows: In the formula, Risk level, For parameter constraint capability, For the model risk boundary, 0.6 and 0.3 represent the preset high and low constraint capability dividing constant thresholds in this embodiment; 0.15 and 0.4 represent the preset low and high model boundary risk dividing constant thresholds. In actual implementation, the above thresholds are pre-calibrated by the historical long-term observation data of farmland micro-meteorology in the target area.

[0097] This embodiment preferably uses principal component analysis (PCA) to extract the slope of the principal axis centerline. PCA uses orthogonal transformations to find the direction that maximizes the data variance, thus more realistically reflecting the core framework of a narrow point set in two-dimensional space. The specific calculation steps are as follows: Extract the set of coordinates of all grid points falling within the global equivalent region in the two-dimensional response surface. First, calculate the centroid coordinates of the point set. And translate all coordinate points to a new coordinate system with the centroid as the origin: Based on the centralized set of coordinate points, construct covariance matrix : Perform eigenvalue decomposition on the covariance matrix and extract the first principal component eigenvector corresponding to the largest eigenvalue. .

[0098] The direction of this feature vector is the principal axis centerline of the equivalent region. The ratio of the x and y components of this feature vector is directly extracted as the slope of the centerline. In the formula, The slope of the principal axis centerline. This represents the total number of grid points within the global equivalent region. and These represent the directional components of the first principal component eigenvector in the canopy drag dimension and the aerodynamic drag dimension, respectively. For coordinate index.

[0099] The final risk classification level, along with the baseline compensation index and local vulnerability index calculated under the specific weather conditions, and the geometric morphology of the equivalent zone extracted based on the boundary contour, are combined with the slope of the principal axis centerline of the narrow equivalent zone extracted using principal component analysis. These five core indicators are then packaged as a complete digital feature map, directly output and mapped to the user interface or the underlying database of the agricultural information management software.

[0100] When faced with special weather patterns such as sunny, hot, and dry conditions or extreme prolonged periods of overcast and rainy weather, it is easy to fall into the trap of a perfectly fitted surface equation but distorted internal parameters. By forcibly delineating risk levels using multiple conditional thresholds, the dangerous inertia of forcibly extrapolating predictions regardless of their accuracy is directly cut off, providing a safe bottom line for irrigation decisions.

[0101] More importantly, the remaining four output indicators constitute a highly penetrating physical origin-tracing matrix: the compensation and vulnerability indicators reveal whether the error presents global ambiguity or local extreme sensitivity; the equivalence zone morphology indicates whether parameter coupling is random diffusion or directional compromise; and the centerline slope pierces the mathematical disguise of the Penman equation. Under this specific meteorological condition, when aerodynamic drag is overestimated due to instrument errors or gust disturbances, canopy drag will inevitably be proportionally underestimated along this slope path, thus deceiving the equation to maintain the surface conservation of latent heat flux. This panoramic quantitative output, from macro-risk warnings to micro-alternative path analysis, completely breaks down the industry barrier of judging the quality of models solely based on statistical errors, achieving a thorough analysis of the parameter equivalence and substitution mechanism.

[0102] Table 1: Spatial Quantitative Assessment Results of Meteorological States As can be seen from the table, by traversing and evaluating the panoramic meteorological conditions, implicit errors are transformed into visualized physical characteristics, which comprehensively improves the engineering reliability of micro-meteorological simulation under complex meteorological conditions. Specifically, under extreme weather conditions such as extreme desertification, drought, and heat waves, this scheme keenly captures the characteristic of a sharp decline in parameter constraint capability due to the loss of physical constraints on parameters by the environment. By forcibly triggering high-risk judgments through multiple conditional thresholds, it directly cuts off the dangerous inertia of forced extrapolation predictions, achieving precise safety risk blocking. At the same time, this scheme deeply understands and quantifies the heterogeneous parameter same-effect mechanism. Utilizing the elongated equivalent zone morphology and characteristic spectral data such as the principal axis centerline slope of up to 0.85, it pierces the mathematical disguise of the Penman equation and thoroughly analyzes the extremely sensitive directional linear substitution relationship between aerodynamic drag and canopy drag. Furthermore, with the help of this quantitative evaluation result, the scheme can accurately remove fitting artifacts and, in reverse, guide the locking of the optimal high-reliability parameter inversion window under low-risk conditions such as mild and sunny spring mornings when compensation indicators are at extremely low levels and the equivalent zone tends to be circular. This allows for the acquisition of reference parameter points with absolute physical truth, completely ensuring the purity of the underlying physical environment.

[0103] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0104] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0105] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0106] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for compensating for aerodynamic drag and canopy drag, and for quantitatively evaluating their equivalence, characterized in that, The specific steps include: Step 1: Obtain meteorological observation data for the target area, including canopy height and observed latent heat flux; Step 2: Calculate reference aerodynamic drag based on the meteorological observation data and the canopy height. Using the reference aerodynamic drag as a fixed parameter and the observed latent heat flux as the target, calculate the reference canopy drag by inversion using the Penman equation. Then, construct reference parameter points based on the reference aerodynamic drag and the reference canopy drag. Step 3: Construct a two-dimensional grid within the preset aerodynamic drag search range and canopy drag search range; calculate the predicted latent heat flux based on the aerodynamic drag, canopy drag and meteorological observation data of the two-dimensional grid; and construct a two-dimensional response surface through the two-dimensional grid and the predicted latent heat flux. Step 4: In the two-dimensional response surface, calculate the latent heat flux error based on the predicted latent heat flux and the observed latent heat flux, select the set of grid points whose latent heat flux error meets the preset relative error threshold as the global equivalent region, calculate the area of ​​the global equivalent region, and calculate the compensation index based on the total area of ​​the search interval of the two-dimensional grid and the area of ​​the global equivalent region. Step 5: Using the reference parameter point as the center, apply a preset proportion of perturbation to the aerodynamic drag and the canopy drag respectively, recalculate the latent heat flux, and calculate the local vulnerability index based on the latent heat flux, the aerodynamic drag and the canopy drag; Step 6: Perform a joint analysis based on the compensation index and the local vulnerability index to obtain a quantitative assessment result.

2. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 1, characterized in that: The meteorological observation data include: net radiation, soil heat flux, air temperature, atmospheric pressure, relative humidity, air density, and wind speed.

3. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 2, characterized in that: The calculation of the reference aerodynamic drag includes: The zero-plane displacement height and surface roughness length are determined based on the canopy height. Atmospheric stability is calculated based on the air temperature and wind speed, and a stability correction function is calculated based on the atmospheric stability and the Monin-Obukhov similarity theory. The reference aerodynamic drag is calculated based on the logarithmic wind speed profile by combining the wind speed, the zero plane displacement height, the surface roughness length, and the stability correction function. The inversion of the reference canopy drag includes: calculating the saturated vapor pressure difference and the slope of the saturated vapor pressure curve based on the air temperature and the relative humidity; and calculating the humidity constant based on the atmospheric pressure. Substitute the saturated vapor pressure difference, the slope of the saturated vapor pressure curve, the net radiation, the soil heat flux, the humidity constant, the air density, and the reference aerodynamic drag into the Penman equation. Use the reference aerodynamic drag as a fixed parameter, the observed latent heat flux as the target value, and the reference canopy drag as the parameter to be inverted. Substitute the parameter to be inverted into the Penman equation to calculate the corresponding latent heat flux. When the relative error between the latent heat flux and the observed latent heat flux is less than or equal to a preset convergence threshold, the corresponding canopy resistance is determined as the reference canopy resistance.

4. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 1, characterized in that: The construction of the two-dimensional mesh includes: The preset aerodynamic drag search range and canopy drag search range are discretized according to the preset step size to form aerodynamic drag nodes and canopy drag nodes; A two-dimensional regular mesh is constructed using the aerodynamic drag nodes as the first coordinate axis and the canopy drag nodes as the second coordinate axis. The aerodynamic drag and canopy drag corresponding to each grid node in the two-dimensional grid are used as input parameters. Combined with the meteorological observation data, they are substituted into the Penman equation to calculate the corresponding predicted latent heat flux. A two-dimensional response surface is constructed based on the aerodynamic drag, canopy drag, and predicted latent heat flux corresponding to each grid node in the two-dimensional grid.

5. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 1, characterized in that: The determination of the global equivalent region includes: Calculate the absolute value of the difference between the predicted latent heat flux and the observed latent heat flux at each grid point; Calculate the ratio of the absolute value to the observed latent heat flux to obtain the latent heat flux error for each grid point; The set of grid points whose latent heat flux error is less than or equal to a preset relative error threshold is divided into a global equivalent region.

6. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 5, characterized in that: The area corresponding to a single grid cell is determined based on the grid step size of the two-dimensional grid; The number of grid cells contained in the global equivalent region is counted, and the area of ​​the global equivalent region is calculated based on the number of grid cells and the area corresponding to a single grid cell. Contour extraction is performed on the boundary of the global equivalent region to obtain the geometric feature parameters of the boundary of the global equivalent region. The geometric feature parameters include at least one of boundary curvature, aspect ratio, and principal axis direction. The ratio of the global equivalent region area to the total area of ​​the search interval is corrected based on the geometric feature parameters to obtain the compensation index.

7. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 1, characterized in that: The calculation of the local vulnerability index includes: Using the logarithmic central difference method, positive and negative perturbations of a preset proportion are applied to the aerodynamic drag and canopy drag of the reference parameter point, respectively, and the corresponding latent heat flux is calculated by substituting them back into the Penman equation. Based on the central difference approximation, the first sensitivity index of the latent heat flux to the aerodynamic drag and the second sensitivity index of the latent heat flux to the canopy drag are calculated respectively. Compare the first sensitivity index with the second sensitivity index, and extract the maximum absolute value of the two as the local vulnerability index.

8. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 3, characterized in that: A meteorological state space is constructed using the net radiation and the saturated water vapor pressure difference as coordinate axes; The meteorological state space is divided into grids according to a preset interval, and the compensation index and the local vulnerability index are mapped as attribute features to the corresponding net radiation and saturated water vapor pressure difference coordinate points in the meteorological state space.

9. The method for compensating for aerodynamic drag and canopy drag and for quantitative evaluation of equivalence according to claim 8, characterized in that: Based on the compensation index and local vulnerability index within each meteorological state grid, regression analysis is performed on each meteorological state grid to determine the parameter constraint capability and model risk boundary corresponding to each meteorological state grid. Based on the parameter constraints and the model risk boundary, the risk level of the corresponding meteorological state is determined, and the compensation index, local vulnerability index, equivalent zone morphology, centerline slope and risk level corresponding to each meteorological state are output as quantitative assessment results.

Images