A weather forecasting method based on an AI weather forecasting model and a WRF model

Abstract

This invention relates to the field of AI weather forecasting technology and discloses a weather forecasting method based on an AI weather forecasting model and a WRF model. The method includes: inputting ERA5 global atmospheric reanalysis gridded data into the AI ​​weather forecasting model to obtain an assimilated background field; constructing an effective multi-source observation vector based on observation data from ground automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data; solving for the optimal atmospheric analysis field; extracting atmospheric state variables from the optimal atmospheric analysis field and writing them into the WRF model's side boundary driving file, driving the WRF model to complete the integration operation, and obtaining a gridded forecast field of meteorological elements. This method enables the WRF boundary condition update frequency to be higher than the traditional 6-hour static scheme, effectively characterizing the rapid evolution of the atmospheric state in the boundary region at the sub-hourly scale, thereby improving the WRF model's forecasting capability for high-impact weather events such as severe convection.

Description

Technical Field

[0001] This invention relates to the field of AI weather forecasting technology, and in particular to a weather forecasting method based on an AI weather forecasting model and a WRF model. Background Technology

[0002] The WRF model in meteorological forecasting is currently the most widely used regional numerical weather prediction system. Its forecast accuracy is highly dependent on the accuracy of the lateral boundary conditions. However, existing technologies usually use the 6-hour interval static field output by the GFS or ECMWF global model to directly drive the WRF lateral boundary after interpolation. This causes the global model error to continue to propagate into the computational domain with the WRF integral, and no real-time observation data can be obtained for constraint correction throughout the entire forecast period.

[0003] In recent years, AI weather forecasting models, represented by Pangu-Weather, have demonstrated their advantages in global atmospheric state prediction due to their high timeliness. However, there is a significant systematic bias between AI model predictions and actual observation data. If AI model predictions are directly used to drive WRF side boundary conditions without observational constraint correction, the impact of background field errors on WRF forecast accuracy cannot be eliminated. Currently, there is no established technical solution for dynamically driving WRF model side boundary conditions after assimilation and correction of AI weather forecasting model predictions using multi-source observation data. Quantitative fusion methods for multi-source heterogeneous observation data and AI-predicted background fields, hourly dynamic boundary condition generation mechanisms, and boundary relaxation forcing strategies all lack complete technical implementation paths, hindering further improvements in the accuracy of regional refined weather forecasts. Summary of the Invention

[0004] This invention provides a weather forecasting method based on an AI weather forecasting model and a WRF model. This invention enables the WRF boundary condition update frequency to be higher than that of the traditional 6-hour static scheme, and effectively characterizes the rapid evolution of the atmospheric state in the boundary area at the sub-hour scale, thereby improving the WRF model's ability to forecast high-impact weather events such as severe convection.

[0005] In a first aspect, the present invention provides a weather forecasting method based on an AI weather forecasting model and a WRF model, wherein the weather forecasting method based on the AI ​​weather forecasting model and the WRF model includes: ERA5 global atmospheric reanalysis gridded data was input into an AI weather forecasting model for autoregressive iterative inference and bias correction to obtain the assimilated background field. An effective multi-source observation vector is constructed based on observation data from ground automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data. Using the assimilated background field as a priori constraint and the effective multi-source observation vector as an observation constraint, the optimal atmospheric analysis field is solved. The atmospheric state variables of each grid point in the WRF model boundary relaxation zone are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary driving file to drive the WRF model to complete the integration operation and obtain the meteorological element grid point forecast field.

[0006] In conjunction with the first aspect, in the first implementation of the first aspect of the present invention, ERA5 global atmospheric reanalysis gridded data is input into an AI weather forecasting model for autoregressive iterative inference and bias correction to obtain an assimilated background field, including: Input the ERA5 global atmospheric reanalysis gridded data into the AI ​​weather forecasting model to perform autoregressive iterative inference and obtain the global atmospheric element gridded prediction field. The mean of the systematic bias of the global atmospheric element grid prediction field is statistically analyzed, and bias correction is performed on the global atmospheric element grid prediction field based on the mean of the systematic bias to obtain the assimilated background field.

[0007] In conjunction with the first aspect, in a second implementation of the first aspect of the present invention, ERA5 global atmospheric reanalysis gridded data is input into an AI weather forecasting model to perform autoregressive iterative inference, resulting in a global atmospheric element gridded prediction field, including: The input state tensor of the AI ​​weather forecasting model is constructed based on the upper atmospheric element field and the surface atmospheric element field in the ERA5 global atmospheric reanalysis grid data. The input state tensor is input into the AI ​​weather forecasting model to perform autoregressive iterative inference, and the atmospheric state tensor output at each step of the inference is used as the input state tensor for the next step of the inference. The process is iterated until the forecast lead time is reached to obtain the global atmospheric element grid prediction field.

[0008] In conjunction with the first aspect, in the third implementation of the first aspect of the present invention, an effective multi-source observation vector is constructed based on ground automatic weather station observation data, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data, including: Interpolation processing is performed on the first observation data of the ground automatic weather station to obtain the ground automatic weather station observation data; The second radiosonde observation data is standardized to obtain radiosonde observation data. CO2 slice cloud detection and data removal are performed on the raw radiance data of geostationary meteorological satellites to obtain geostationary meteorological satellite radiance data; The signal-to-noise ratio threshold test and data removal were performed on the raw VAD wind profile data of the Doppler radar to obtain the Doppler radar VAD wind profile data. An effective multi-source observation vector is constructed based on the observation data from the ground automatic weather station, the radiosonde observation data, the radiance data from the geostationary meteorological satellite, and the Doppler radar VAD wind profile data.

[0009] In conjunction with the first aspect, in the fourth implementation of the first aspect of the present invention, after constructing the effective multi-source observation vector, it further includes: The first instrument error variance of the ground automatic weather station, the second instrument error variance of the radiosonde, the third instrument error variance of the geostationary meteorological satellite, and the fourth instrument error variance of various instruments of the Doppler radar were obtained respectively. Based on the first instrument error variance, the second instrument error variance, the third instrument error variance, and the fourth instrument error variance, the diagonal elements of each observation element in the effective multi-source observation vector are assigned values ​​respectively to construct the observation error covariance matrix.

[0010] In conjunction with the first aspect, in the fifth implementation of the first aspect of the present invention, the optimal atmospheric analysis field is solved using the assimilated background field as a priori constraint and the effective multi-source observation vector as an observation constraint, including: A background constraint term is constructed based on the assimilated background field, and an observation constraint term is constructed based on the effective multi-source observation vector and the observation error covariance matrix. The background constraint term and the observation constraint term are summed to obtain the first cost function. The atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable to obtain the second cost function. The optimal atmospheric analysis field is obtained by solving the second cost function through iterative minimization.

[0011] In conjunction with the first aspect, in the sixth implementation of the first aspect of the present invention, the atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable to obtain a second cost function, and the optimal atmospheric analysis field is obtained by iteratively minimizing the second cost function, including: Based on the background error covariance matrix, the atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable, resulting in a second cost function with the incremental control variable as the independent variable; wherein, the second cost function is: ; The second cost function value. Incremental control variables, Background constraint term in the second cost function. Assimilation of the background field, Background error covariance matrix The square root operator, Nonlinear observation operator : Effective multi-source observation vector Observation error covariance matrix The inverse matrix; The second cost function is iteratively minimized with respect to the incremental control variables to obtain the optimal incremental control variables. The optimal incremental control variables are then inversely transformed using the background error covariance matrix and superimposed onto the assimilation background field to obtain the optimal atmospheric analysis field.

[0012] In conjunction with the first aspect, in the seventh implementation of the first aspect of the present invention, atmospheric state variables of each grid point in the WRF model boundary relaxation region are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary drive file, driving the WRF model to complete the integration operation, and obtaining the meteorological element grid point forecast field, including: The atmospheric state variables of each grid point in the boundary relaxation region of the WRF model are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary drive file; The WRF model is driven by the side boundary driving file of the WRF model to complete the integral operation and output the meteorological element grid forecast field of surface temperature, surface wind speed, surface air pressure, precipitation, cloud cover and convective available potential energy.

[0013] In conjunction with the first aspect, in the eighth implementation of the first aspect of the present invention, extracting the atmospheric state variables of each grid point in the WRF model boundary relaxation region from the optimal atmospheric analysis field and writing them into the WRF model side boundary drive file includes: In the optimal atmospheric analysis field, global latitude and longitude grid data are interpolated to the WRF model, and logarithmic pressure linear interpolation is performed with the pressure value corresponding to each η-mixed σ-pressure coordinate layer as the interpolation node, and then transformed to the η-mixed σ-pressure coordinate layer of the WRF model. The atmospheric state variables and corresponding trend terms of each grid point in the WRF model boundary relaxation region are extracted by the η-mixed σ-pressure coordinate layer of the WRF model, and the atmospheric state variables and the trend terms are written into the WRF model side boundary driving file.

[0014] In conjunction with the first aspect, in the ninth implementation of the first aspect of the present invention, the WRF model is driven to complete the integral operation based on the WRF model side boundary driving file, and outputs a gridded forecast field of meteorological elements including surface air temperature, surface wind speed, surface air pressure, precipitation, cloud cover, and convective available potential energy, including: Based on the WRF model side boundary driving file, the spatial attenuation coefficient of each grid point in the WRF model boundary relaxation region is calculated, and the atmospheric dynamic equations of the WRF model are constructed according to the spatial attenuation coefficient and the atmospheric state variables. The atmospheric dynamic equations of the WRF model are driven by a third-order Runge-Kutta time integration scheme to perform numerical integration on the ARW dynamic core of the WRF model, outputting a gridded forecast field of meteorological elements such as surface temperature, surface wind speed, surface air pressure, precipitation, cloud cover, and convective available potential energy.

[0015] The technical solution provided by this invention combines AI weather forecast model predictions with multi-source observation data and dynamically drives the WRF model's side boundary conditions. It replaces the traditional 6-hour interval global model static output field with the hourly global atmospheric state variable gridded prediction field output by the Pangu-Weather model. Furthermore, it performs grid-by-grid bias correction on the prediction field by statistically analyzing the systematic bias mean of historical forecast-analysis field difference samples grouped by variable, isobaric layer, and season, eliminating the systematic bias of the AI ​​model predictions and making the assimilated background field closer to the real atmospheric state. By performing targeted quality control on four types of heterogeneous observation data—from automatic weather stations, radiosondes, geostationary meteorological satellites, and Doppler radar VAD wind profiles—and constructing an iterative minimization solution for the cost function, the optimal atmospheric analysis gridded field constrained by real-time multi-source observations is used as the driving source for the WRF side boundary conditions. This achieves quantitative correction of the AI ​​predicted background field by observational information, effectively suppressing the propagation and accumulation of background field errors during the WRF integration process. After coordinate transformation, the optimal atmospheric analysis grid field is written into the WRF model's side boundary condition file at hourly update intervals. Combined with the Logistic function-type spatial decay relaxation forcing mechanism, the WRF boundary condition update frequency is higher than that of the traditional 6-hour static scheme. The rapid evolution of the atmospheric state in the boundary region at the sub-hourly scale is effectively characterized, thereby improving the WRF model's forecasting capability for high-impact weather events such as severe convection.

[0016] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention are realized and obtained in accordance with the structures particularly pointed out in the description, claims and drawings.

[0017] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of an embodiment of the weather forecasting method based on the AI ​​weather forecasting model and the WRF model in this invention. Figure 2 This is a schematic diagram of an embodiment of the AI ​​weather forecasting model performing autoregressive iterative inference and bias correction in this invention. Figure 3This is a schematic diagram of an embodiment of the present invention that drives a WRF model to complete the integral operation; Figure 4 This is a schematic diagram of the standard isobaric layer vertical structure in an embodiment of the present invention; Figure 5 This is a schematic diagram of the logarithmic pressure linear interpolation process in an embodiment of the present invention; Figure 6 This is a schematic diagram of the global latitude and longitude grid of the optimal atmospheric analysis field in an embodiment of the present invention; Figure 7 This is a schematic diagram illustrating the process of constructing the atmospheric dynamic equations of the WRF model in an embodiment of the present invention. Detailed Implementation

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

[0020] The terms "comprising" and "having," and any variations thereof, used in the embodiments of this invention are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the steps or units listed, but may optionally include other steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices.

[0021] To facilitate understanding of this embodiment, a weather forecasting method based on an AI weather forecasting model and a WRF model, as disclosed in this embodiment of the invention, will first be described in detail. For example... Figure 1 As shown, this method includes the following steps: 101. Input the ERA5 global atmospheric reanalysis gridded data into the AI ​​weather forecasting model for autoregressive iterative inference and bias correction to obtain the assimilated background field; 102. Construct an effective multi-source observation vector based on ground automatic weather station observation data, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data; 103. Solve for the optimal atmospheric analysis field by using the assimilated background field as a priori constraint and the effective multi-source observation vector as an observation constraint; 104. Extract the atmospheric state variables of each grid point in the relaxation zone of the WRF model boundary from the optimal atmospheric analysis field and write them into the WRF model side boundary driving file. Drive the WRF model to complete the integration operation and obtain the meteorological element grid point forecast field.

[0022] like Figure 2 In one specific embodiment, the process of performing step 101 may specifically include the following steps: 201. Input the ERA5 global atmospheric reanalysis grid data into the AI ​​weather forecast model and perform autoregressive iterative inference to obtain the global atmospheric element grid prediction field; 202. Calculate the mean systematic bias of the global atmospheric element grid prediction field, and perform bias correction on the global atmospheric element grid prediction field based on the mean systematic bias to obtain the assimilated background field.

[0023] Specifically, at the forecast start time, ERA5 global atmospheric reanalysis grid data is acquired, and the upper-air field and surface field are uniformly organized according to the input structure of the AI ​​weather forecast model. The upper-air field uses 13 standard isobaric layers, namely 1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100 and 50 hPa. Five types of atmospheric state variables, namely air temperature, zonal wind speed, meridional wind speed, specific humidity and geopotential height, are input into each isobaric layer. Four types of near-surface variables, namely 10-meter zonal wind speed, 10-meter meridional wind speed, 2-meter air temperature and sea level pressure, are input into the surface field. The overall grid resolution is 721×1440, forming the atmospheric state tensor at the current time t and sending it into the AI ​​weather forecast model for forward inference. After completing the calculation of the current input, the AI ​​weather forecast model outputs the atmospheric state tensor at the predicted time t+Δt. The time step Δt is fixed at 1 hour. The current output is then used as the input for the next time step, and the model continues to perform hourly rolling autoregressive iterative inference, continuously generating a global atmospheric element grid prediction field set for the next 0–72 hours. The prediction field covers 13 types of meteorological elements, including 2-meter air temperature, 10-meter zonal wind speed, 10-meter meridional wind speed, surface air pressure, sea level air pressure, air temperature at each isobaric layer, zonal wind speed at each isobaric layer, meridional wind speed at each isobaric layer, specific humidity at each isobaric layer, geopotential height at each isobaric layer, vertical velocity at each isobaric layer, surface sensible heat flux, and surface latent heat flux. During the model training phase, hourly global grid fields from 1979 to 2018 (40 years in total) are used as training samples, and (X... t ,X t+Δt Approximately 350,400 training sample pairs are organized in a paired manner. The training loss function uses dimension-weighted mean square error, and its expression is as follows: ; in, This represents the training loss value. Indicates the index of meteorological variables. Represents a global grid index. Indicates the first The latitude of each grid point This represents the values ​​of the actual analytical field at the corresponding variables and grid points. This represents the model's predicted value. This represents the latitude weighting factor, used to compensate for the uneven sample density caused by grid area shrinkage in high-latitude regions. During training, the optimizer uses AdamW, with an initial learning rate set to... The weight decay coefficient was set to 0.1, the training batch size to 8, and the training epochs to 100. The hardware consisted of 8 A100 GPUs, each with 80GB of VRAM, and the total training time was approximately 14 days. After model construction was completed, during the operational phase, a systematic bias estimation was performed on the global atmospheric element gridded prediction field output by the AI ​​weather forecast model to mitigate the impact of long-term statistical errors on the assimilation of the background field quality. Therefore, the data from the past 30 days of simultaneous daily data collection was retrieved. The grid-by-grid difference samples between the AI ​​weather forecast model predictions and the corresponding ERA5 real-time analysis fields were analyzed and grouped by variable, isobaric layer, and season. The seasons were divided into spring (March to May), summer (June to August), autumn (September to November), and winter (December to February). The mean of the systematic bias was calculated using the following formula: ; in, Representing variables isobaric layer ,season mean of systematic deviation under the condition, This indicates the number of days in the statistical sample, with a value of 30. Indicates the first The predicted values ​​of the AI ​​weather forecast model for each sample day and time. Indicates the first The true values ​​of the ERA5 analysis field for each sample day and time period are given. To characterize the dispersion of the background error, the variance of the background error is calculated simultaneously, and its expression is as follows: ; in, This represents the background error variance corresponding to the variable and the isobaric layer. A 30-day rolling statistical window is used because a window length of less than 7 days results in insufficient sample size, which can easily lead to excessively large variance in the bias estimate. Conversely, a window length greater than 60 days introduces cross-seasonal climate differences, which is not conducive to reflecting the error characteristics under the current weather conditions. Therefore, a 30-day window strikes a balance between statistical stability and timeliness representativeness. After completing the statistical analysis of the systematic bias mean, bias correction is performed on the current global atmospheric element gridded prediction field for each variable, level, and grid point. The correction result is then used as the assimilated background field. The correction expression is: ,in, This represents the assimilation background field after bias correction. This represents the gridded prediction field of global atmospheric elements output by the AI ​​weather forecasting model at the current moment. This represents the mean of systematic deviations obtained in advance by pre-statistics based on variables, isobars, and seasons.

[0024] In one specific embodiment, the process of inputting ERA5 global atmospheric reanalysis gridded data into an AI weather forecasting model to perform autoregressive iterative inference and obtain a global atmospheric element gridded prediction field can specifically include the following steps: The input state tensor of the AI ​​weather forecasting model is constructed based on the upper atmospheric element field and the surface atmospheric element field in the ERA5 global atmospheric reanalysis grid data. The input state tensor is fed into the AI ​​weather forecast model to perform autoregressive iterative inference, and the atmospheric state tensor output by each inference step is used as the input state tensor for the next inference step. The process is iterated until the forecast lead time, and the global atmospheric element grid prediction field is obtained.

[0025] Specifically, in the forecast initiation phase, the complete atmospheric state information of the current time t of the ERA5 global atmospheric reanalysis grid data is read, and the ERA5 global atmospheric reanalysis grid data is split into two parts, the upper atmospheric element field and the surface atmospheric element field, according to the established data organization method of the AI ​​weather forecast model. The upper atmospheric element field is used to describe the three-dimensional thermodynamic structure of the free atmosphere, and the surface atmospheric element field is used to describe the atmospheric state near the ground and sea level. Then, variable filtering, hierarchical alignment, grid registration and tensor encapsulation are performed on the two parts of the data respectively to jointly construct the input state tensor of the AI ​​weather forecast model. The tensor dimensions of the upper-level atmospheric element field are set as follows: 13 standard isobars × 5 types of variables × 721 zonal grid points × 1440 meridional grid points. The 13 standard isobars are 1000 hPa, 925 hPa, 850 hPa, 700 hPa, 600 hPa, 500 hPa, 400 hPa, 300 hPa, 250 hPa, 200 hPa, 150 hPa, 100 hPa, and 50 hPa. The 5 types of upper-level variables are air temperature T, zonal wind speed u, meridional wind speed v, specific humidity q, and geopotential height z. The tensor dimensions of the surface atmospheric element field are set as follows: 4 types of variables × 721 zonal grid points × 1440 meridional grid points. The 4 types of surface variables are: 10-meter zonal wind speed u. 10 10-meter meridional wind speed v 10 2 meters temperature T 2m And sea level pressure (MSL). Through the above organization method, the upper-air atmospheric element field and the surface atmospheric element field are uniformly encapsulated into an atmospheric state tensor at the same moment. The input state tensor is fed into the AI ​​weather forecasting model for forward inference, and the AI ​​weather forecasting model then calculates the atmospheric state tensor for the current moment. The atmospheric state is used for feature encoding, spatiotemporal correlation extraction, and joint mapping of high and low level variables to output the next forecast time. The atmospheric state tensor, where the time step is... The atmospheric state tensor for the next forecast time is fixed at 1 hour and denoted as . In order to achieve continuous and timely forecasting, after obtaining... Then, the atmospheric state tensor output from this inference is directly used as the input state tensor for the next round of inference, continuing the forward calculation for the next hour, forming an autoregressive iterative chain of "input-inference-output-backfill-re-inference". That is, the first inference uses the input state tensor constructed from ERA5 global atmospheric reanalysis grid data; the second inference uses the atmospheric state tensor obtained from the first inference as input; the third inference uses the atmospheric state tensor obtained from the second inference as input; and subsequent inferences are performed in the same loop until the set forecast lead time is reached. Step size With a forecast period of 1 hour and a lead time covering the next 0 to 72 hours, a set of prediction fields for 72 forecast times can be generated through continuous rolling inference. During the autoregressive iteration process, the dimension of the atmospheric state tensor output at each step remains consistent with that of the input state tensor, thus eliminating the need to modify the data structure between adjacent inference steps. As the iteration progresses, the model outputs hourly gridded predictions of meteorological elements across the globe for the corresponding time, forming a global atmospheric element gridded prediction field.

[0026] In one specific embodiment, the process of performing step 102 may specifically include the following steps: Interpolation processing is performed on the first observation data of the ground automatic weather station to obtain the ground automatic weather station observation data; The second radiosonde observation data is standardized to obtain radiosonde observation data. CO2 slice cloud detection and data removal are performed on the raw radiance data of geostationary meteorological satellites to obtain geostationary meteorological satellite radiance data; The signal-to-noise ratio threshold test and data removal were performed on the raw VAD wind profile data of the Doppler radar to obtain the Doppler radar VAD wind profile data. An effective multi-source observation vector is constructed based on observation data from ground-based automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data.

[0027] Specifically, while the AI ​​weather forecast model generates the assimilation background field, the acquisition, quality control, and unified compilation of multi-source real-time observation data are initiated in parallel. Hourly first observation data from automatic weather stations are received, including surface temperature, relative humidity, station air pressure, zonal wind speed, meridional wind speed, surface precipitation, visibility, and surface dew point temperature. First, extreme value checks and time series continuity checks are performed on the single-station measured data. Measured values ​​with a deviation of more than 10K between adjacent time intervals for temperature and more than 20m / s between adjacent time intervals for wind speed are marked as invalid and discarded. For the measured data after preliminary quality control, for missing values, neighboring station observations are searched within a 150km radius of each station, and Cressman distance-weighted interpolation is performed to fill the gaps. The interpolation weights are: ; in, This indicates the search radius, with a value of 150km. Indicates the first The station and the first The distances between adjacent stations are expressed in km. Preliminary quality-controlled measured data and interpolated values ​​are retained as observation data from surface automatic weather stations. For the second observation data of radiosondes, L-band radiosonde vertical profile data from 08:00 and 20:00 Beijing time are read. Five types of profile elements—air temperature, dew point temperature, zonal wind speed, meridional wind speed, and geopotential height—are organized around 13 standard isobaric layers. During standardization, temperature vertical lapse rate and wind speed vertical shear tests are first performed. When the temperature vertical gradient between adjacent radiosondes meets the requirements... When the vertical shear of zonal or meridional wind speed between adjacent layers exceeds 20 m / s / km, the corresponding layer's wind field data is also set to missing, resulting in radiosonde observation data with unified layers and outliers removed. For the raw radiance data from the geostationary meteorological satellite, the raw radiance field of the FY-4B geostationary meteorological satellite's GIIRS hyperspectral interferometer is read. The raw radiance field covers 1650 spectral channels in long-wave infrared and mid-wave infrared, and is set at a channel brightness temperature of 11μm BT. 11 As a CO2 slice cloud detection indicator, when BT 11At K < 220 K, the corresponding pixel is determined to be affected by high-level ice clouds or deep convective cloud tops. All 1650 channels of radiance data for that pixel are marked as invalid and removed, retaining only the radiance data from geostationary meteorological satellites corresponding to clear-sky pixels. For the raw Doppler radar VAD wind profile data, the radial velocity data at each elevation angle of the C-band Doppler weather radar are first read, and the vertical profile of the horizontal wind field is obtained through VAD algorithm inversion. Based on this, a signal-to-noise ratio threshold check and subsequent removal processing are performed. Radial data with a signal-to-noise ratio below 10 dB are directly removed. Simultaneously, velocity fuzziness de-folding processing is performed on the VAD inversion results, and profile level data with a VAD fitting residual root mean square exceeding 3 m / s are removed, resulting in Doppler radar VAD wind profile data that passes the validity screening. After interpolation verification, hierarchical standardization, cloud detection and removal, and signal-to-noise ratio threshold verification were completed for the four types of observation data, the ground automatic weather station observation data, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data were sequentially spliced ​​according to a unified observation coding order to construct an effective multi-source observation vector y that passed quality control. obs ∈R M Where M represents the total number of valid observation data of the four types, forming a unified observation constraint input.

[0028] In one specific embodiment, after constructing the effective multi-source observation vector, the method further includes: The first instrument error variance of the ground automatic weather station, the second instrument error variance of the radiosonde, the third instrument error variance of the geostationary meteorological satellite, and the fourth instrument error variance of various instruments of the Doppler radar were obtained respectively. Based on the variances of the first, second, third, and fourth instrument errors, the diagonal elements of each observation element in the effective multi-source observation vector are assigned values ​​to construct the observation error covariance matrix.

[0029] Specifically, the observation error covariance matrix is ​​in diagonal form, used to characterize the comprehensive dispersion of the measurement errors and representative errors of various observation data, and serves as the quantitative basis for the weights of various observation constraints in the assimilation calculation. The first instrument error variance of the automatic weather station, the second instrument error variance of radiosondes, the third instrument error variance of geostationary meteorological satellites, and the fourth instrument error variance of various Doppler radar instruments are obtained respectively, and error variance lookup relationships are established item by item according to different observation types and different observation elements. For the observation data of the automatic weather station, a fixed prior error variance assignment method is adopted, where the first instrument error variance of the surface temperature observation element T is taken as 1.0K. 2 The variance of the first instrument error for the zonal and meridional wind speed observation elements u / v is taken as 1.5m. 2 / s 2The first instrument error variance of the atmospheric pressure observation element P at the station is taken as 2.0 hPa. 2 For radiosonde observation data, the variance of the second instrument error for the temperature observation element T is taken as 0.5K. 2 The variance of the second instrument error for the zonal and meridional wind speed observation elements u / v is taken as 1.0 m. 2 / s 2 For geostationary meteorological satellite radiance data, the variance of the third instrument error for the atmospheric top brightness temperature (BT) observation element is taken as 1.2K. 2 For Doppler radar VAD wind profile data, the fourth instrument error variance of the zonal and meridional wind speed observation elements u / v is taken as 2.0 m. 2 / s 2 All error variance values ​​are set according to the nominal accuracy of the various observation instruments and in conjunction with the operational error configuration of the corresponding observation type in the meteorological operational assimilation system. This reflects the measurement noise level of different instruments and the comprehensive uncertainty introduced by the mismatch between observation representativeness and model resolution. Based on the arrangement order of the effective multi-source observation vectors, the variance is assigned diagonally to each observation element in the effective multi-source observation vectors. Since the effective multi-source observation vector y... obs It is formed by sequentially splicing observation data from ground-based automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data. Therefore, the diagonal elements of the observation error covariance matrix R are also written in the same order. When the current observation element comes from a ground-based automatic weather station and the observation type is temperature, the corresponding diagonal element is assigned a value of 1.0K. 2 When the current observed element originates from a ground-based automatic weather station and the observation type is zonal or meridional wind speed, the corresponding diagonal element is assigned a value of 1.5 m. 2 / s 2 When the current observed element originates from a ground-based automatic weather station and the observation type is station pressure, the corresponding diagonal element is assigned a value of 2.0 hPa. 2 When the current observed element originates from radiosonde and the observation type is temperature, the corresponding diagonal element is assigned a value of 0.5K. 2 When the current observation element originates from radiosonde and the observation type is zonal or meridional wind speed, the corresponding diagonal element is assigned a value of 1.0 m. 2 / s 2 When the current observed element originates from a geostationary meteorological satellite and the observation type is brightness temperature, the corresponding diagonal element is assigned a value of 1.2 K. 2 When the current observed element originates from the Doppler radar VAD wind profile and the observation type is zonal or meridional wind speed, the corresponding diagonal element is assigned a value of 2.0 m. 2 / s 2By assigning values ​​element-wise, each valid observation in the observation vector can find a unique corresponding error variance position in the observation error covariance matrix. This allows observation information from different sources, with different dimensions, and different levels of precision to participate in weight calculation in a unified matrix form within the assimilation cost function. Since this scheme adopts the three-dimensional variational assimilation standard assumption, the error correlation between different observations is ignored when constructing the observation error covariance matrix. Only the distribution of the error variance of each observation element on the main diagonal is retained, while off-diagonal elements are uniformly set to zero, forming a diagonal matrix form of the observation error covariance matrix R. The physical meaning of this approach is that it assumes that the observation errors of various types of observations—including automatic weather stations, radiosondes, geostationary meteorological satellites, and Doppler radar—are independent of each other, and the cross-error correlation between different observation elements within different observation sources is not explicitly modeled. Therefore, it can reduce the matrix storage size and inversion complexity during the assimilation calculation process, while simultaneously meeting the computational efficiency and stability requirements of operational three-dimensional variational systems. In this construction method, if the total length of the effective multi-source observation vector is M, then the observation error covariance matrix can be represented as an M×M diagonal matrix, where the main diagonal elements are composed of the variances of the four types of instrument errors in the order of the observed elements, and the remaining positions of the matrix are all zero. Therefore, the final result includes the effective multi-source observation vector y that passes quality control. obs ∈R M And the corresponding observation error covariance diagonal matrix R.

[0030] In one specific embodiment, the process of performing step 103 may specifically include the following steps: The background constraint term is constructed based on the assimilated background field, and the observation constraint term is constructed based on the effective multi-source observation vector and the observation error covariance matrix. The background constraint term and the observation constraint term are summed to obtain the first cost function. The atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable to obtain the second cost function. The optimal atmospheric analysis field is obtained by solving the second cost function through iterative minimization.

[0031] Specifically, the bias-corrected assimilation background field is introduced as a priori state constraint into the three-dimensional variational assimilation framework, where the assimilation background field vector is denoted as x. b x bThis represents the background state obtained at the current assimilation time after systematic bias correction based on the prediction results of the AI ​​weather forecast model. The background state covers all grid points, all vertical levels, and corresponding meteorological variables within the WRF simulation domain. Simultaneously, the atmospheric analysis state vector to be solved is denoted as x, representing the optimal atmospheric state to be obtained under the combined effects of observational and background constraints. The background error covariance matrix is ​​denoted as B, which characterizes the spatial discreteness of the background field error and the statistical correlation between different state components. Based on the above definitions, a background constraint term is constructed. The physical meaning of the background constraint term is to measure the deviation of the atmospheric analysis state vector x from the assimilated background field vector x. b The degree of deviation should be weighted according to the statistical structure of the background error covariance matrix B to avoid the analysis field deviating too far from the background field. Then, an effective multi-source observation vector y, composed of data from ground automatic weather stations, radiosonde observations, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data, is introduced. obs , where y obs This represents the set of all valid observations retained after quality control; simultaneously, an observation error covariance matrix R is introduced, which characterizes the magnitude of the error of each observation element; further, an observation operator H is defined, which represents the transformation relationship mapping the atmospheric analysis state vector x in the model space to the observation space. Based on y obs R and H construct the observation constraint term. The physical meaning of the observation constraint term is to measure the degree of consistency between the simulated observations obtained from the analysis field mapping and the actual observations, and it is weighted according to the observation error covariance matrix. The first cost function is: ; in, Represents the atmospheric analysis state vector to be determined. The first cost function, the first term This is a background constraint term used to limit the deviation between the analysis field and the background field. The second term... The observation constraint term limits the deviation between the simulated and actual observations of the analysis field. Therefore, by summing the background constraint term and the observation constraint term, the two types of constraints, "approaching the background field as closely as possible" and "fitting the observed values ​​as closely as possible," are unified into a single optimization objective. Since the first cost function directly uses... The background error covariance matrix is ​​the independent variable. The relevant structures will directly participate in the optimization process, therefore incremental control variables are introduced. Perform a standardized transformation on the state variables, and the transformation relationship is written as follows: ,in, Represents the background error covariance matrix The square root operator is used to map standardized control variables back to the actual state space; This represents the incremental control variable, used to describe the expression of the analysis increment in the background error standardized space. After substituting the above transformation into the first cost function, the background constraint term changes from its original form... Transform into At the same time, the state variables in the observation constraint terms are also replaced. This yields incremental control variables. The second cost function for the independent variable: ; in, This represents the value of the second cost function. This represents the background constraint term in the control variable space. This represents the observation space residual vector. The second cost function is then subjected to an adjustment with respect to the incremental control variables. The iterative minimization method is used, employing the L-BFGS-B quasi-Newton method. The maximum number of iterations is set to 150, and the gradient convergence criterion is set to... ,in, This represents the second cost function relative to the incremental control variable at the current iteration step. gradient vector, This represents the L2 norm of the current gradient vector. This represents the gradient vector at the initial iteration time. This represents the L2 norm of the initial gradient vector. The optimal incremental control variable is obtained when the second cost function satisfies the convergence condition. Then, the optimal incremental control variables are mapped back to the actual state space through inverse transformation and superimposed onto the assimilated background field to obtain the optimal atmospheric analysis field. ,in, This represents the optimal atmospheric analysis field. This represents the optimal incremental control variable that minimizes the second cost function. Thus, the optimal atmospheric analysis field maintains the large-scale continuous structure of the assimilated background field while incorporating real-time constraint information provided by effective multi-source observation vectors, resulting in an analysis that satisfies both background statistical regularity and observation consistency requirements.

[0032] In one specific embodiment, the process of transforming the atmospheric analysis state vector in the first cost function into incremental control variables to obtain the second cost function, and then solving the second cost function through iterative minimization to obtain the optimal atmospheric analysis field, can specifically include the following steps: Based on the background error covariance matrix, the atmospheric analysis state vector to be determined in the first cost function is transformed into incremental control variables, resulting in a second cost function with the incremental control variables as independent variables; whereby the second cost function is: ; The second cost function value. Incremental control variables, Background constraint term in the second cost function. Assimilation of the background field, Background error covariance matrix The square root operator, Nonlinear observation operator : Effective multi-source observation vector Observation error covariance matrix The inverse matrix; The second cost function is iteratively minimized with respect to the incremental control variables to obtain the optimal incremental control variables. The optimal incremental control variables are then inversely transformed using the background error covariance matrix and superimposed onto the assimilation background field to obtain the optimal atmospheric analysis field.

[0033] Specifically, based on the background error covariance matrix For the atmospheric analysis state vector to be determined in the first cost function Perform control variable transformation, where the assimilation background field is denoted as... , representing the prior state vector formed after prediction and bias correction by the AI ​​weather forecast model, and the square root operator of the background error covariance matrix is ​​denoted as . , represents the mapping operator that satisfies the background error statistical decomposition relation, and the incremental control variable is denoted as . , representing the expression of the analytical increment in the standardized error space, therefore, the atmospheric analytical state vector to be determined and the incremental control variables satisfy the following relationship. This relationship. The physical meaning of the above relationship is that it represents the analytical increment in the actual state space. It is decomposed into two parts: "background error statistical structure" and "dimensionless control variables". The background error statistical structure consists of... The magnitude of the control variable is borne by Taking responsibility. Substituting the transformation relationship into the first cost function, we obtain the incremental control variable. The second cost function is the independent variable. Since the second cost function has been transformed into a function about the incremental control variable... Since this is an unconstrained minimization problem, the second cost function can be directly solved by iterative minimization. In each iteration, the value of the second cost function and gradient information are calculated based on the current incremental control variable, and the search direction and step size are updated using the L-BFGS-B quasi-Newton method, so that the second cost function gradually decreases along the descent direction until the preset convergence condition is met or the maximum number of iterations is reached. When the iteration terminates, the optimal incremental control variable is obtained. In obtaining Then, the square root operator of the background error covariance matrix is ​​used to perform an inverse transformation on the optimal incremental control variable, restoring the optimal control quantity in the standardized space to the analysis increment in the actual state space. This is then superimposed with the assimilated background field to obtain the optimal atmospheric analysis field. ,in, This represents the optimal atmospheric analysis field. This represents the optimal incremental control variable. This represents the analytical increment obtained after the inverse transformation.

[0034] like Figure 3 In one specific embodiment, the process of performing step 104 may specifically include the following steps: 301. Extract the atmospheric state variables of each grid point in the boundary relaxation region of the WRF model from the optimal atmospheric analysis field and write them into the WRF model side boundary drive file; 302. Based on the WRF model side boundary driving file, drive the WRF model to complete the integral operation and output the meteorological element grid forecast field of surface temperature, surface wind speed, surface air pressure, precipitation, cloud cover and convective available potential energy.

[0035] Specifically, the optimal atmospheric analysis field The side boundary driving information can be directly identified and invoked in the WRF model, where, This represents the optimal analysis result obtained after three-dimensional variational assimilation of AI weather forecast model predictions and multi-source real-time observation data, with a spatial fundamental resolution of [missing value]. The vertical direction includes 13 standard isobaric layers and surface meteorological elements. Therefore, before writing the data into the WRF model, the optimal atmospheric analysis field is first analyzed using the metgrid module in WPS. Two types of coordinate transformations were performed: horizontal and vertical. Horizontally, bicubic spline interpolation was used to map the global latitude and longitude grid field to the Lambert conformal conic projection grid of the WRF simulation domain. The horizontal resolution of the WRF simulation domain was 3 km, and the grid size was 750×600. Vertically, based on the n-mixed σ-pressure coordinate structure of the WRF model, state variables such as temperature, zonal wind speed, meridional wind speed, specific humidity, and geopotential height from 13 standard isobaric layers were mapped to 50 vertical isobaric layers using logarithmic barometric linear interpolation. noodle The densification layer also needs to be vertically extrapolated and corrected based on actual surface air pressure to ensure that the near-surface thermal and dynamic structure matches the complex terrain conditions. After coordinate transformation, atmospheric state variables are extracted from each grid point of the relaxation region around the WRF simulation domain from the transformed optimal analysis field. The width of the relaxation region is 5 grids, corresponding to a width of 15km at 3km resolution, which is used to achieve a smooth transition between the outer boundary constraints and the inner free integration region during model integration. For each hour... The atmospheric state values ​​of each grid point within the boundary relaxation region are denoted as follows: ,in, Indicates time The boundary state quantities include driving variables such as temperature, wind field, humidity, and pressure at the boundary grid points, and will... This is written into the corresponding time-boundary state field of the WRF model's side boundary driving file wrfbdy_d01. Simultaneously, to enable the WRF model to complete continuous-time evolution between adjacent integer boundary times, the time trend term of the boundary state variables is calculated based on the atmospheric state variables of two adjacent boundary times. The formula is as follows: ; in, Indicates time The boundary state quantity time trend term, Indicates the next hour. Boundary state quantities, The boundary update interval is 3600s. After the time trend term is synchronously written into the boundary hyperpotential field of the corresponding time in the wrfbdy_d01 file, the WRF model enables the recovery of boundary constraints at any intermediate time between two integer boundary conditions through linear time interpolation, forming a boundary-driven link that is dynamically updated hourly. After the side boundary driving file is written, the wrfbdy_d01 file is used to drive the WRF model to complete the integration operation. The WRF model adopts ARW Dynamics Core version 4.5, using an Arakawa-C staggered grid in the horizontal direction and a 50-layer η-mixed σ-pressure coordinate system in the vertical direction. In terms of physical parameterization, the Morrison two-parameter ice phase scheme is used for microphysical processes, the MYJ scheme is used for boundary layer processes, the Noah-MP scheme is used for land surface processes, and the RRTMG scheme is used for long-wave and short-wave radiation. The convection process is explicitly analyzed at a resolution of 3km. During mode integration, WRF does not directly and rigidly replace the boundary state variables with internal grid points. Instead, it applies a gradually decaying relaxation forcing term to each boundary grid point within the four-sided boundary relaxation regions to suppress boundary abrupt changes and avoid generating spurious gravity waves. The corresponding boundary relaxation forcing expression is: ; in, Indicates the first The forcing term imposed on each boundary relaxed lattice point, This represents the grid point index within the boundary relaxation region, with values ​​ranging from 1 to 5. This represents the boundary state quantity at the most recent integer time. This represents the boundary time trend term at the corresponding moment. This represents the time difference between the current integration time and the nearest integer boundary time, with a value ranging from 0 to 3600 seconds. This indicates that the current integration time of the WRF model is at time n. Mode state values ​​at each boundary relaxation grid point This represents the boundary relaxation coefficient, used to control the strength of the external boundary constraints on different boundary grid points. The relaxation coefficient adopts a Logistic function-type spatial decay form, and its expression is: ; in, Represents the natural constant. The time corresponds to the outermost boundary grid point. The time corresponds to the innermost grid point in the relaxation region. According to this formula, the outermost boundary grid point has a larger relaxation weight, while the relaxation weight of the innermost boundary grid point is significantly reduced. This allows the assimilated boundary field to smoothly transition from the outside to the WRF free integration domain, maintaining the constraint of boundary-driven evolution on the model while reducing the risk of numerical instability caused by boundary perturbations propagating inward. In the specific time-progression process, the WRF model adopts a third-order Runge-Kutta time integration scheme with a large step size. Take 18s, and use a small step size for explicit integral of the sound wave term splitting. Take 3 seconds, WRF self-reporting time Continuous integration to h, and reads the dynamically updated wrfbdy_d01 boundary condition file every 3600 seconds. As the integration progresses, the model, under the combined influence of microphysics, boundary layer, radiation, land surface, and dynamic core, generates a high-resolution gridded forecast field of meteorological elements covering the Chinese region hourly. In the final output elements, surface temperature is calculated from 2-meter temperature... The unit is K; surface wind speed is expressed as zonal wind speed per 10 meters. and 10-meter meridional wind speed Synthetic characterization, in which, and The units are all m / s; the ground air pressure is from The unit is Pa; precipitation is expressed as hourly cumulative precipitation. The unit for cloud cover is mm; cloud cover is expressed as a cloud cover fraction. This indicates that the value ranges from 0 to 1; the convective effective potential energy is determined by... It indicates that the unit is J / kg.

[0036] In one specific embodiment, the process of extracting atmospheric state variables from each grid point in the WRF model boundary relaxation region from the optimal atmospheric analysis field and writing them into the WRF model side boundary drive file can specifically include the following steps: In the optimal atmospheric analysis field, global latitude and longitude grid data are interpolated to the WRF model, and logarithmic pressure linear interpolation is performed with the pressure value corresponding to each η-mixed σ-pressure coordinate layer as the interpolation node, and then transformed to the η-mixed σ-pressure coordinate layer of the WRF model. Atmospheric state variables and corresponding trend terms are extracted from each grid point in the WRF model boundary relaxation region using the η-mixed σ-pressure coordinate layer of the WRF model, and then written into the WRF model side boundary driving file.

[0037] Specifically, coordinate transformation and grid remapping were performed on the optimal atmospheric analysis field for the WRF model grid system. The optimal atmospheric analysis field adopts a global 0.25°×0.25° latitude and longitude grid organization and includes upper atmospheric state information from 13 standard isobaric layers. Therefore, the global latitude and longitude grid data was transformed into the regional projection grid and vertical η-mixed σ-pressure coordinate system used by the WRF model using the metgrid module in WPS. In the horizontal direction, the global latitude and longitude grid field was interpolated to the Lambert conformal conic projection grid of the WRF simulation domain. The center of the WRF simulation domain was set to 105°E and 35°N, the horizontal resolution was set to 3km, and the horizontal grid size was set to 750×600. Bicubic spline interpolation was used to complete the horizontal resampling of the scalar field and wind field. At the same time, height correction was performed on the terrain-blocked areas to ensure that the transformed atmospheric field is consistent with the complex regional terrain. After horizontal resampling, a coordinate transformation is performed in the vertical direction. The WRF model uses η-hybrid σ-pressure coordinates in the vertical direction, and the reference pressure profile is written as follows: ,in, Represents the η coordinate layer The corresponding reference pressure value, in hPa; The dimensionless η-mixing σ-pressure coordinate layer in the vertical direction of the WRF model represents the WRF model. This represents the ground reference pressure, with a value of 1013.25 hPa. The top reference pressure of the model is represented by 50 hPa. Based on the aforementioned reference pressure profile, 50 reference pressure nodes corresponding to each of the n layers are obtained. Then, using the reference pressure values ​​corresponding to the 50 n layers as interpolation nodes, logarithmic barometric linear interpolation is performed on atmospheric state variables such as air temperature, zonal wind speed, meridional wind speed, specific humidity, and geopotential height on 13 standard isobaric layers, mapping them to the 50 η-mixed σ-baric coordinate layers of the WRF model. The near-surface densification region is set at... Within the area, the near-surface densification zone comprises 20 layers, with the lowest layer approximately 20 meters thick. Furthermore, additional layers based on actual ground air pressure are added within the near-surface layer. The vertical extrapolation correction, where, This represents the actual surface air pressure, expressed in hPa. Vertical extrapolation correction is used to improve the matching accuracy of near-surface thermal and dynamic fields in complex terrain areas. After horizontal interpolation, projection transformation, generation of η-layer pressure nodes, and logarithmic pressure linear interpolation, the optimal atmospheric analysis field is uniformly converted into a regional gridded state field that can be directly identified by the WRF model. Atmospheric state variables of each grid point in the four boundary relaxation zones are extracted from the converted regional gridded state field. The width of the boundary relaxation zone is set to 5 grids, corresponding to a width of 15km under a 3km grid resolution. The boundary relaxation zone covers all boundary grids on the east, west, south, and north sides of the WRF simulation domain, and is applied at each hour. Extract the atmospheric state values ​​of each grid point within the boundary relaxation region. ,in, Indicates time The boundary state quantities include the air temperature, wind field, humidity, geopotential height, and related driving variables of each η layer at the boundary grid points. The units depend on the variable type: K for air temperature, m / s for wind speed, kg / kg for specific humidity, and gpm for geopotential height. Write the boundary state variables fields at the corresponding time step into the WRF model's side boundary driving file wrfbdy d01. To enable continuous-time interpolation between two adjacent integer boundary fields in the WRF model, and to calculate the trend term corresponding to the boundary state variables, the trend term is represented in the form of a time trend term, written as: ; in, Indicates time The boundary state quantity trend term, which is the time trend term of the boundary state quantity, is expressed in units of the corresponding state variable units divided by seconds. Indicates the next hour. Boundary state quantities; Indicates the current hour. Boundary state quantities; This represents the boundary update interval, with a value of 3600s. The calculated boundary state variables are then written into the boundary trend term field of the corresponding time step in the wrfbdy d01 file. This allows the WRF model to perform linear time interpolation based on the current point-in-time boundary state variables and the boundary trend term within the integration step, forming a dynamic boundary driving file that can directly drive the evolution of the WRF model's side boundaries.

[0038] Figure 4This is a schematic diagram of the vertical structure of the standard isobaric layer in the optimal atmospheric analysis field. The vertical axis in the figure corresponds to the standard isobaric layers with pressure values ​​of 50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, and 1000 hPa from top to bottom. The top is labeled "Model Top" corresponding to a model top pressure of 50 hPa, and "Ground" corresponding to a ground pressure of 1000 hPa. This vertical structure serves as the source coordinate layer for logarithmic pressure linear interpolation.

[0039] Figure 5 This diagram illustrates the process of performing logarithmic barometric linear interpolation from the standard isobaric layer to the WRF model's η-mixed σ-baric coordinate layer. Solid circles represent grid points in the source layer (standard isobaric layer), solid triangles represent grid points in the target layer (η-coordinate layer), and dashed arrows represent the mapping relationship of interpolating atmospheric state variables from the source layer to the target layer using ln(p) as the interpolation node. One source layer can provide interpolation contributions to multiple target layers, demonstrating that the η-coordinate layer has a higher vertical resolution than the standard isobaric layer, especially with a denser distribution in the near-surface layer.

[0040] Figure 6 This is a schematic diagram of the global latitude and longitude grid distribution of the optimal atmospheric analysis field, including a global 0.25°×0.25° resolution regular latitude and longitude grid. The grid spacing in the horizontal and vertical directions is uniform, corresponding to a spatial resolution of 0.25° in the latitudinal and longitudinal directions, respectively. The total number of grid points is 721 latitude grid points × 1440 longitude grid points. This regular latitude and longitude grid field is used as the source grid points for horizontal interpolation. After bicubic spline interpolation, it is transformed to a 3 km resolution grid under the Lambert conformal projection coordinates of the WRF model.

[0041] In one specific embodiment, the process of executing the WRF model to complete the integral operation based on the WRF model side boundary driving file and outputting the meteorological element grid forecast field of surface air temperature, surface wind speed, surface air pressure, precipitation, cloud cover and convective available potential energy can specifically include the following steps: Based on the WRF model side boundary driving file, the spatial attenuation coefficient of each grid point in the WRF model boundary relaxation region is calculated, and the atmospheric dynamic equations of the WRF model are constructed according to the spatial attenuation coefficient and atmospheric state variables. The atmospheric dynamic equations of the WRF model are driven by a third-order Runge-Kutta time integration scheme to perform numerical integration on the ARW dynamic core of the WRF model, outputting a gridded forecast field of meteorological elements such as surface temperature, surface wind speed, surface pressure, precipitation, cloud cover, and convective available potential energy.

[0042] Specifically, the boundary state variables and boundary trend terms at each hour are read from the WRF model. Within a 5-grid boundary relaxation region surrounding the simulation domain, a boundary constraint strength that decreases with spatial location is established for each boundary grid point. This gradually transitions the outer boundary drive to the inner free integration region, preventing abrupt boundary changes from directly entering the model and inducing spurious gravity waves or numerical oscillations. The boundary relaxation region is then used to... The spatial attenuation coefficient of each grid point is denoted as ,in, Indicates the first The relaxation weight coefficients corresponding to each boundary relaxation grid point are dimensionless. This represents the grid point number in the boundary relaxation region, with values ​​ranging from 1 to 5. Corresponding to the outermost boundary grid point, This corresponds to the innermost relaxed lattice point. The spatial decay coefficient is calculated using a Logistic function, with the following expression: ; in, Let represent the natural constant. According to this formula, the outermost boundary grid points have larger relaxation weights, while the relaxation weights decrease as the grid points move closer to the model interior. Therefore, the boundary forcing exhibits a spatially decreasing distribution, forming a smooth constraint band from strong to weak. Based on the boundary state variables, the boundary trend term, and the spatial decay coefficient, a boundary relaxation forcing term is constructed within the boundary relaxation region. This boundary relaxation forcing term is denoted as . Its expression is: ; in, Indicates the first The atmospheric state forcing term for each boundary relaxed grid point has the same unit as the corresponding state variable; Indicates the current hour. The units of boundary state quantities depend on the type of variable. When the variable is temperature, the unit is K; when the variable is wind speed, the unit is m / s; when the variable is air pressure, the unit is Pa; and when the variable is specific humidity, the unit is kg / kg. This represents the boundary trend term at the current hour, which is the rate of change of the boundary state quantity with respect to time, and the unit is the corresponding variable unit divided by seconds; This represents the time difference between the current integration time and the nearest integer boundary time, expressed in seconds. This indicates that the current integration time of the WRF model is at time n. The model state values ​​at each boundary relaxation grid point. The physical meaning of the aforementioned forcing term is to first recover the boundary target state that should be reached at the current integration time using the boundary state variables and the boundary overpotential term, then subtract the target state from the current model state to form the boundary residual, and finally use the spatial attenuation coefficient to spatially weight the boundary residual, so that the grid points outside the boundary relaxation region are more closely aligned with the state provided by the driving file, while the grid points inside retain greater model free evolution capability. The boundary relaxation forcing term is incorporated into the atmospheric dynamic equations of the ARW dynamic core of the WRF model. The WRF model uses ARW dynamic core version 4.5, whose atmospheric dynamic equations include the mass conservation equation, momentum equation, thermodynamic equation, and water vapor continuity equation, and are solved discretized on the Arakawa-C staggered grid using finite difference form. In the boundary relaxation region, various atmospheric state variables, such as temperature, wind field, water vapor, and air pressure, are simultaneously affected by internal dynamic processes, physical processes, and boundary relaxation forcing terms during numerical integration. Therefore, the spatial attenuation coefficient is understood as the weight modulation of boundary constraints in the dynamic equations, while the boundary state variables are understood as the target states of external driving sources in the dynamic equations. After the boundary forcing term is added, the WRF model no longer relies on a static 6-hour boundary update method. Instead, it reads new boundary state variables and boundary trend terms every 3600 seconds based on the hourly updated WRF model side boundary driving file, thereby continuously correcting the boundary evolution trajectory and ensuring that the boundary conditions of the WRF regional model remain time-synchronized with the optimal atmospheric analysis field assimilated and corrected by the AI ​​meteorological forecast model and multi-source observations. After completing the equations construction, a third-order Runge-Kutta time integration scheme is used to drive the ARW dynamic core of the WRF model to perform numerical integration. The WRF model time advancement uses a third-order Runge-Kutta explicit integration scheme, with the large step size denoted as... The value is 18s, where, The time step size of the master dynamic equation is expressed in seconds (s). For the fast acoustic term, a split explicit small-step scheme is adopted, with the small step size denoted as . The value is 3s, where, This represents the acoustic wave substep size, measured in seconds (s). The ratio of the large step size to the small step size is 6:1, which meets the ARW power core's requirements for acoustic wave splitting integral stability. Simultaneously, the horizontal diffusion coefficient is denoted as... The value is 250m 2 / s, where This represents the diffusion coefficient in the horizontal direction of the model, in units of m. 2 / s; vertical mixing is diagnosed using the turbulent kinetic energy equation in the MYJ boundary layer scheme. The model is based on the reporting time... Start continuous integration until h, where The time in time indicates the start time of the model, expressed in hours (h). As the third-order Runge-Kutta time integral progresses cyclically, the ARW dynamic core sequentially calculates the contributions of processes such as advection, pressure gradient force, gravity, turbulent diffusion, microphysical transformation, radiative heating, boundary relaxation forcing, and surface exchange to each state variable at each time step. These contributions are then integrated and updated into the atmospheric state field at the next integration time, forming a continuous three-dimensional numerical forecast covering the entire simulation domain. After integration, the WRF model outputs gridded forecast fields of meteorological elements at hourly intervals, where surface temperature is calculated from 2-meter temperature... express, The unit is K; ground wind speed is the zonal wind speed at 10 meters per second. and 10-meter meridional wind speed Common characteristics, among which, and The units for both are m / s, and the surface wind speed can be obtained by combining the two; the surface air pressure is... express, The unit is Pa; precipitation is calculated from hourly cumulative precipitation. express, The unit is mm; cloud cover is expressed as cloud cover fraction. express, It is a dimensionless quantity, ranging from 0 to 1; the convective effective potential energy is... express, The unit is J / kg.

[0043] Figure 7 This is a flowchart illustrating the process of constructing the atmospheric dynamic equations based on the WRF model's side boundary driving file. The process, from top to bottom, is as follows: Calculate the spatial attenuation coefficients of each grid point within the boundary relaxation region based on the WRF model's side boundary driving file; construct relaxation forcing terms in each equation based on the spatial attenuation coefficients and atmospheric state variables; append the relaxation forcing terms to the three branch equations: the mass continuity equation, the zonal momentum equation, and the meridional momentum equation; and merge these three branches into the thermodynamic equation and the water vapor continuity equation to complete the construction of the complete WRF model atmospheric dynamic equations, including the boundary relaxation forcing terms.

[0044] This invention applies targeted quality control to four types of heterogeneous observation data: ground-based automatic weather stations, radiosondes, geostationary meteorological satellites, and Doppler radar VAD wind profiles. It constructs an iterative minimization solution for the cost function and uses the optimal atmospheric analysis grid field, constrained by real-time multi-source observations, as the driving source for the WRF side boundary conditions. This achieves quantitative correction of the AI-predicted background field by observational information, effectively suppressing the propagation and accumulation of background field errors during the WRF integration process. The optimal atmospheric analysis grid field, after coordinate transformation, is written into the WRF model's side boundary condition file at hourly update intervals. Combined with a Logistic function-type spatial decay relaxation forcing mechanism, the WRF boundary condition update frequency is significantly higher than the traditional 6-hour static scheme. This effectively represents the rapid sub-hourly evolution of atmospheric conditions in the boundary region, thereby improving the WRF model's forecasting capability for high-impact weather events such as severe convection.

[0045] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0046] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0047] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A meteorological prediction method based on an AI meteorological prediction model and a WRF model, characterized in that, include: ERA5 global atmospheric reanalysis gridded data was input into an AI weather forecasting model for autoregressive iterative inference and bias correction to obtain the assimilated background field. An effective multi-source observation vector is constructed based on observation data from ground automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data. The optimal atmospheric analysis field is solved by using the assimilated background field as a priori constraint and the effective multi-source observation vector as an observation constraint. The atmospheric state variables of each grid point in the WRF model boundary relaxation zone are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary driving file to drive the WRF model to complete the integration operation and obtain the meteorological element grid point forecast field. 2.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 1, wherein, ERA5 global atmospheric reanalysis gridded data was input into an AI weather forecasting model for autoregressive iterative inference and bias correction to obtain the assimilated background field, including: Input the ERA5 global atmospheric reanalysis gridded data into the AI ​​weather forecasting model to perform autoregressive iterative inference and obtain the global atmospheric element gridded prediction field. The mean of the systematic bias of the global atmospheric element grid prediction field is statistically analyzed, and bias correction is performed on the global atmospheric element grid prediction field based on the mean of the systematic bias to obtain the assimilated background field. 3.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 2, characterized in that, By inputting ERA5 global atmospheric reanalysis gridded data into an AI weather forecasting model and performing autoregressive iterative inference, a global atmospheric element gridded prediction field is obtained, including: The input state tensor of the AI ​​weather forecasting model is constructed based on the upper atmospheric element field and the surface atmospheric element field in the ERA5 global atmospheric reanalysis grid data. The input state tensor is input into the AI ​​weather forecasting model to perform autoregressive iterative inference, and the atmospheric state tensor output at each step of the inference is used as the input state tensor for the next step of the inference. The process is iterated until the forecast lead time is reached to obtain the global atmospheric element grid prediction field. 4.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 1, wherein, An effective multi-source observation vector is constructed based on observation data from ground-based automatic weather stations, radiosonde observation data, geostationary meteorological satellite radiance data, and Doppler radar VAD wind profile data, including: Interpolation processing is performed on the first observation data from the automatic weather station to obtain the observation data from the automatic weather station. The second radiosonde observation data is standardized to obtain radiosonde observation data. CO2 slice cloud detection and data removal are performed on the raw radiance data of geostationary meteorological satellites to obtain geostationary meteorological satellite radiance data; The signal-to-noise ratio threshold test and data removal were performed on the raw VAD wind profile data of the Doppler radar to obtain the Doppler radar VAD wind profile data. An effective multi-source observation vector is constructed based on the observation data from the ground automatic weather station, the radiosonde observation data, the geostationary meteorological satellite radiance data, and the Doppler radar VAD wind profile data. 5.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 4, characterized in that, After constructing an effective multi-source observation vector, the following steps are also included: The first instrument error variance of the ground automatic weather station, the second instrument error variance of the radiosonde, the third instrument error variance of the geostationary meteorological satellite, and the fourth instrument error variance of various instruments of the Doppler radar were obtained respectively. Based on the first instrument error variance, the second instrument error variance, the third instrument error variance, and the fourth instrument error variance, the diagonal elements of each observation element in the effective multi-source observation vector are assigned values ​​respectively to construct the observation error covariance matrix. 6.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 5, characterized in that, Using the assimilation background field as a priori constraint and the effective multi-source observation vector as an observation constraint, the optimal atmospheric analysis field is solved, including: A background constraint term is constructed based on the assimilated background field, and an observation constraint term is constructed based on the effective multi-source observation vector and the observation error covariance matrix. The background constraint term and the observation constraint term are summed to obtain the first cost function. The atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable to obtain the second cost function. The optimal atmospheric analysis field is obtained by solving the second cost function through iterative minimization. 7.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 6, characterized in that, The atmospheric analysis state vector to be determined in the first cost function is transformed into incremental control variables to obtain the second cost function. The optimal atmospheric analysis field is obtained by iteratively minimizing the second cost function, including: Based on the background error covariance matrix, the atmospheric analysis state vector to be determined in the first cost function is transformed into an incremental control variable, resulting in a second cost function with the incremental control variable as the independent variable; wherein, the second cost function is: ； : second cost function value, : increment control variable, : background constraint term in the second cost function, : assimilate background field, : background error covariance matrix : square root operator of, : nonlinear observation operator, : effective multi-source observation vector, : inverse of observation error covariance matrix : inverse of The second cost function is iteratively minimized with respect to the incremental control variables to obtain the optimal incremental control variables. The optimal incremental control variables are then inversely transformed using the background error covariance matrix and superimposed onto the assimilation background field to obtain the optimal atmospheric analysis field. 8.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 1, wherein, The atmospheric state variables of each grid point in the WRF model boundary relaxation region are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary drive file. This drives the WRF model to complete the integration run, obtaining the meteorological element grid point forecast field, including: The atmospheric state variables of each grid point in the boundary relaxation region of the WRF model are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary drive file; The WRF model is driven by the side boundary driving file of the WRF model to complete the integral operation and output the meteorological element grid forecast field of surface temperature, surface wind speed, surface air pressure, precipitation, cloud cover and convective available potential energy. 9.The weather forecasting method based on the AI weather forecasting model and the WRF model according to claim 8, characterized in that, The atmospheric state variables of each grid point in the WRF model boundary relaxation region are extracted from the optimal atmospheric analysis field and written into the WRF model side boundary drive file, including: In the optimal atmospheric analysis field, global latitude and longitude grid data are interpolated to the WRF model, and logarithmic pressure linear interpolation is performed with the pressure value corresponding to each η-mixed σ-pressure coordinate layer as the interpolation node, and then transformed to the η-mixed σ-pressure coordinate layer of the WRF model. The atmospheric state variables and corresponding trend terms of each grid point in the WRF model boundary relaxation region are extracted by the η-mixed σ-pressure coordinate layer of the WRF model, and the atmospheric state variables and the trend terms are written into the WRF model side boundary driving file.

10. The weather forecasting method based on the AI ​​weather forecasting model and the WRF model according to claim 9, characterized in that, The WRF model is driven by the side boundary driving file of the WRF model to complete the integration operation and output a gridded forecast field of meteorological elements, including surface air temperature, surface wind speed, surface air pressure, precipitation, cloud cover, and convective available potential energy. Based on the WRF model side boundary driving file, the spatial attenuation coefficient of each grid point in the WRF model boundary relaxation region is calculated, and the atmospheric dynamic equations of the WRF model are constructed according to the spatial attenuation coefficient and the atmospheric state variables. The atmospheric dynamic equations of the WRF model are driven by a third-order Runge-Kutta time integration scheme to perform numerical integration on the ARW dynamic core of the WRF model, outputting a gridded forecast field of meteorological elements such as surface temperature, surface wind speed, surface air pressure, precipitation, cloud cover, and convective available potential energy.

Images