A satellite observation operator construction method and system for a micro-gridded earth system

Abstract

The application discloses a satellite observation operator construction method and system for a differentiable earth system, comprising the following steps: decomposing a satellite radiation transmission physical process into a plurality of differentiable calculation nodes, constructing a forward calculation graph, calculating an atmospheric state tensor, a prediction factor tensor, a total layer optical thickness tensor, a transmittance tensor, a top layer radiance tensor and an observed brightness temperature tensor through the differentiable calculation nodes, constructing a backward propagation calculation graph according to the forward calculation graph, and encapsulating the forward calculation graph, the backward propagation calculation graph and a standardized calling interface into a satellite observation operator. The method can not only realize end-to-end gradient propagation of an observation physical process, output observation and gradient information of the observation on an earth system state variable, but also has good interpretability and can be directly applied to a satellite observation operator construction system.

Description

Technical Field

[0001] This invention relates to the field of differentiable physical process representation technology, and in particular to a method and system for constructing satellite observation operators for differentiable Earth systems. Background Technology

[0002] Satellite microwave observations are a core observational constraint for numerical weather prediction and data assimilation systems. Operational systems rely on fast radiative transfer models (such as RTTOV and CRTM) to forward model state variables such as atmospheric temperature and humidity into simulated brightness temperatures, driving variational assimilation or ensemble filtering to correct the model initial fields. With the rise of cutting-edge research areas such as differentiable Earth systems and AI-NWP, observation operators must share a unified computational graph with dynamic processes and neural network components and achieve end-to-end gradient propagation, highlighting the increasing architectural bottlenecks of traditional observation operators.

[0003] Existing fast radiative transfer modes suffer from three fundamental drawbacks: First, gradient propagation relies on manually written tangent linear and adjoint modes, resulting in long development cycles, high maintenance costs, and the need to rewrite adjoint code when sensors are updated or physical parameters are adjusted, making it difficult to support the end-to-end gradient chain of observation-state-dynamics; second, the computational architecture is mainly based on CPU profile-by-profile serial processing, which cannot fully utilize the parallel computing power of heterogeneous hardware such as GPUs / TPUs, making it difficult to cope with the orders-of-magnitude increase in the number of channels and observation frequency of next-generation spaceborne instruments; finally, traditional modes exist as independent software packages and cannot be directly embedded into the unified computation graph of deep learning frameworks, leading to a break in the gradient chain between the observational physics process and the AI ​​model. Therefore, this invention proposes a method and system for constructing satellite observation operators for differentiable Earth systems. By decomposing the radiative transfer physical process into multiple differentiable computational nodes, a forward computation graph is constructed and an automatic differentiation mechanism is used to generate a backpropagation computation graph. Finally, the forward computation graph, the backpropagation computation graph, and the standardized calling interface are encapsulated into a satellite observation operator. This method achieves differentiable reconstruction while maintaining the parameterized system of radiative transfer physics. Through automatic differentiation, it realizes end-to-end gradient propagation and heterogeneous parallel computation of the observation physical process. It can serve as a core differentiable observation layer to solve the gradient chain breakage problem between existing operators and differentiable systems. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for constructing satellite observation operators for differentiable Earth systems.

[0005] To achieve the above objectives, the present invention is implemented according to the following technical solution: This invention includes the following steps: The physical process of satellite radiative transfer is decomposed into multiple differentiable computational nodes, and a forward computation graph is constructed. Acquire user layer observation data, organize user layer atmospheric profile data into user layer profile tensors in batch dimension, construct sparse interpolation weight matrix and vertically interpolate user layer profile tensors to coefficient layer grid to obtain coefficient layer atmospheric state tensors. The prediction factor tensor is constructed based on the atmospheric state tensor of the coefficient layer, and the total layer optical thickness is obtained by performing tensor shrinkage and summation operations in combination with the coefficient tensor. The cumulative optical thickness tensor is obtained by performing tensor quantization and cumulative summation on the total layer optical thickness tensor. The cumulative optical thickness tensor is then back-interpolated from the coefficient layer grid to the user layer grid using a conformal differentiable interpolation operator, and finally the user layer transmittance tensor is obtained by exponential transformation. The source function values ​​of each user layer are calculated and tensorized radiance integrals are performed in combination with the user layer transmittance tensor to obtain the top layer radiance tensor. An analytically differentiable inverse Planck transform is then performed on the top layer radiance tensor to obtain the observed brightness temperature tensor. Based on the forward computation graph, a backpropagation computation graph is constructed, and the forward computation graph, the backpropagation computation graph, and the standardized calling interface are encapsulated into a satellite observation operator; The user-layer observation data includes atmospheric profile data, surface boundary data, and observation geometric parameters; The coefficient layer state data includes atmospheric state tensor, reference profile, and observed zenith angle; The differentiable computation nodes include vertical interpolation nodes, prediction factor construction nodes, optical thickness calculation nodes, transmittance solution nodes, radiation integration nodes, and observation transformation nodes. The shape of the tensor calculation includes batch dimension, channel dimension, layer dimension, and prediction factor dimension; The forward computation graph is an ordered combination of differentiable computational nodes connected by tensor data streams, used to output the observed brightness temperature tensor. The backpropagation computation graph is obtained by backpropagation based on the differentiable computation nodes and is used to output the gradient information of the observed brightness temperature tensor on the user layer observation data. The standardized calling interface is used to provide a forward inference interface, a gradient query interface, and a diagnostic output interface, and to return intermediate physical diagnostic quantities as needed; the intermediate physical diagnostic quantities include layer optical thickness, layer transmittance, cumulative optical thickness, weighting function, contribution of each radiation component, prediction factor value, and layer source function value.

[0006] Furthermore, the method for obtaining the atmospheric state tensor of the coefficient layer includes: Organize the user-layer atmospheric profile data into user-layer profile tensors in the batch dimension; In the logarithmic coordinate system of user layer air pressure and coefficient layer air pressure, the interpolation interval is determined, and non-zero weights are assigned to coefficient layer positions located within the target interval, while zero weights are assigned to coefficient layer positions located outside the target interval, thus constructing a sparse interpolation weight matrix. The profile tensor is vertically interpolated to the coefficient layer grid in the batch dimension by matrix multiplication to obtain the atmospheric state tensor of the coefficient layer; the mapping process is fully parallel and differentiable in the batch dimension; the weights of the vertical interpolation are determined by the sparse interpolation weight matrix.

[0007] Furthermore, the method for obtaining the accumulated optical thickness includes: The slant path factor is calculated from the observed zenith angle, and the expression is: ; in For slant path factor, To observe the zenith angle; The prediction factor vector is constructed based on the atmospheric state tensor of the coefficient layer and the reference profile. The specific steps are as follows: The temperature deviation factor is determined by the deviation between the temperature of the coefficient layer and the reference temperature profile; the temperature ratio factor is determined by the ratio of the temperature of the coefficient layer to the reference temperature profile; and the water vapor ratio factor is determined by the ratio of the humidity of the coefficient layer to the reference humidity profile. The forecast factor vector is composed of the temperature deviation factor, the temperature ratio factor, the water vapor ratio factor, and the slant path factor. The forecast factor vector is calculated according to the type of absorption component; the types of absorption components include mixed gas absorption, water vapor line absorption, and water vapor continuum absorption. The forecast factor vectors are combined into a forecast factor tensor according to the channel, layer number, and absorption component type. Tensor contraction is then performed on the forecast factor tensor and coefficient tensor along the forecast factor dimension to obtain the component layer optical thickness. Finally, the components are summed along the absorption component dimension to obtain the total layer optical thickness. The expression is as follows: ; ; in For the optical thickness of the layer, For the absorption component type index, For channel index, For hierarchical indexing, The number of forecast factors in the forecast factor vector. For the forecast factor index, The absorption coefficient is... Forecasting factors, To accumulate optical thickness, The optical thickness of the same layer of mixed gas absorption components. The optical thickness of the same water vapor line absorption component layer. The optical thickness of the absorption component layer of the same water vapor continuum.

[0008] Furthermore, the method for obtaining the user layer transmittance tensor includes: First, tensor quantization cumulative summation is performed on the total layer optical thickness along the layer dimension to obtain the cumulative optical thickness of the coefficient layer; the cumulative optical thickness is monotonically non-decreasing along the layer dimension. The cumulative optical thickness tensor is back-interpolated from the coefficient layer grid to the user layer pressure grid using a conformal differentiable interpolation operator, and then the user layer transmittance tensor is obtained by exponential transformation. The conformal differentiable interpolation operator is the Rochon interpolation operator, specifically interpolating in various pressure ranges. The slope is determined using a cubic Hermite polynomial, and the node slope is determined according to the shape-preserving rule: when the slope of the current interval secant is opposite in sign to the slope of the previous interval secant or either is zero, the node slope is 0; otherwise, the node slope is the weighted harmonic average of the slopes of the current interval secant and the slope of the previous interval secant. The endpoint slope satisfies the one-sided constraint: the one-sided constraint is that the endpoint slope does not exceed three times the slope of the adjacent secant. The exponential transformation expression is: ; in For the user layer Node transmittance, For the user layer The cumulative optical thickness of the nodes is back-interpolated to the user layer pressure grid using a conformal differentiable interpolation operator.

[0009] Furthermore, the method for obtaining the observed brightness temperature tensor includes: Within the channel spectral response function range and the operating temperature range, the exact Planck function is least-squares fitted to a quadratic polynomial of temperature to obtain the band-corrected Planck function, expressed as: ; in To correct the Planck function for the band, For temperature, For channel index, , , For channel The band correction coefficients are predetermined in one step by weighted fitting of the spectral response function; The source function values ​​of each user layer are calculated using the band-corrected Planck function. The tensor radiance tensor of the top layer is obtained by performing tensor radiance integral along the layer dimension using the user layer transmittance tensor. Perform an analytically differentiable inverse Planck transform on the top-level radiance tensor to obtain the observed brightness temperature tensor; The steps of the analytically differentiable and inverse Planck transformation are as follows: First, solve the quadratic equation. Then take the real root that makes the brightness temperature fall within the physical temperature range. As the observed brightness temperature, the band correction coefficient is included. As it approaches 0, it degenerates into a linear inverse transform. The discriminant remains positive within the operating temperature range. The inverse transform is given in analytical closed form and is differentiable throughout.

[0010] Furthermore, the satellite observation operator outputs the observed brightness-temperature tensor and gradient information; the gradient information includes the scalar loss gradient, the Jacobian matrix, and higher-order gradients.

[0011] Secondly, a satellite observation operator construction system for differentiable Earth systems includes: Computational graph construction module: used to decompose the physical process of satellite radiative transfer into multiple differentiable computational nodes and construct a forward computational graph; Data input module: used to acquire user layer observation data and organize user layer atmospheric profile data into user layer profile tensors in batch dimension; The node execution module is used to construct a sparse interpolation weight matrix and vertically interpolate the user layer profile tensor to the coefficient layer grid to obtain the coefficient layer atmospheric state tensor; it is used to construct a forecast factor tensor based on the coefficient layer atmospheric state tensor, and perform tensor shrinkage and summation operations on the coefficient tensor to obtain the total layer optical thickness; it is used to perform tensor quantization cumulative summation on the total layer optical thickness tensor to obtain the cumulative optical thickness tensor, and use a conformal differentiable interpolation operator to interpolate the cumulative optical thickness tensor from the coefficient layer grid to the user layer grid, and then perform an exponential transformation to obtain the user layer transmittance tensor; it is used to calculate the source function values ​​of each user layer and perform tensor quantization radiance integral on the user layer transmittance tensor to obtain the top layer radiance tensor, and perform an analytically differentiable inverse Planck transform on the top layer radiance tensor to obtain the observed brightness temperature tensor; Automatic differentiation module: used to construct a backpropagation computation graph based on the forward computation graph; Operator encapsulation module: used to encapsulate the forward computation graph, the backpropagation computation graph, and the standardized calling interface into satellite observation operators.

[0012] The beneficial effects of this invention are: This invention provides a method and system for constructing satellite observation operators for differentiable Earth systems. Compared with existing technologies, this invention has the following technical advantages: This invention reconstructs the radiative transfer physical process into a differentiable observable operator in a unified computational graph, and generates a backpropagation computational graph based on an automatic differentiation mechanism, realizing end-to-end gradient propagation between observation, state, and dynamic processes, eliminating the need to manually write tangent linear modes and adjoint modes. This invention achieves differentiable reconfiguration while maintaining the physical parameterization system. Unlike purely data-driven neural network proxy models, it can output intermediate physical diagnostic quantities on demand through standardized calling interfaces, supporting error tracing and mechanism diagnosis. This invention replaces the sequential loop of profile-by-profile with batch tensor operations, reconstructing the physical process of observation into parallel tensor calculations along the batch dimension, channel dimension, layer dimension, and prediction factor dimension. The computational throughput scales linearly with the batch size and can be executed seamlessly on heterogeneous computing devices such as CPU / GPU / TPU. This invention supports multiple calling modes, including scalar loss gradient, Jacobian matrix, and higher-order gradient. The Jacobian matrix mode can automatically select forward or backward mode based on the relative size of the number of channels and the dimension of the state variables to optimize computational complexity. The higher-order gradient mode supports automatic solving of the second derivative and Hessian matrix. This invention provides forward inference, gradient query, and diagnostic output interfaces through standardized encapsulation of a unified differentiable observation operator. It can be directly embedded as a differentiable physical layer into a deep neural network, as a forward operator into a variational data assimilation system, or as a physical constraint term into a deep learning loss function, exhibiting good modularity and reusability. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating the steps of a satellite observation operator construction method for differentiable Earth systems according to the present invention. Figure 2 The brightness temperature statistics are for the differentiable model of the docking satellite observation operator and the traditional RTTOV model in this embodiment of the invention. Figure 3 This is a comparison of the temperature Jacobian profiles of the differentiable model of the docking satellite observation operator and the traditional RTTOV model in this embodiment of the invention. Figure 4 This is a comparison of the water vapor Jacobian profile between the differentiable model of the docking satellite observation operator and the traditional RTTOV model in this embodiment of the invention. Detailed Implementation

[0014] The present invention will be further described below through specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0015] This invention discloses a method and system for constructing satellite observation operators for differentiable Earth systems, comprising the following steps: like Figure 1 As shown, this embodiment includes the following steps: The physical process of satellite radiative transfer is decomposed into multiple differentiable computational nodes, and a forward computation graph is constructed. Acquire user layer observation data, organize user layer atmospheric profile data into user layer profile tensors in batch dimension, construct sparse interpolation weight matrix and vertically interpolate user layer profile tensors to coefficient layer grid to obtain coefficient layer atmospheric state tensors. The prediction factor tensor is constructed based on the atmospheric state tensor of the coefficient layer, and the total layer optical thickness is obtained by performing tensor shrinkage and summation operations in combination with the coefficient tensor. The cumulative optical thickness tensor is obtained by performing tensor quantization and cumulative summation on the total layer optical thickness tensor. The cumulative optical thickness tensor is then back-interpolated from the coefficient layer grid to the user layer grid using a conformal differentiable interpolation operator, and finally the user layer transmittance tensor is obtained by exponential transformation. The source function values ​​of each user layer are calculated and tensorized radiance integrals are performed in combination with the user layer transmittance tensor to obtain the top layer radiance tensor. An analytically differentiable inverse Planck transform is then performed on the top layer radiance tensor to obtain the observed brightness temperature tensor. Based on the forward computation graph, a backpropagation computation graph is constructed, and the forward computation graph, the backpropagation computation graph, and the standardized calling interface are encapsulated into a satellite observation operator; The user-layer observation data includes atmospheric profile data, surface boundary data, and observation geometric parameters; The coefficient layer state data includes atmospheric state tensor, reference profile, and observed zenith angle; The differentiable computation nodes include vertical interpolation nodes, prediction factor construction nodes, optical thickness calculation nodes, transmittance solution nodes, radiation integration nodes, and observation transformation nodes. The shape of the tensor calculation includes batch dimension, channel dimension, layer dimension, and prediction factor dimension; The forward computation graph is an ordered combination of differentiable computational nodes connected by tensor data streams, used to output the observed brightness temperature tensor. The backpropagation computation graph is obtained by backpropagation based on the differentiable computation nodes and is used to output the gradient information of the observed brightness temperature tensor on the user layer observation data. The standardized calling interface is used to provide a forward inference interface, a gradient query interface, and a diagnostic output interface, and to return intermediate physical diagnostic quantities as needed; the intermediate physical diagnostic quantities include layer optical thickness, layer transmittance, cumulative optical thickness, weighting function, contribution of each radiation component, prediction factor value, and layer source function value.

[0016] In this embodiment, the method for obtaining the atmospheric state tensor of the coefficient layer includes: Organize the user-layer atmospheric profile data into user-layer profile tensors in the batch dimension; In the logarithmic coordinate system of user layer air pressure and coefficient layer air pressure, the interpolation interval is determined, and non-zero weights are assigned to coefficient layer positions located within the target interval, while zero weights are assigned to coefficient layer positions located outside the target interval, thus constructing a sparse interpolation weight matrix. The profile tensor is vertically interpolated to the coefficient layer grid in the batch dimension by matrix multiplication to obtain the atmospheric state tensor of the coefficient layer; the mapping process is fully parallel and differentiable in the batch dimension; the weights of the vertical interpolation are determined by the sparse interpolation weight matrix. In practical assessments, the physical processes of satellite radiative transfer include microwave radiative transfer under the clear-sky pure absorption approximation, as well as radiative transfer processes with cloud and precipitation extensions. The observed physical processes are graphically represented while retaining the physical parameterization system and computational structure of the radiative transfer model. This physical parameterization system is based on a coefficient system provided by the fast radiative transfer model or a coefficient system generated by regression from a line-by-line integral model, including but not limited to layer optical thickness parameterization, absorption coefficient parameterization, continuum absorption parameterization, and source function parameterization. The computational operations in the computational graph nodes correspond to the specific mathematical expressions of the physical processes, ensuring the physical consistency and interpretability of the differentiable observable operators. The computational operations of the differentiable observable operators constitute the automatic differentiation mechanism of the tensor computational framework. The traditional radiative transfer model reconstructs serial cyclic computation based on the profile into parallel tensor computation along the batch dimension. The multidimensional scalar operations in the observation physics process are rewritten into tensor shrinking, tensor broadcasting, and tensor element-wise operations, forming a modern tensor computation framework. The batch parallel execution capability of the modern tensor computation framework is used to achieve synchronous computation of multiple samples, multiple channels, and multiple layers, so that the computational throughput of the observation operator expands linearly with the batch size. Tensor-based parallel computing supports execution on heterogeneous computing devices. The shapes of tensor-based parallel computing include profile tensors of [batch size, number of layers], coefficient tensors of [number of channels, number of layers, number of forecast factors], and output brightness-temperature tensors of [batch size, number of channels]. These heterogeneous computing devices include, but are not limited to, central processing units (CPUs), graphics processing units (GPUs), and tensor processors (TPUs). Through the device abstraction layer of the tensor computing framework, differentiable observable operators can be seamlessly migrated between different hardware platforms. The computation graph and its parameters are automatically deployed to the target device without manual intervention in data transfer and device scheduling.

[0017] In this embodiment, the method for obtaining the accumulated optical thickness includes: The slant path factor is calculated from the observed zenith angle, and the expression is: ; in For slant path factor, To observe the zenith angle; The prediction factor vector is constructed based on the atmospheric state tensor of the coefficient layer and the reference profile. The specific steps are as follows: The temperature deviation factor is determined by the deviation between the temperature of the coefficient layer and the reference temperature profile; the temperature ratio factor is determined by the ratio of the temperature of the coefficient layer to the reference temperature profile; and the water vapor ratio factor is determined by the ratio of the humidity of the coefficient layer to the reference humidity profile. The forecast factor vector is composed of the temperature deviation factor, the temperature ratio factor, the water vapor ratio factor, and the slant path factor. The forecast factor vector is calculated according to the type of absorption component; the types of absorption components include mixed gas absorption, water vapor line absorption, and water vapor continuum absorption. The forecast factor vectors are combined into a forecast factor tensor according to the channel, layer number, and absorption component type. Tensor contraction is then performed on the forecast factor tensor and coefficient tensor along the forecast factor dimension to obtain the component layer optical thickness. Finally, the components are summed along the absorption component dimension to obtain the total layer optical thickness. The expression is as follows: ; ; in For the optical thickness of the layer, For the absorption component type index, For channel index, For hierarchical indexing, The number of forecast factors in the forecast factor vector. For the forecast factor index, The absorption coefficient is... Forecasting factors, To accumulate optical thickness, The optical thickness of the same layer of mixed gas absorption components. The optical thickness of the same water vapor line absorption component layer. The optical thickness of the absorption component layer of the same water vapor continuum.

[0018] In this embodiment, the method for obtaining the user layer transmittance tensor includes: First, tensor quantization cumulative summation is performed on the total layer optical thickness along the layer dimension to obtain the cumulative optical thickness of the coefficient layer; the cumulative optical thickness is monotonically non-decreasing along the layer dimension. The cumulative optical thickness tensor is back-interpolated from the coefficient layer grid to the user layer pressure grid using a conformal differentiable interpolation operator, and then the user layer transmittance tensor is obtained by exponential transformation. The conformal differentiable interpolation operator is the Rochon interpolation operator, specifically interpolating in various pressure ranges. The slope is determined using a cubic Hermite polynomial, and the node slope is determined according to the shape-preserving rule: when the slope of the current interval secant is opposite in sign to the slope of the previous interval secant or either is zero, the node slope is 0; otherwise, the node slope is the weighted harmonic average of the slopes of the current interval secant and the slope of the previous interval secant. The endpoint slope satisfies the one-sided constraint: the one-sided constraint is that the endpoint slope does not exceed three times the slope of the adjacent secant. The exponential transformation expression is: ; in For the user layer Node transmittance, For the user layer The cumulative optical thickness of the nodes is back-interpolated to the user layer pressure grid using a conformal differentiable interpolation operator.

[0019] In this embodiment, the method for obtaining the observed brightness temperature tensor includes: Within the channel spectral response function range and the operating temperature range, the exact Planck function is least-squares fitted to a quadratic polynomial of temperature to obtain the band-corrected Planck function, expressed as: ; in To correct the Planck function for the band, For temperature, For channel index, , , For channel The band correction coefficients are predetermined in one step by weighted fitting of the spectral response function; The source function values ​​of each user layer are calculated using the band-corrected Planck function. The tensor radiance tensor of the top layer is obtained by performing tensor radiance integral along the layer dimension using the user layer transmittance tensor. Perform an analytically differentiable inverse Planck transform on the top-level radiance tensor to obtain the observed brightness temperature tensor; The steps of the analytically differentiable and inverse Planck transformation are as follows: First, solve the quadratic equation. Then take the real root that makes the brightness temperature fall within the physical temperature range. As the observed brightness temperature, the band correction coefficient is included. As it approaches 0, it degenerates into a linear inverse transform. The discriminant remains positive within the operating temperature range. The inverse transform is given in analytical closed form and is differentiable throughout.

[0020] In this embodiment, the satellite observation operator outputs the observed brightness-temperature tensor and gradient information; the gradient information includes the scalar loss gradient, the Jacobian matrix, and higher-order gradients. In actual assessment, relying on the automatic differentiation mechanism of the tensor computation framework, the backpropagation computation graph is automatically constructed based on the forward computation graph to realize the automatic solution of the gradient of the observed brightness temperature on the atmospheric state variables. The automatic differentiation mechanism replaces the manually written tangent linear mode and adjoint mode in the traditional model, and transforms the gradient solution process from program development to automatic derivation of the computation graph. Scalar loss gradient is suitable for scenarios where a scalar loss function is constructed from observations. It obtains the gradient tensor of the loss function with respect to all state variables through a single backpropagation. Jacobian matrix is ​​suitable for scenarios that require the complete Jacobian matrix of the observation operator. It calculates the Jacobian matrix of the observation with respect to all state variables in batches through the vector-Jacobi product interface, and automatically selects forward or backward mode based on the relative size of the number of channels and the dimension of the state variables to optimize computational complexity. Higher-order gradient supports automatic solution of second derivatives and Hessian matrix, providing higher-order gradient information for second-order optimization algorithms, observation sensitivity analysis, and uncertainty quantification. The encapsulated satellite observation operator is linked to differentiable Earth system models, Earth fundamental models, end-to-end AI numerical weather prediction systems, and physically constrained deep learning frameworks through a standardized integration interface, achieving the following functions: It is embedded as a physical layer in a deep neural network to participate in end-to-end training; it is embedded as a forward modeling operator in a variational data assimilation system to participate in cost function minimization iteration; it is embedded as a physical constraint term in a deep learning loss function to provide physical consistency constraints; and it serves as a physical diagnostic tool to output intermediate physical diagnostic quantities such as layer optical thickness, layer transmittance, cumulative optical thickness, weight function, contribution of each radiation component, forecast factor value, and layer source function value. Among them, differentiable Earth system models include differentiable atmospheric circulation models, differentiable land surface models, and differentiable ocean models; Earth fundamental models include Earth system models based on the Transformer architecture, Earth system representation models based on graph neural networks, and Earth system evolution models based on neural operators; end-to-end AI numerical weather prediction systems include differentiable forecast-observation closed-loop systems, end-to-end variational assimilation systems, and neural network weather prediction systems; and physically constrained deep learning frameworks include physical information neural networks, physically constrained generative models, and physically guided reinforcement learning frameworks. The satellite observation operator achieves differentiable reconstruction while maintaining the physical structure and parameterization system of radiative transfer. Unlike the pure data-driven proxy model that replaces the physical process with a neural network, the differentiable observation operator has physical consistency, cross-scenario generalization ability, interpretability of internal calculation process, scalability of new sensors, and numerical stability in long-term operation. Satellite observation operators are not only applicable to microwave observations, but also to infrared (replacing infrared band absorption coefficient parameterization and source function parameterization), hyperspectral (expanding channel dimension to high spectral resolution and adapting to corresponding spectral response function processing), radar (introducing scattering parameterization and active radiative transfer calculation), GNSS radio occultation (adapting to bending angle or refractive index observation operators), lidar (introducing backscattering coefficient and extinction coefficient parameterization), and multi-source heterogeneous observations (combining differentiable operators of multiple observation types under a unified computation graph framework to construct a multi-source joint differentiable observation layer), supporting the expansion of a unified differentiable observation layer for different observation types.

[0021] In addition, another set of data examples is provided to compare the effects of the satellite observation operator construction method for differentiable Earth systems of the present invention. The specific implementation of the present invention is illustrated by taking the MWHS-2 (15-channel) microwave humidity meter carried by the Fengyun-3 series satellite as an example, but the present invention is not limited thereto. The method is also applicable to other spaceborne passive remote sensing instruments in the microwave and submillimeter wave bands. By sequentially performing the above steps, a satellite observation operator for differentiable Earth systems can be constructed. To verify the correctness of the forward calculation of this method, under the condition of 15 channels and multi-layer standard profile samples in MWHS-2, the forward brightness temperature of this method is compared with the results of the traditional fast radiative transfer mode RTTOV under the same coefficient system. Figure 2 The statistical results of brightness temperature deviation and standard deviation of the differentiable model of this invention and the traditional RTTOV model on all 15 channels are presented, among which... Figure 2 (a) Brightness temperature deviation, Figure 2 (b) represents the standard deviation of the brightness temperature. The results show that the forward brightness temperature of this method is highly consistent with that of the traditional mode across all channels, with an average channel deviation of less than 0.02K and a standard deviation of less than 0.015K, which proves the accuracy of this method in reproducing the physical behavior of the fast radiative transfer mode. To verify the accuracy of the automatic differential gradient solution proposed in this method, the Jacobian of the brightness temperature versus temperature profile obtained by automatic differentiation is compared with the output results of the traditional Jacobian model of RTTOV. Figure 3 The results of temperature Jacobian profile comparison of all 15 channels of MWHS-2 are shown. The solid line in the figure is the output of the traditional Jacobian mode of RTTOV, the dashed line is the output of the automatic differentiation of the differentiable model of this invention, and the gray area is the difference between the two. The results show that the two are highly consistent in terms of magnitude and vertical structure, with a maximum difference of no more than 0.5%, proving that the automatic differentiation mechanism of this method can accurately replace the traditional manually written adjoint code and provide reliable gradient information for temperature profile sensitivity analysis. The Jacobian of the brightness temperature versus water vapor profile obtained by this method based on automatic differentiation is compared with the output results of the traditional RTTOV Jacobian model. Figure 4 The results of the comparison of water vapor Jacobian profiles in channels 11 to 15 of MWHS-2 are presented. Channels 11 to 15 of MWHS-2 are located near the 183.31 GHz water vapor absorption line and are highly sensitive to the vertical distribution of atmospheric humidity. They are key channels for verifying the accuracy of automatic differential humidity gradient solution. The comparative results show that the water vapor Jacobian obtained by the automatic differentiation of the differentiable model of this invention is highly consistent with the output of the traditional RTTOV Jacobian model in terms of both vertical structure and magnitude, verifying the accuracy of the automatic differentiation mechanism of this invention in calculating humidity sensitivity. This result is of great significance for data assimilation of water vapor variables in differentiable Earth systems, physical constraint deep learning, and end-to-end intelligent assimilation applications.

[0022] Secondly, a satellite observation operator construction system for differentiable Earth systems includes: Computational graph construction module: used to decompose the physical process of satellite radiative transfer into multiple differentiable computational nodes and construct a forward computational graph; Data input module: used to acquire user layer observation data and organize user layer atmospheric profile data into user layer profile tensors in batch dimension; The node execution module is used to construct a sparse interpolation weight matrix and vertically interpolate the user layer profile tensor to the coefficient layer grid to obtain the coefficient layer atmospheric state tensor; it is used to construct a forecast factor tensor based on the coefficient layer atmospheric state tensor, and perform tensor shrinkage and summation operations on the coefficient tensor to obtain the total layer optical thickness; it is used to perform tensor quantization cumulative summation on the total layer optical thickness tensor to obtain the cumulative optical thickness tensor, and use a conformal differentiable interpolation operator to interpolate the cumulative optical thickness tensor from the coefficient layer grid to the user layer grid, and then perform an exponential transformation to obtain the user layer transmittance tensor; it is used to calculate the source function values ​​of each user layer and perform tensor quantization radiance integral on the user layer transmittance tensor to obtain the top layer radiance tensor, and perform an analytically differentiable inverse Planck transform on the top layer radiance tensor to obtain the observed brightness temperature tensor; Automatic differentiation module: used to construct a backpropagation computation graph based on the forward computation graph; Operator encapsulation module: used to encapsulate the forward computation graph, the backpropagation computation graph, and the standardized calling interface into satellite observation operators.

[0023] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for constructing satellite observation operators for differentiable Earth systems, characterized in that, Includes the following steps: S1. Decompose the physical process of satellite radiative transfer into multiple differentiable computational nodes and construct a forward computation graph; S2. Obtain user layer observation data, organize user layer atmospheric profile data into user layer profile tensors in batch dimension, construct sparse interpolation weight matrix and vertically interpolate user layer profile tensors to coefficient layer grid to obtain coefficient layer atmospheric state tensors. S3. Construct a forecast factor tensor based on the atmospheric state tensor of the coefficient layer, and perform tensor shrinkage and summation operations in combination with the coefficient tensor to obtain the total layer optical thickness. S4. Perform tensor quantization and cumulative summation on the total layer optical thickness tensor to obtain the cumulative optical thickness tensor. Use a conformal differentiable interpolation operator to interpolate the cumulative optical thickness tensor from the coefficient layer grid to the user layer grid, and then perform an exponential transformation to obtain the user layer transmittance tensor. S5. Calculate the source function values ​​of each user layer and perform tensorized radiative integral with the user layer transmittance tensor to obtain the top layer radiance tensor. Perform analytical differentiable inverse Planck transform on the top layer radiance tensor to obtain the observed brightness temperature tensor. S6. Construct a backpropagation computation graph based on the forward computation graph, and encapsulate the forward computation graph, the backpropagation computation graph, and the standardized calling interface into a satellite observation operator; The user-layer observation data includes atmospheric profile data, surface boundary data, and observation geometric parameters; The differentiable computation nodes include vertical interpolation nodes, prediction factor construction nodes, optical thickness calculation nodes, transmittance solution nodes, radiation integration nodes, and observation transformation nodes. The shape of the tensor calculation includes batch dimension, channel dimension, layer dimension, and prediction factor dimension; The forward computation graph is an ordered combination of differentiable computational nodes connected by tensor data streams, used to output the observed brightness temperature tensor. The backpropagation computation graph is obtained by backpropagation based on the differentiable computation nodes and is used to output the gradient information of the observed brightness temperature tensor on the user layer observation data. The standardized calling interface is used to provide a forward inference interface, a gradient query interface, and a diagnostic output interface, and to return intermediate physical diagnostic quantities as needed; the intermediate physical diagnostic quantities include layer optical thickness, layer transmittance, cumulative optical thickness, weighting function, contribution of each radiation component, prediction factor value, and layer source function value.

2. The method for constructing satellite observation operators for differentiable Earth systems according to claim 1, characterized in that, The method for obtaining the atmospheric state tensor of the coefficient layer includes: Organize the user-layer atmospheric profile data into user-layer profile tensors in the batch dimension; In the logarithmic coordinate system of user layer air pressure and coefficient layer air pressure, the interpolation interval is determined, and non-zero weights are assigned to coefficient layer positions located within the target interval, while zero weights are assigned to coefficient layer positions located outside the target interval, thus constructing a sparse interpolation weight matrix. The profile tensor is vertically interpolated to the coefficient layer grid in the batch dimension by matrix multiplication to obtain the atmospheric state tensor of the coefficient layer; the mapping process is fully parallel and differentiable in the batch dimension; the weights of the vertical interpolation are determined by the sparse interpolation weight matrix.

3. The method for constructing satellite observation operators for differentiable Earth systems according to claim 1, characterized in that, The method for obtaining the cumulative optical thickness includes: The slant path factor is calculated from the observed zenith angle, and the expression is: ； in For slant path factor, To observe the zenith angle; The prediction factor vector is constructed based on the atmospheric state tensor of the coefficient layer and the reference profile. The specific steps are as follows: The temperature deviation factor is determined by the deviation between the temperature of the coefficient layer and the reference temperature profile; the temperature ratio factor is determined by the ratio of the temperature of the coefficient layer to the reference temperature profile; and the water vapor ratio factor is determined by the ratio of the humidity of the coefficient layer to the reference humidity profile. The forecast factor vector is composed of the temperature deviation factor, the temperature ratio factor, the water vapor ratio factor, and the slant path factor. The forecast factor vector is calculated according to the type of absorption component; the types of absorption components include mixed gas absorption, water vapor line absorption, and water vapor continuum absorption. The forecast factor vectors are combined into a forecast factor tensor according to the channel, layer number, and absorption component type. Tensor contraction is then performed on the forecast factor tensor and coefficient tensor along the forecast factor dimension to obtain the component layer optical thickness. Finally, the components are summed along the absorption component dimension to obtain the total layer optical thickness. The expression is as follows: ； ； in For the optical thickness of the layer, For the absorption component type index, For channel indexing, For hierarchical indexing, The number of forecast factors in the forecast factor vector. For the forecast factor index, The absorption coefficient is... Forecasting factors, To accumulate optical thickness, The optical thickness of the same layer of mixed gas absorption components. The optical thickness of the same water vapor line absorption component layer. The optical thickness of the absorption component layer of the same water vapor continuum.

4. The method for constructing satellite observation operators for differentiable Earth systems according to claim 1, characterized in that, The method for obtaining the user layer transmittance tensor includes: First, tensor quantization cumulative summation is performed on the total layer optical thickness along the layer dimension to obtain the cumulative optical thickness of the coefficient layer; the cumulative optical thickness is monotonically non-decreasing along the layer dimension. The cumulative optical thickness tensor is back-interpolated from the coefficient layer grid to the user layer pressure grid using a conformal differentiable interpolation operator, and then the user layer transmittance tensor is obtained by exponential transformation. The conformal differentiable interpolation operator is the Rochon interpolation operator, specifically interpolating in various pressure ranges. The slope is determined using a cubic Hermite polynomial, and the node slope is determined according to the shape-preserving rule: when the slope of the current interval secant is opposite in sign to the slope of the previous interval secant or either is zero, the node slope is 0; otherwise, the node slope is the weighted harmonic average of the slopes of the current interval secant and the slope of the previous interval secant. The endpoint slope satisfies the one-sided constraint: the one-sided constraint is that the endpoint slope does not exceed three times the slope of the adjacent secant. The exponential transformation expression is: ； in For the user layer Node transmittance, For the user layer The cumulative optical thickness of the nodes is back-interpolated to the user layer pressure grid using a conformal differentiable interpolation operator.

5. The method for constructing satellite observation operators for differentiable Earth systems according to claim 1, characterized in that, The method for obtaining the observed brightness temperature tensor includes: Within the channel spectral response function range and the operating temperature range, the exact Planck function is least-squares fitted to a quadratic polynomial of temperature to obtain the band-corrected Planck function, expressed as: ； in To correct the Planck function for the band, For temperature, For channel indexing, , , For channel The band correction coefficients are predetermined in one step by weighted fitting of the spectral response function; The source function values ​​of each user layer are calculated using the band-corrected Planck function. The tensor radiance tensor of the top layer is obtained by performing tensor radiance integral along the layer dimension using the user layer transmittance tensor. Perform an analytically differentiable inverse Planck transform on the top-level radiance tensor to obtain the observed brightness temperature tensor; The steps of the analytically differentiable and inverse Planck transformation are as follows: First, solve the quadratic equation. Then take the real root that makes the brightness temperature fall within the physical temperature range. As the observed brightness temperature, the band correction coefficient is included. As it approaches 0, it degenerates into a linear inverse transform. The discriminant remains positive within the operating temperature range. The inverse transform is given in analytical closed form and is differentiable throughout.

6. The method for constructing satellite observation operators for differentiable Earth systems according to claim 1, characterized in that: The satellite observation operator outputs the observed brightness-temperature tensor and gradient information; the gradient information includes the scalar loss gradient, the Jacobian matrix, and higher-order gradients.

7. A satellite observation operator construction system for differentiable Earth systems, used to perform the method according to any one of claims 1-6, characterized in that, include: Computational graph construction module: used to decompose the physical process of satellite radiative transfer into multiple differentiable computational nodes and construct a forward computational graph; Data input module: used to acquire user layer observation data and organize user layer atmospheric profile data into user layer profile tensors in batch dimension; Node execution module: Used to construct a sparse interpolation weight matrix and vertically interpolate the user layer profile tensor to the coefficient layer grid to obtain the coefficient layer atmospheric state tensor; used to construct a forecast factor tensor based on the coefficient layer atmospheric state tensor, and perform tensor shrinkage and summation operations on the coefficient tensor to obtain the total layer optical thickness; used to perform tensor quantization cumulative summation on the total layer optical thickness tensor to obtain the cumulative optical thickness tensor, and use a conformal differentiable interpolation operator to interpolate the cumulative optical thickness tensor from the coefficient layer grid to the user layer grid, and then perform an exponential transformation to obtain the user layer transmittance tensor; The source function values ​​of each user layer are used to calculate the source function values ​​and perform tensorized radiance integrals in combination with the user layer transmittance tensor to obtain the top layer radiance tensor. An analytically differentiable inverse Planck transform is then performed on the top layer radiance tensor to obtain the observed brightness temperature tensor. Automatic differentiation module: used to construct a backpropagation computation graph based on the forward computation graph; Operator encapsulation module: used to encapsulate the forward computation graph, the backpropagation computation graph, and the standardized calling interface into satellite observation operators.

Images