Satellite-ground collaborative convection inception identification and forecasting method based on physical identification

Abstract

The invention discloses a satellite-ground collaborative convection inception identification and forecasting method based on physical identification, and belongs to the technical field of convection inception identification and forecasting based on artificial intelligence. The prediction performance of the artificial intelligence prediction model can be effectively improved through more accurate prediction input samples, and the empty report rate and the false report rate of the prediction model are reduced. According to the method, based on multi-source observation and forecast data of satellites, ground-based radars and the like, high space-time matching is carried out on a satellite-based CI objective judgment result, a radar CI recognition result and ground minute-level rainfall data, and a more accurate high-quality convection inception annotation data set is formed; on the basis, a CI environment parameter matrix based on physical judgment is obtained by combining a numerical model product and utilizing a deep learning neural network model, and is fused with satellite parameters, model precipitation and the like based on a CI annotation data set, so that satellite-ground collaborative high-quality forecasting factors are provided for CI intelligent forecasting model training, and the empty report rate and the false alarm rate of CI are effectively reduced.

Description

Technical Field

[0001] The present invention belongs to the technical field of incipient convection recognition and prediction based on artificial intelligence, and in particular to a satellite-ground collaborative incipient convection recognition and prediction method based on physical identification. Background Art

[0002] The emergence of Convective Initiation (CI) is a sign of the beginning of severe convective activity. The identification and prediction of CI plays a key role in the early warning and prediction of severe convection. How to capture early information in the convection generation stage and accurately and effectively warn of CI is a problem of particular concern in nowcasting. Due to the small spatial scale and rapid development of new convection, it is difficult to effectively monitor and warn of CI by conventional meteorological observations. Usually, weather radars and geostationary meteorological satellites with high temporal and spatial resolution are used to evaluate the precursor signals of CI occurrence conditions, thereby providing strong support for the fixed-time and fixed-point CI nowcasting.

[0003] CI algorithms based on radar data mainly use radar echoes to reflect the intensity and structural changes of convective clouds. Usually, the occurrence and location of CI are determined by setting the radar reflectivity factor threshold ≥35dBz. Representative thunderstorm identification algorithms based on radar data technology abroad include radar echo correlation TREC algorithm, storm unit identification and tracking SCIT algorithm, thunderstorm nowcasting TSPN algorithm, etc. Short-term convective severe weather warning systems have also been developed in China, such as the storm nowcasting SWIFT system of the Guangdong Provincial Meteorological Bureau, the heavy rain nowcasting small vortex SWIRLS system of the Hong Kong Observatory, China, and the short-term nowcasting warning service SWAN system of the China Meteorological Administration.

[0004] CI algorithms based on satellite data mainly use satellite brightness temperature to reflect the temperature and height changes of convective cloud clusters, set thresholds and use machine learning methods to determine CI. Abroad, Maddox first proposed the brightness temperature and area threshold standards for mesoscale convective complexes based on satellite data. After that, representative algorithms based on satellite data technology also emerged, such as the nowcasting RDT algorithm for thunderstorm systems, the warning decision support WDSS-Ⅱ system, the satellite convection analysis and tracking SATCAST method, the convective incipient nowcasting UWCI method, and the RDCA algorithm for cumulus areas.

[0005] However, the above CI algorithms based on radar or satellite still have some problems. For example: (1) There are blind spots, insufficient coverage, and terrain influence in the radar scanning angle, which result in some convective clouds not being effectively captured or misjudged; (2) In addition, considering the influence of factors such as different regions, seasons, and underlying surfaces on CI characteristics, the setting and calculation methods of CI thresholds have certain differences and are difficult to unify; (3) Satellite brightness temperature is affected by factors such as cloud top phase, cloud microphysics, and cloud coverage, which result in some convective clouds not being effectively identified or misjudged; (4) There are certain differences between different satellite data, and the setting of CI thresholds is inconsistent; (5) Convection identification based on radar reflectivity factor threshold method is applied in an existing precipitation environment, resulting in a short advance forecast time, and there is still room for improvement in timeliness. Therefore, if CI identification based on radar or satellite is relied on alone, there will be a high false alarm rate, and the CI training samples constructed in this way will lead to low accuracy of CI forecast results based on machine learning methods. Summary of the invention

[0006] Purpose of the invention: The purpose of the present invention is to identify and forecast the initiation of convection in advance. It provides a satellite-ground collaborative convection initiation identification and forecasting method based on Fengyun satellites and integrating multi-source ground observation data. It can solve the problems of underreporting and false reporting in existing convection initiation forecasting methods to a certain extent.

[0007] Technical solution: A method for identifying and predicting the initial formation of satellite-ground coordinated convection based on physical identification of the present invention comprises the following steps:

[0008] Step 1: Collect and obtain satellite data in the target area, including geostationary meteorological satellite observation data, ground-based weather radar mosaic data, minute-level rainfall data from ground-based rain gauges, and high-resolution model numerical forecast products;

[0009] Step 2: Use the single-channel threshold method to extract convective cloud cluster information from satellite data. Based on the area overlap method, use the dynamic area threshold to determine the overlap rate, determine the tracking targets in turn, and carry out a series of tracking and numbering of convective cloud clusters. By calculating the continuous cooling rate of convective cloud clusters with the same number, determine the convective primary cloud clusters, and obtain the satellite-based convective primary identification results.

[0010] Step 3: read the radar echo data at a certain moment in the target area, obtain the cloud clusters in the grid connected domain whose radar echo is greater than the set threshold, and traverse to determine whether the cloud clusters meet the convection initiation conditions. If so, they are identified as convection initiation clouds; if not, they are not identified as convection initiation clouds;

[0011] Step 4: Record the cloud clusters identified as the initial convection, use the area overlap method and the cross-correlation method based on the sliding window method to track the recorded initial convection cloud clusters, and determine whether the tracked cloud clusters meet the conditions for continued tracking. If so, continue tracking in the next time period; if not, stop tracking;

[0012] Step 5: traverse all records, determine whether each record meets the conditions for retaining the record, and after determining the records to be retained, complete the radar-based convection initiation identification and tracking algorithm to obtain the radar-based convection initiation identification result;

[0013] Step 6: extract the radar-based and satellite-based convection primary recognition results and the corresponding minute-level rainfall data in the target area, and perform spatiotemporal matching. If the centroid latitude and longitude point of the satellite-based convection primary recognition result is within a preset range and before and after a preset time, the radar-based convection primary recognition result is marked as 1, otherwise it is marked as 0 and stored in the convection primary matching preprocessing database as the condition one result;

[0014] Step 7: determine whether there are stations with minute-level precipitation greater than a preset threshold within the movement range of the primary cloud cluster at the time of the satellite-based primary convection identification result triggering and within a certain preset time thereafter. If there is, mark it as 1, otherwise mark it as 0, and store it in the convection primary matching preprocessing database as the result of condition 2;

[0015] Step 8: Read the discrimination results of the convection initiation matching preprocessing database. Any event that meets one of the conditions 1 and 2 is identified as a convection initiation event. If not, the satellite-based convection initiation identification results are eliminated, thereby forming an accurate convection initiation data set sample.

[0016] Step 9: Based on the spatiotemporal information of convection initiation in the convection initiation dataset sample, using the regional numerical forecast model product, select or calculate and derive the physical quantities closely related to the occurrence and development of convection initiation within a certain preset range within the corresponding time, merge and classify them according to the season, perform feature analysis and extraction, convert them into a physical quantity matrix and input them into the LightGBM neural network model, so as to obtain the relevant physical quantities of each optimal weight, that is, the environmental parameter matrix based on physical judgment;

[0017] Step 10: Obtain the corresponding satellite parameters according to the convection incipient annotation data set, integrate the environmental parameter matrix based on physical identification and ground precipitation data to form a training data sample set, input it into the XGBOOST machine learning model, and build a convection incipient intelligent forecast model through training and model tuning;

[0018] Step 11: Input the real-time satellite data channel factors or parameters, environmental physical quantities and precipitation data provided by the model as prediction factors into the convection incipient intelligent forecast model to realize the convection incipient forecast and warning.

[0019] Further, step 1 is specifically as follows: the geostationary meteorological satellite observation data is the full disk and regional observation data of the FY-4A satellite, with a spatial resolution of 4km*4km, including 0.64μm reflectivity, 3.75μm brightness temperature, 6.25μm brightness temperature, 7.1μm brightness temperature, 8.5μm brightness temperature, 10.8μm brightness temperature, 12.0μm brightness temperature, 13.5μm brightness temperature, 30min 10.8μm brightness temperature change rate, 10.8μm and 6.2μm brightness temperature difference, 12μm and 10.8μm brightness temperature difference, 8.5μm and 10.8μm brightness temperature difference, 10.8μm and 6.2μm brightness temperature change rate, 12μm and 10 .8μm brightness temperature difference change rate, 8.5μm and 10.8μm brightness temperature difference change rate; the ground-based weather radar combined reflectivity is the 10-minute radar area puzzle data with a resolution of 0.01°×0.01°; the minute-by-minute rainfall data of the ground-based rain gauge, including the data of meteorological observation stations and regional automatic observation stations after quality control; the spatial resolution of the numerical model forecast product is 3km, and the selected physical quantities include 850hPa temperature, CAPE value, 850hPa divergence, 200hPa divergence, K index, 500hPa relative humidity, 850hPa relative humidity, convective inhibition energy CIN, and 0-6km vertical wind shear.

[0020] Furthermore, step 2 specifically includes the following steps:

[0021] Step 2.1, matching the full disk and target area satellite data in time to form a continuously collected satellite observation data sequence;

[0022] Step 2.2: For each satellite data at each time point, use the brightness temperature value of the 10.8 μm channel to retain the area below the set threshold to determine the convective cloud range;

[0023] Step 2.3: Use the connected domain segmentation algorithm to preliminarily extract the convective cloud cluster information;

[0024] Step 2.4, detect the area of ​​convective clouds, remove clouds with an area smaller than the minimum area threshold and larger than the maximum area threshold, and number the clouds;

[0025] Step 2.5, calculate the area overlap rate of the numbered cloud clusters before and after. For convective cloud clusters of different scales, set the dynamic area overlap rate threshold. When the threshold is reached, it is determined whether the next cloud cluster is developed from the previous cloud cluster. Match the two cloud clusters before and after, and follow the same steps to match the next time, thus forming the time series data of the same cloud cluster. The calculation formula of the overlap rate is as follows:

[0026]

[0027] Among them C i , C j is the cloud cluster at the previous and next time points, i, j are the cloud cluster numbers, and S is the cloud cluster area;

[0028] Step 2.6: Use the time series data of the cloud clusters with the same number to calculate the continuous cloud top cooling rate of the cloud clusters that last for more than 30 minutes. Mark the cloud clusters that reach the cooling rate threshold as nascent clouds and continue to track them.

[0029] Furthermore, step 3 specifically includes the following steps:

[0030] Step 3.1, read the radar combined reflectivity data at the corresponding time, set the grid points less than 30 dBz in the radar echo data to 0, and set the grid points greater than or equal to 30 dBz to 1;

[0031] Step 3.2, specify the eight grid points above, below, left, right and diagonally of the grid point as adjacent grid points, and use the eight-neighborhood connected domain algorithm to obtain all connected domains in the radar echo data;

[0032] Step 3.3, traverse all connected domains. If the grid point where the connected domain is located is greater than 30 dBz, then the connected domain is considered to be a cloud cluster greater than 30 dBz. Otherwise, the connected domain is not considered to be a cloud cluster greater than 30 dBz. After the traversal is completed, the cloud cluster position of the grid point connected domain where the radar echo is greater than 30 dBz can be obtained;

[0033] Step 3.4: Determine whether the cloud cluster meets the conditions for the initiation of convection. If the cloud cluster meets all the above conditions, it is determined to be a convective initiation cloud cluster; if the cloud cluster does not meet any of the above conditions, it is determined not to be a convective initiation cloud cluster.

[0034] Furthermore, in step 3.4, the convection initiation condition is specifically: the cloud cluster area is 60km 2 Up to 1000km 2 There are 8 or more grid points exceeding 35dBz in the cloud cluster; the number of echo grid points above 35dBz in the cloud cluster is 1.45 times less than the number of echo grid points above 35dBz within 15km of the location of each cloud cluster at the previous time; the central longitude and latitude of the cloud cluster are outside the range of the cloud cluster tracked at the previous time.

[0035] Furthermore, step 4 specifically includes the following steps:

[0036] Step 4.1, read the radar echo data at the next moment, and obtain the cloud clusters with a value greater than 30 dBz in the radar echo data at the next moment;

[0037] Step 4.2, determine the cloud group to be tracked according to the records, and obtain the longitude and latitude of the cloud group to be tracked and the cloud group with a signal greater than 30 dBz in the radar echo data at the next moment;

[0038] Step 4.3, determine the overlap between the cloud to be tracked and each cloud with a position greater than 30 dBz in the radar echo data at the next moment. If the overlap between a certain cloud and the cloud to be tracked is greater than 25%, then determine that the cloud is the tracking cloud of the cloud to be tracked at the next moment. If there are multiple cloud clusters with a position overlap of greater than 25% with the cloud to be tracked, then select the cloud with the largest position overlap as the tracking cloud of the cloud to be tracked.

[0039] Step 4.4, if the overlap between the position of no cloud and the cloud to be tracked is greater than 25%, then calculate the average radius of the cloud to be tracked, obtain the grid points with a radius of 25 km around the average radius of the cloud to be tracked in the radar echo data at the next moment, slide and extract the data window with the size of the average radius of the cloud to be tracked in these grid points, and calculate the correlation between the data window and the data within the average radius of the cloud to be tracked. The operation formula is:

[0040]

[0041] Among them, I t (i, j) and I t+Δt (i, j) represents the pixel value of the image at position (i, j) at time t and t+△t respectively; μ t and μ t+△t They are image I t and I t+△t The mean of ; n and m are the number of rows and columns of the image respectively;

[0042] If there is no grid point with correlation greater than 0.35 in the sliding window, the cloud is judged to be dissipated and no longer tracked. Otherwise, the grid point with the highest dBz in the data window with the highest correlation in the sliding window is selected, and the grid point with a correlation greater than 30 dBz is obtained through the connected domain as the tracking cloud cluster to track the cloud cluster;

[0043] Step 4.5: Determine whether the tracking cloud satisfies the conditions for continuing tracking. If the cloud satisfies all the above conditions, it is determined to continue tracking; if the tracking cloud does not meet any of the above conditions, the tracking of the cloud is stopped.

[0044] Furthermore, in step 4.5, the conditions for continuing tracking are as follows: the cloud area to be tracked is within 30km 2 Up to 2000km 2 There are grid points with a brightness exceeding 35 dBz in the tracked cloud cluster; the tracked cloud cluster does not merge into other cloud clusters within 30 minutes; the central longitude and latitude of the cloud cluster are outside the range of the cloud cluster tracked at the last time.

[0045] Furthermore, in step 5, it is determined whether each record meets the conditions for retaining the record, and the specific determination conditions are: the cloud cluster duration is less than 30 minutes; the cloud cluster merges into other cloud clusters within 30 minutes; two newly formed cloud clusters are close and the cloud cluster is small; the newly formed cloud cluster is close to the cloud cluster being tracked; the maximum radar echo of the cloud cluster during the entire tracking process is less than 40dBz; if the cloud cluster meets any of the above conditions, the record of the cloud cluster is deleted.

[0046] Furthermore, step 9 specifically includes the following steps:

[0047] Step 9.1, obtain the time and longitude and latitude position of CI in the sample of the convection primary data set, and extract the physical quantity data of 500hPa temperature, 850hPa temperature, 700hPa temperature dew point difference, CAPE value, 850hPa divergence, 200hPa divergence, K index, 500hPa relative humidity, 850hPa relative humidity, convection inhibition energy CIN, and 0-6km vertical wind shear in the regional numerical forecast model products within 32km in the corresponding time. Among them, the physical quantities of K index, divergence, and vertical wind shear are calculated by formulas. The specific calculation formulas are as follows:

[0048] K Index:

[0049] The K index is a meteorological unit used to analyze the stability of the atmosphere. The first term in the K index calculation formula represents the temperature lapse rate, the second term represents the low-level water vapor conditions, and the third term represents the saturation level in the middle layer. The larger the K index, the more unstable the stratification. Its specific expression is:

[0050] K=(T 850 -T 500 )+T d850 -(TT d ) 700

[0051] Where T 850 , T 700 and T 500 Represents the temperature at 850hPa, 700hPa and 500hPa respectively, T d700 and T d850 Represents the dew point temperature at 700hPa and 850hPa respectively;

[0052] Divergence:

[0053] Divergence indicates the rate of change of unit volume during fluid movement. The area where the fluid is concentrated during movement is called convergence, and the area where the fluid is dispersed during movement is called divergence. When the divergence value is negative, it indicates convergence, which is conducive to the development and enhancement of convective weather systems. When it is positive, it indicates divergence, which is not conducive to the development of weather systems.

[0054]

[0055] Where u is the wind speed in the x-axis direction, v is the wind speed in the y-axis direction, the x-axis represents the tangent direction of the latitude circle, where the positive direction is from west to east, and the y-axis represents the tangent direction of the meridian circle, where the positive direction is from south to north;

[0056] 0-6km vertical wind shear:

[0057] Vertical wind shear refers to the rate of change of wind speed in the vertical direction of the atmosphere, which is calculated by the difference in wind speed and wind direction components at two different altitudes;

[0058]

[0059] Where u6 and u0 represent the zonal wind speed at 6km and 0km respectively; v6 and v0 represent the meridional wind speed at 6km and 0km respectively;

[0060] Then, the data integrity is checked and supplemented, and the processed data are merged and classified according to season;

[0061] Step 9.2: Perform feature processing on the physical quantities before inputting the model:

[0062] Firstly, based on CI events in different seasons, the meteorological analysis method was used to screen out environmental physical parameters and their threshold ranges that were highly correlated with CI.

[0063] Then, the variance selection method is used to calculate the variance of highly correlated physical quantities. The features with larger variance thresholds are selected through the threshold, and preliminary weights are assigned to various features according to the calculation results to avoid excessive similarity of features, resulting in too low feature values. The variance calculation formula is as follows:

[0064]

[0065] Among them, x i is each observation in the sample, is the sample mean, n is the number of observations in the sample;

[0066] Finally, through the method of correlation coefficient calculation and chi-square test, the correlation coefficient value of each feature to the target value is counted, and then the secondary weight is assigned according to the strength of the correlation, combined with the divergent weight, the final weight value of each feature is integrated; Pearson correlation coefficient:

[0067]

[0068] Among them, x i and i represent the observed values ​​of two variables X and Y respectively, and denote the means of X and Y respectively, and n is the number of observations;

[0069] Spearman correlation coefficient:

[0070]

[0071] Among them, x i and i represent the observed values ​​of two variables X and Y respectively, and represent the means of X and Y respectively;

[0072] Step 9.3: Input the physical quantities after feature processing into the LightGBM neural network model, perform feature analysis and extraction on the data, and then obtain the optimal weights of each related physical quantity, that is, the environmental parameter matrix based on physical judgment.

[0073] Furthermore, step 10 specifically includes the following steps:

[0074] Step 10.1, obtain the corresponding satellite parameters according to the convection primary annotation data set, integrate the environmental parameter matrix based on physical identification and ground precipitation data, and form a sample set for building the CI intelligent forecast model. The satellite parameters specifically include:

[0075] 0.64μm albedo, used to represent effective particle radius or optical thickness;

[0076] 3.75 μm, used to represent high and low albedo targets, respectively;

[0077] 6.25 μm, used to represent high-level water vapor;

[0078] 7.1 μm, used to represent water vapor in the middle layer;

[0079] 8.5 μm, used to represent total water vapor and clouds;

[0080] 10.8μm, used to indicate cloud top height;

[0081] 12.0 μm, used to represent clouds and total water vapor;

[0082] 13.5 μm, used to represent small particle aerosols;

[0083] 10.8-6.2 μm brightness temperature difference, used to indicate the height of the cloud top relative to the upper troposphere;

[0084] 10.8-12.0μm brightness temperature difference, used to represent the split window cloud top phase;

[0085] 8.5-10.8μm brightness temperature difference, used to indicate the cloud top phase;

[0086] 2*10.8-12.0-8.5μm brightness temperature difference, used to indicate the cloud top phase;

[0087] 10.8μm brightness temperature temporal change rate, used to indicate cloud lift;

[0088] The time rate of change of the brightness temperature difference between 10.8 and 6.2 μm is used to indicate the lift of clouds;

[0089] The time change rate of the brightness temperature difference between 10.8 and 12.0 μm is used to indicate the phase change of the cloud top;

[0090] The time change rate of the brightness temperature difference between 8.5 and 10.8 μm is used to indicate the change of cloud top freezing;

[0091] Step 10.2, split the CI sample set formed in step 10.1 into a training set and a test set at a ratio of 4:1, and resample the training samples using the SMOTE method;

[0092] Step 10.3: Substitute the XGBOOST model for training, use MSE as the loss function, adjust the learning rate parameters, use the gradient descent method for multiple rounds of training, compare and test the results, obtain the optimal XGBOOST model and save it to form a convection incipient intelligent forecasting model.

[0093] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:

[0094] (1) The satellite-based convection initiation identification algorithm in the present invention adopts a dynamic area threshold to determine the overlap rate based on the traditional area overlap method. That is, different overlap thresholds are used according to the size of the front and rear cloud clusters to determine the tracking cloud cluster target in turn. This algorithm can comprehensively consider the situations of large clouds merging into small clouds, small clouds merging into large clouds, and large clouds splitting, and effectively improve the CI identification effect based on satellite data.

[0095] (2) The prediction accuracy of the artificial intelligence forecasting model is more dependent on the quality of the training samples. Based on the above more accurate satellite CI recognition results, the present invention further performs high-temporal and spatial matching with radar-based CI recognition results, minute-level rainfall data, etc., to form a more accurate and high-quality convection primary annotation dataset. This dataset is used as a training sample for the CI intelligent forecasting model, which can reduce noise and bias in model training and better improve the prediction performance of the model.

[0096] (3) Based on the spatiotemporal information of CI in the samples of the primary convection data set, the present invention uses regional numerical forecast model products to select or calculate and derive physical quantities closely related to the occurrence and development of CI within the corresponding spatiotemporal range, merge and classify them according to the season, perform feature analysis and extraction, and convert them into a physical quantity matrix for input into the LightGBM neural network model to form an optimal combination of physical quantities of the CI environment based on physical judgment. Combined with high-quality primary convection labeled data sets, model precipitation, etc., this can provide high-quality forecast factors for the CI intelligent forecast model. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1 It is a schematic diagram of the process of the present invention;

[0098] Figure 2 It is an objective CI identification result based on satellite data;

[0099] Figure 3 It is the objective CI identification result based on radar data;

[0100] Figure 4 It is a process of matching satellite-based CI identification results with radar CI identification results, precipitation, etc. to form a convection primary annotation data set;

[0101] Figure 5 This is a forecast case of the convective inception intelligent forecast model based on physical identification, where (a) is a CI event identified by the satellite on August 2, 2023, and (b) is the radar echo and radar CI identification algorithm result of the corresponding event. DETAILED DESCRIPTION

[0102] The technical solution of the present invention is further described below in conjunction with the accompanying drawings.

[0103] like Figure 1 As shown, the present invention proposes a method for early warning of convective initiation based on multi-source data and physical identification to construct training samples. This method mainly matches the satellite identification results of convective initiation with the radar identification results of convective initiation and process precipitation to obtain a high-quality convective initiation data set, and combines the convective initiation physical identification results to construct a model training data set, and then intelligently predicts the convective initiation through XGBOOST machine learning modeling.

[0104] 1) Data preprocessing: preprocess the 4km resolution L1 data of the geostationary orbit radiation imager of Fengyun-4A, including radiation calibration, projection conversion, channel separation, etc. The full disk and China regional satellite data are fused in time to form a continuously sampled data sequence;

[0105] 2) Thin clouds and mature clouds are removed by using the threshold set by the brightness temperature of the 10.8 μm channel, and convective cloud targets that may develop are retained. By calculating the overlap rate of the convective cloud clusters in the two satellite cloud images, a dynamic area overlap threshold is used to determine whether the cloud clusters in the previous and next sequences are the same cloud cluster evolution, and number them;

[0106] 3) Calculate the continuous cooling rate of the cloud cluster with the same number in the continuous sequence of satellite cloud images, mark the convective cloud cluster with a cooling rate of 4K / 15min and a development duration of more than 30min as convective incipient, and continue to track it until the cloud cluster merges or disappears;

[0107] 4) Through the above steps, a CI tag data sequence based on satellite data is formed ( Figure 2 );

[0108] 5) Read radar echo data, obtain cloud clusters with grid connection domains greater than 30 dBz, use area overlap method and cross-correlation method based on sliding window method to track recorded convective primary cloud clusters, and judge whether the tracked cloud clusters meet the conditions for retaining records, complete radar-based convective primary identification and tracking, and obtain radar-based CI identification results ( Figure 3 );

[0109] 6) Perform spatiotemporal matching on the satellite-based CI tag data sequence and the radar-based CI recognition result. If there is a radar-based CI mark within 20km of the longitude and latitude of the satellite-based CI centroid and within 30 minutes before and after, the satellite-based CI tag is confirmed to be accurate, and the area is marked as 1, otherwise it is marked as 0, and the matching result is stored in the CI matching preprocessing database as condition 1;

[0110] 7) Determine the stations with precipitation ≥ 0.1 mm within the CI cloud movement range based on the satellite CI triggering time and within 60 minutes thereafter. If there is, mark it as 1, otherwise mark it as 0, and store it in the CI matching preprocessing database as the result of condition 2;

[0111] 8) Based on steps 6 and 7, all events that meet condition 1 or condition 2 are identified as CI events. If they do not meet the conditions, the satellite CI results are removed to form a labeled convection primary label dataset. Figure 4 ;

[0112] 9) Using regional numerical forecast model products, derive physical quantities closely related to the occurrence and development of CI within a range of 32 km within 1 hour before CI occurs, perform feature analysis and extraction, and generate the optimal weight combination of physical quantities with high correlation characteristics covering a certain spatiotemporal window of CI occurrence (i.e., the environmental parameter matrix based on physical identification), match it with CI label data, ground observation precipitation data, etc., to form a training data sample set (Table 1), input it into the XGBOOST machine learning model, and build a convection incipient intelligent forecast model through training and model tuning;

[0113] 10) The real-time satellite data channel factors, environmental physical quantities and precipitation data provided by the model are used as prediction factors and input into the convection incipient intelligent forecast model to achieve CI forecast and warning. Take a CI event on August 2, 2023 as an example ( Figure 5 ), the convective initial intelligent forecast model issued a CI warning in southwestern Hubei at 11:38 Beijing time, Figure 5 (a) The red label area, at 11:42 Beijing time, a radar echo greater than 20dBz appeared at this location Figure 5 In the black box area (b), at 11:48, the radar recognition algorithm identified this echo as the incipient convection.

[0114] Table 1 Environmental physical quantity factors

[0115]

[0116]

Claims

1. A method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification, characterized in that: The steps include: Step 1: Collect and obtain satellite data in the target area, including geostationary meteorological satellite observation data, ground-based weather radar mosaic data, minute-level rainfall data from ground-based rain gauges, and high-resolution model numerical forecast products; Step 2: Use the single-channel threshold method to extract convective cloud cluster information from satellite data. Based on the area overlap method, use the dynamic area threshold to determine the overlap rate, determine the tracking targets in turn, carry out a series of tracking and numbering of convective cloud clusters, and determine the convective primary cloud clusters by calculating the continuous cooling rate of the convective cloud clusters with the same number, and obtain the satellite-based convective primary CI identification results. Step 3: read the radar echo data at a certain moment in the target area, obtain the cloud clusters whose radar echoes are greater than the set threshold grid connected domain, and traverse to determine whether the cloud clusters meet the convection initiation conditions. If so, they are identified as convection initiation clouds; if not, they are not identified as convection initiation clouds; Step 4: Record the cloud clusters identified as the initial convection, use the area overlap method and the cross-correlation method based on the sliding window method to track the recorded initial convection cloud clusters, and determine whether the tracked cloud clusters meet the conditions for continued tracking. If so, continue tracking in the next time period; if not, stop tracking; Step 5: traverse all records, determine whether each record meets the conditions for retaining the record, and after determining the records to be retained, complete the radar-based convection initiation identification and tracking algorithm to obtain the radar-based convection initiation identification result; Step 6: extract the radar-based and satellite-based convection primary recognition results and the corresponding minute-level rainfall data in the target area, and perform spatiotemporal matching. If the centroid longitude and latitude point of the satellite-based convection primary recognition result is within a preset range and before and after a preset time, the radar-based convection primary recognition result is marked as 1, otherwise it is marked as 0 and stored in the CI matching preprocessing database as the condition one result; Step 7: Determine whether there are stations with minute-level precipitation greater than a preset threshold within the movement range of the CI cloud group at the time of triggering the satellite-based convection primary identification result and within a certain preset time thereafter. If there is, mark it as 1, otherwise mark it as 0, and store it in the CI matching preprocessing database as the result of condition 2; Step 8: Read the discrimination results of the CI matching preprocessing database. All events that meet one of the conditions 1 and 2 are identified as CI events. If not, the satellite-based convective primary CI identification results are eliminated, thereby forming an accurate convective primary data set sample. Step 9: Based on the spatiotemporal information of CI in the sample of the convection primary data set, the regional numerical forecast model products are used to select or calculate and derive the physical quantities closely related to the occurrence and development of CI within a certain preset range within the corresponding time, merge and classify them according to the season, and perform feature analysis and extraction, and convert them into a physical quantity matrix to input into the LightGBM neural network model, so as to obtain the relevant physical quantities of each optimal weight, that is, the environmental parameter matrix based on physical judgment; Step 10: Obtain the corresponding satellite parameters according to the convection incipient annotation data set, integrate the environmental parameter matrix based on physical identification and ground precipitation data to form a training data sample set, input it into the XGBOOST machine learning model, and build a convection incipient intelligent forecast model through training and model tuning; Step 11: Input the real-time satellite data channel factors or parameters, environmental physical quantities and precipitation data provided by the model as prediction factors into the convection incipient intelligent forecast model to realize the convection incipient forecast and warning.

2. According to the physical identification-based method for identifying and predicting the incipient satellite-ground coordinated convection according to claim 1, it is characterized in that: Step 1 is specifically as follows: the geostationary meteorological satellite observation data is the full disk and regional observation data of the FY-4A satellite, with a spatial resolution of 4km*4km, including 0.64μm reflectivity, 3.75μm brightness temperature, 6.25μm brightness temperature, 7.1μm brightness temperature, 8.5μm brightness temperature, 10.8μm brightness temperature, 12.0μm brightness temperature, 13.5μm brightness temperature, 30min 10.8μm brightness temperature change rate, 10.8μm and 6.2μm brightness temperature difference, 12μm and 10.8μm brightness temperature difference, 8.5μm and 10.8μm brightness temperature difference, 10.8μm and 6.2μm brightness temperature difference change rate, 12μm and 10.8μm brightness temperature difference, The brightness temperature difference change rate of 8.5μm and 10.8μm respectively; the combined reflectivity of the ground-based weather radar is the radar area mosaic data for 10 minutes with a resolution of 0.01°×0.01°; the minute-by-minute rainfall data of the ground-based rain gauge includes the data of meteorological observation stations and regional meteorological automatic observation stations after quality control; the spatial resolution of the numerical model forecast product is 3km, and the selected physical quantities include 850hPa temperature, CAPE value, 850hPa divergence, 200hPa divergence, K index, 500hPa relative humidity, 850hPa relative humidity, convective inhibition energy CIN, and 0-6km vertical wind shear.

3. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 1 is characterized in that: Step 2 specifically includes the following steps: Step 2.1, matching the full disk and target area satellite data in time to form a continuously collected satellite observation data sequence; Step 2.2: For each satellite data at each time point, use the brightness temperature value of the 10.8 μm channel to retain the area below the set threshold to determine the convective cloud range; Step 2.3: Use the connected domain segmentation algorithm to preliminarily extract the convective cloud cluster information; Step 2.4, detect the area of ​​convective clouds, remove clouds with an area smaller than the minimum area threshold and larger than the maximum area threshold, and number the clouds; Step 2.5, calculate the area overlap rate of the numbered cloud clusters before and after. For convective cloud clusters of different scales, set the dynamic area overlap rate threshold. When the threshold is reached, it is determined whether the next cloud cluster is developed from the previous cloud cluster. Match the two cloud clusters before and after, and follow the same steps to match the next time, thus forming the time series data of the same cloud cluster. The calculation formula of the overlap rate is as follows: Among them C i , C j is the cloud cluster at the previous and next time points, i, j are the cloud cluster numbers, and S is the cloud cluster area; Step 2.6: Use the time series data of the cloud clusters with the same number to calculate the continuous cloud top cooling rate of the cloud clusters that last for more than 30 minutes. Mark the cloud clusters that reach the cooling rate threshold as nascent clouds and continue to track them.

4. The method for identifying and predicting the initial formation of satellite-ground coordinated convection based on physical identification according to claim 1 is characterized in that: Step 3 specifically includes the following steps: Step 3.1, read the radar combined reflectivity data at the corresponding time, set the grid points less than 30 dBz in the radar echo data to 0, and set the grid points greater than or equal to 30 dBz to 1; Step 3.2, specify the eight grid points above, below, left, right and diagonally of the grid point as adjacent grid points, and use the eight-neighborhood connected domain algorithm to obtain all connected domains in the radar echo data; Step 3.3, traverse all connected domains. If the grid point where the connected domain is located is greater than 30 dBz, then the connected domain is considered to be a cloud cluster greater than 30 dBz. Otherwise, the connected domain is not considered to be a cloud cluster greater than 30 dBz. After the traversal is completed, the cloud cluster position of the grid point connected domain where the radar echo is greater than 30 dBz can be obtained; Step 3.4: Determine whether the cloud cluster meets the conditions for the initiation of convection. If the cloud cluster meets all the above conditions, it is determined to be a convective initiation cloud cluster; if the cloud cluster does not meet any of the above conditions, it is determined not to be a convective initiation cloud cluster.

5. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 4 is characterized in that: In step 3.4, the initial convection conditions are as follows: the cloud cluster area is 60 km 2 Up to 1000km 2 There are 8 or more grid points exceeding 35dBz in the cloud cluster; the number of echo grid points above 35dBz in the cloud cluster is 1.45 times less than the number of echo grid points above 35dBz within 15km of the location of each cloud cluster at the previous time; the central longitude and latitude of the cloud cluster are outside the range of the cloud cluster tracked at the previous time.

6. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 1 is characterized in that: Step 4 specifically includes the following steps: Step 4.1, read the radar echo data at the next moment, and obtain the cloud clusters with a value greater than 30 dBz in the radar echo data at the next moment; Step 4.2, determine the cloud group to be tracked according to the records, and obtain the longitude and latitude of the cloud group to be tracked and the cloud group with a signal greater than 30 dBz in the radar echo data at the next moment; Step 4.3, determine the overlap between the cloud to be tracked and each cloud with a position greater than 30 dBz in the radar echo data at the next moment. If the overlap between a certain cloud and the cloud to be tracked is greater than 25%, then determine that the cloud is the tracking cloud of the cloud to be tracked at the next moment. If there are multiple cloud clusters with a position overlap of greater than 25% with the cloud to be tracked, then select the cloud with the largest position overlap as the tracking cloud of the cloud to be tracked. Step 4.4, if the overlap between the position of no cloud and the cloud to be tracked is greater than 25%, then calculate the average radius of the cloud to be tracked, obtain the grid points with a radius of 25 km around the average radius of the cloud to be tracked in the radar echo data at the next moment, slide and extract the data window with the size of the average radius of the cloud to be tracked in these grid points, and calculate the correlation between the data window and the data within the average radius of the cloud to be tracked. The operation formula is: Among them I t (i, j) and I t+Δt (i, j) represents the pixel value of the image at position (i, j) at time t and t+△t respectively; μ t and μ t+△t They are image I t and I t+△t The mean of ; n and m are the number of rows and columns of the image respectively; If there is no grid point with correlation greater than 0.35 in the sliding window, the cloud is judged to be dissipated and no longer tracked. Otherwise, the grid point with the highest dBz in the data window with the highest correlation in the sliding window is selected, and the grid point with a correlation greater than 30 dBz is obtained through the connected domain as the tracking cloud cluster to track the cloud cluster; Step 4.5: Determine whether the tracking cloud satisfies the conditions for continuing tracking. If the cloud satisfies all the above conditions, it is determined to continue tracking; if the tracking cloud does not meet any of the above conditions, the tracking of the cloud is stopped.

7. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 6 is characterized in that: In step 4.5, the conditions for continuing tracking are: the cloud area is within 30km 2 Up to 2000km 2 There are grid points with a brightness exceeding 35 dBz in the tracked cloud cluster; the tracked cloud cluster does not merge into other cloud clusters within 30 minutes; the central longitude and latitude of the cloud cluster are outside the range of the cloud cluster tracked at the last time.

8. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 1 is characterized in that: In step 5, it is determined whether each record meets the conditions for retaining the record, and the specific determination conditions are: the cloud cluster duration is less than 30 minutes; the cloud cluster merges into other cloud clusters within 30 minutes; two newly formed cloud clusters are close and the cloud cluster is small; the newly formed cloud cluster is close to the cloud cluster being tracked; the maximum radar echo of the cloud cluster during the entire tracking process is less than 40dBz; if the cloud cluster meets any of the above conditions, the record of the cloud cluster is deleted.

9. The method for identifying and predicting the initial formation of satellite-ground coordinated convection based on physical identification according to claim 1 is characterized in that: Step 9 specifically includes the following steps: Step 9.1, obtain the time and longitude and latitude position of CI in the sample of the convection primary data set, and extract the physical quantity data of 500hPa temperature, 850hPa temperature, 700hPa temperature dew point difference, CAPE value, 850hPa divergence, 200hPa divergence, K index, 500hPa relative humidity, 850hPa relative humidity, convection inhibition energy CIN, and 0-6km vertical wind shear in the regional numerical forecast model products within 32km in the corresponding time. Among them, the physical quantities of K index, divergence, and vertical wind shear are calculated by formulas. The specific calculation formulas are as follows: K Index: The K index is a meteorological unit used to analyze the stability of the atmosphere. The first term in the K index calculation formula represents the temperature lapse rate, the second term represents the low-level water vapor conditions, and the third term represents the saturation level in the middle layer. The larger the K index, the more unstable the stratification. Its specific expression is: K=(T 850 -T 500 )+T d850 -(T-T d ) 700 Where T 850 , T 700 and T 500 Represents the temperature at 850hPa, 700hPa and 500hPa respectively, T d700 and T d850 Represents the dew point temperature at 700hPa and 850hPa respectively; Divergence: Divergence indicates the rate of change of unit volume during fluid movement. The area where the fluid is concentrated during movement is called convergence, and the area where the fluid is dispersed during movement is called divergence. When the divergence value is negative, it indicates convergence, which is conducive to the development and enhancement of convective weather systems. When it is positive, it indicates divergence, which is not conducive to the development of weather systems. Where u is the wind speed in the x-axis direction, v is the wind speed in the y-axis direction, the x-axis represents the tangent direction of the latitude circle, where the positive direction is from west to east, and the y-axis represents the tangent direction of the meridian circle, where the positive direction is from south to north; 0-6km vertical wind shear: Vertical wind shear refers to the rate of change of wind speed in the vertical direction of the atmosphere, which is calculated by the difference in wind speed and wind direction components at two different altitudes; Where u6 and u0 represent the zonal wind speed at 6km and 0km respectively; v6 and v0 represent the meridional wind speed at 6km and 0km respectively; Then the data integrity is checked and supplemented, and the processed data are merged and classified according to season; Step 9.2: Perform feature processing on the physical quantities before inputting the model: Firstly, based on CI events in different seasons, the meteorological analysis method was used to screen out environmental physical parameters and their threshold ranges that were highly correlated with CI. Then, the variance selection method is used to calculate the variance of highly correlated physical quantities. The features with larger variance thresholds are selected through the threshold, and preliminary weights are assigned to various features according to the calculation results to avoid excessive similarity of features, resulting in too low feature values. The variance calculation formula is as follows: Where xi is each observation in the sample, is the sample mean, n is the number of observations in the sample; Finally, through the method of correlation coefficient calculation and chi-square test, the correlation coefficient value of each feature to the target value is counted, and then the secondary weight is assigned according to the strength of the correlation, combined with the divergent weight, to integrate the final weight value of each feature; Pearson correlation coefficient: Among them, x i and i represent the observed values ​​of two variables X and Y respectively, and denote the means of X and Y respectively, and n is the number of observations; Spearman correlation coefficient: Among them, x i and i represent the observed values ​​of two variables X and Y respectively, and represent the means of X and Y respectively; Step 9.3: Input the physical quantities after feature processing into the LightGBM neural network model, perform feature analysis and extraction on the data, and then obtain the optimal weights of each related physical quantity, that is, the environmental parameter matrix based on physical judgment.

10. The method for identifying and predicting the incipient satellite-ground coordinated convection based on physical identification according to claim 1, characterized in that: Step 10 specifically includes the following steps: Step 10.1, obtain the corresponding satellite parameters according to the convection primary annotation data set, integrate the environmental parameter matrix based on physical identification and ground precipitation data, and form a sample set for building the CI intelligent forecast model. The satellite parameters specifically include: 0.64μm albedo, used to represent effective particle radius or optical thickness; 3.75 μm, used to represent high and low albedo targets, respectively; 6.25 μm, used to represent high-level water vapor; 7.1 μm, used to represent water vapor in the middle layer; 8.5 μm, used to represent total water vapor and clouds; 10.8μm, used to indicate cloud top height; 12.0 μm, used to represent clouds and total water vapor; 13.5 μm, used to represent small particle aerosols; 10.8-6.2 μm brightness temperature difference, used to indicate the height of the cloud top relative to the upper troposphere; 10.8-12.0μm brightness temperature difference, used to represent the split window cloud top phase; 8.5-10.8μm brightness temperature difference, used to indicate the cloud top phase; 2*10.8-12.0-8.5μm brightness temperature difference, used to indicate the cloud top phase; 10.8μm brightness temperature temporal change rate, used to indicate cloud lift; The time rate of change of the brightness temperature difference between 10.8 and 6.2 μm is used to indicate the lift of clouds; The time change rate of the brightness temperature difference between 10.8 and 12.0 μm is used to indicate the phase change of the cloud top; The time change rate of the brightness temperature difference between 8.5 and 10.8 μm is used to indicate the change of cloud top freezing; Step 10.2, split the CI sample set formed in step 10.1 into a training set and a test set at a ratio of 4:1, and resample the training samples using the SMOTE method; Step 10.3: Substitute the XGBOOST model for training, use MSE as the loss function, adjust the learning rate parameters, use the gradient descent method for multiple rounds of training, compare and test the results, obtain the optimal XGBOOST model and save it to form a convection incipient intelligent forecasting model.

Images