Double-layer cloud wind direction retrieval algorithm based on double-layer cloud top height information

Abstract

The present application belongs to the technical field of satellite cloud wind retrieval, and specifically relates to a double-layer cloud wind retrieval algorithm based on double-layer cloud top height information. The present application utilizes the double-layer cloud (upper ice cloud and lower water cloud) cloud top height information retrieved from the thermal infrared multi-channel observation of Japanese HImawari-8 satellite, and according to the pixel offset of the double-layer cloud at continuous time, the optical flow field at the corresponding time is estimated by using the dense optical flow method, and then the cloud wind vector is converted. The present application firstly obtains the double-layer cloud wind information, and the cloud wind retrieval test for the 7th tropical storm "Mulan" in 2022 shows that the spiral center and wind field difference of the hollow structure typhoon at different height layers can be accurately retrieved, which has important significance for the typhoon path and intensity prediction and the assimilation application of cloud wind data in numerical prediction.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of satellite cloud wind retrieval, and particularly relates to a double-layer cloud wind retrieval algorithm based on double-layer cloud top height information. BACKGROUND

[0002] Numerical weather prediction is an important means of modern weather prediction, and its essence is to numerically solve the state of atmospheric motion on the basis of given initial value and boundary condition from the physical law inside the atmosphere. With the development of numerical prediction model and high-performance supercomputer, the physical process considered by the model is gradually improved, and the computing power and running efficiency have also achieved great leap. Therefore, the initial value problem is the most core problem in the current numerical weather prediction. Data assimilation can fuse the model background field and multi-source observation data, and provide more accurate initial field for numerical weather prediction, which is one of the important means to improve the prediction ability of the model.

[0003] Among numerous meteorological observation data, conventional meteorological observation is difficult to provide effective initial field information on small and medium scale prediction due to the limitation of spatial and temporal resolution. Satellite can provide more wide coverage and higher spatial and temporal resolution atmospheric observation, and is widely used in global and regional numerical prediction business systems. The cloud wind retrieved by satellite observation and retrieval can not only effectively supplement the wind field data in areas without observation stations such as plateau and ocean, but also can be applied in numerical prediction initial field to improve the prediction accuracy of the model. Cloud wind can accurately and clearly show the changes of weather system and surrounding environment field, and has important indicative significance for rainstorm area analysis and prediction, typhoon path and intensity prediction, etc.

[0004] As early as in the 1970s, American scholars began to manually calculate cloud wind by using visible light cloud image. After decades of development, the cloud wind retrieval algorithm has made great progress, and the cross-correlation coefficient method and the optical flow method based on continuous time cloud tracking have been proposed. Subsequent scholars of various countries have continuously improved the two kinds of retrieval algorithms by combining fast Fourier transform, image features and neural network algorithm. The improved maximum correlation coefficient method is simple to implement and has improved computing efficiency, but it cannot meet the correlation coefficient condition on some deformed, nascent and decaying cloud clusters, thus greatly limiting the accuracy and applicability of cloud wind. The optical flow method can infer the instantaneous speed of pixel motion by studying the change of image gray value at continuous time, and is widely used in the research fields of convection tracking and satellite calculation of atmospheric motion vector. Compared with the correlation coefficient method, the algorithm is not limited by cloud shape change, has high operation efficiency, and can obtain pixel-by-pixel and refined speed vector.

[0005] The cloud wind derived by the optical flow method can effectively improve the typhoon path prediction and intensity prediction in data assimilation. However, the cloud wind derived by the satellite can only reflect the movement of the uppermost cloud, and the wind field information of the lower layer cannot be obtained. From the actual application effect, the cloud top cloud wind available at present is limited in the improvement effect of the model, and the cloud wind derivation algorithm and assimilation application still need to be improved. With the increase of the spatial and spectral resolution of the satellite and the development of artificial intelligence technology, the satellite can derive the cloud top information of two layers of ice cloud and water cloud in the double-layer cloud, thereby providing great technical support for the double-layer cloud wind derivation. SUMMARY

[0006] In order to solve the single direction of the existing satellite cloud wind derivation, the double-layer cloud wind derivation algorithm based on the double-layer cloud top height information is provided by using the cloud top height information of the double-layer cloud (upper ice cloud and lower water cloud).

[0007] The double-layer cloud wind derivation algorithm based on the double-layer cloud top height information provided by the application has the following specific steps:

[0008] (1) Image preprocessing, which is to process the cloud top height images at different times into a unified gray value matrix according to the same threshold range. Specifically, the maximum cloud top height of the upper layer is fixed as 18 kilometers, and the maximum cloud top height of the lower double-layer water cloud is fixed as 8 kilometers; the minimum cloud top height is set as 0.5 kilometers. The pixels with a cloud height less than 0.5 kilometers are set as a default value, and the pixels with an upper cloud height greater than 18 kilometers (a lower cloud height greater than 8 kilometers) are all set as 18 kilometers. After this preprocessing, the numerical range of the upper cloud height is all between 0.5 and 18, and the numerical range of the lower cloud height is all between 0.5 and 8. In the upper and lower layers of the cloud, the cloud height is mapped to the gray value range of 0-255 through the maximum and minimum thresholds of the two-layer cloud height, as shown in formula (1);

[0009]

[0010] (2) Calculate the quadratic polynomial expansion coefficient of each cloud pixel. For the gray matrix of each center pixel and the set neighborhood window (PEW) (the window size is generally (2n+1) x (2n+1), and n is an integer) around it, the gray matrix can be regarded as a function f(x, y) of the spatial position of the independent variable, and the function is expanded by the quadratic polynomial through formula (2):

[0011]

[0012] Wherein, C=r1,

[0013] r1 is a constant term of polynomial expansion, and r2, r3, r4, r5 and r6 are coefficients of each term of polynomial expansion;

[0014] The expansion coefficients A, B, and C are estimated using weighted least squares. Since not every pixel's sample error in the neighborhood has the same influence on the center point, a two-dimensional Gaussian distribution is used in the function to assign weights to each.

[0015] (3) Calculate the displacement of each pixel. Based on the grayscale images of two consecutive time points, assuming the brightness of the center pixel at time T1 is f1(Xd) and the value after moving a distance d at time T2 is f2(X), if the brightness is constant, then we have formula (3):

[0016] f1(Xd)=C1+B1 T (Xd)+A1(Xd) T (Xd)=f2(X)=C2+B2 T X+A2X T X, (3)

[0017] Where A1, B1, C1, and A2, B2, C2 are the expansion coefficients of f1(Xd) and f2(X), respectively;

[0018] By rearranging the coefficients according to formulas (2) and (3), we can obtain:

[0019] A1 = A2, B2 = B1 - 2A1d, C2 = A1d T d+B1 T d+c1, (4)

[0020]

[0021] The displacement d can be solved if the above conditions are met. However, in practice, it is difficult to find a pixel that meets all the conditions, so an intermediate value is used. Substituting the estimated values ​​of A1 and A2 into the second constraint condition in formula (4), we get:

[0022]

[0023] d can be calculated using formula (6). However, the result obtained by this method contains too much noise. Therefore, a minimum variance objective function can be constructed within a larger local optimization window (LOW), as shown in formula (7):

[0024]

[0025] in, I indicates the element's position within the LOW window, w I d(X) represents the weight of each element in LOW with respect to the center element. d(X) is obtained by finding the minimum value of the objective function.

[0026] (4) The basic assumption of the dense optical flow method is that the image grayscale is continuously differentiable. Therefore, it may have a large error when the image motion is large. To solve this problem, this invention introduces the idea of ​​pyramid layering to achieve "coarse to fine" iteration for the optimal displacement d of the objective function. First, the displacement d is calculated starting from the coarsest scale. i To obtain a first-order but reasonable displacement estimate, and then according to formula (8), the scaling factor s and the displacement d estimated in the previous layer are adjusted. i The residual L is passed to the next layer at a finer scale to obtain a more accurate estimate.

[0027] d i+1 =s×d i +L, (8).

[0028] (5) Finally, based on the temporal and spatial resolution of the satellite images, the pixel displacements in the image coordinates are converted into velocity components (u, v) in geographic coordinates. For example, if the estimated displacement d has two directional components Δx and Δy, which actually represent the change in the pixel's position in the matrix, the converted displacements in geographic coordinates should be Δx×Δlon and Δy×Δlon, where Δlon and Δlat represent the actual geographic spatial resolution corresponding to a pixel in the satellite image. After converting to actual geographic displacement, dividing by the time difference between the two images (the temporal resolution of the satellite images) yields the velocity components (u, v) in geographic coordinates.

[0029] Determining the cloud top height using cloud-guided winds is the second most critical issue in traditional cloud-guided wind algorithms and a major source of error. This invention employs a DNN multilayer cloud inversion model within the optical flow method to determine the cloud top height, which is more accurate than the height obtained using the water vapor-infrared brightness temperature intercept method in other studies, especially in two-layer cloud scenarios.

[0030] The beneficial effects of this invention are as follows:

[0031] The inversion algorithm of this invention was applied to Typhoon Mulan, the 7th typhoon of the Pacific typhoon season in August 2022. The results show that for the hollow typhoon Mulan, the cloud-leading winds retrieved from the upper cloud top height are largely consistent with the results retrieved from the infrared brightness temperature channel; however, there is a significant offset between the cloud vortex center and the optimal path location. In contrast, the wind field spiral center retrieved from the lower cloud top height is largely consistent with the optimal path location at the corresponding time. The wind field inversion using double-layer cloud-leading winds will effectively improve the path and intensity forecasts of irregular typhoons. Attached Figure Description

[0032] Figure 1 This is a flowchart illustrating the dual-layer cloud wind inversion algorithm of the present invention. Detailed Implementation

[0033] The method of the present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0034] Figure 1 This invention demonstrates the cloud-guided wind inversion algorithm flow. The invention utilizes L1-level infrared multi-channel brightness temperature data from the Japanese Himawari-8 AHI sensor, with a spatial resolution of 5 km and a temporal resolution of 10 minutes. The multi-layer cloud detection and remote sensing inversion algorithm can obtain cloud phase and cloud top height data with the same spatiotemporal resolution, selecting moments where multi-layer clouds exist (in this invention, Typhoon Mulan is selected as an example). Since the optical flow method is based on pixel grayscale value changes, the inverted cloud height data undergoes quality control and standardization processing to ensure that the same cloud height value has the same grayscale value at different times.

[0035] In the upper-layer cloud quality control stage, for cloud heights less than 0.5km, all pixels are set to the default value (which is the minimum grayscale value in the optical flow algorithm), and for cloud heights greater than 18km, all pixel values ​​are calibrated to 18km. The cloud height data at all times are processed into a grayscale matrix using the two thresholds of 0.5 and 18.

[0036] In the lower-level cloud quality control stage, cloud heights below 0.5km are also treated as invalid values. Two thresholds, 0.5 and 8, are used to process all cloud height data into a grayscale matrix.

[0037] The processed grayscale matrix from two consecutive time points is input into the dense optical flow algorithm, with a neighborhood window size of 11×11. To better remove inversion noise and ensure the continuity of the wind field, a 50×50 window is used to smooth the displacement field. To uphold the fundamental assumption of continuous differentiability of image grayscale, a three-layer pyramid approach is introduced for iterative calculation from coarse to fine. After calculation using the dense optical flow method, the displacement change vector of each pixel based on the two consecutive time points can be obtained. To convert this into the U and V wind fields required for practical applications, the distance change vector needs to be converted into a velocity vector based on the satellite's spatial resolution of 5km and temporal resolution of 10min. Finally, combined with the cloud top height, the cloud-guided winds at different altitudes on the ice and water cloud layers can be obtained.

[0038] In a specific embodiment, the cloud-leading winds over Tropical Storm "Mulan" on August 9, 2022, were inverted. The inversion results of upper-level cloud top height and cloud-leading winds show that the cloud-leading wind inversion algorithm based on cloud top height can accurately invert the structure of the typhoon rotating cyclonically around its center, with the wind field at the eye area significantly smaller than the surrounding area. Furthermore, the cloud top height of the northern cloud system of "Mulan" is significantly smaller than that of the southern cloud system, and both the cloud area and wind field exhibit a "non-circular" shape. The spiral center of the upper-level "Mulan" is marked with a blue dot. Comparing the typhoon's position in the optimal path at the current moment, there is a significant difference between the spiral center of the upper-level cloud-leading winds and the position of the typhoon eye at the surface. The upper-level eye area is located northwest of the actual typhoon eye. The inversion results of lower-level cloud top height and cloud-leading winds show that low-level clouds exist at the typhoon's spiral center and in the outer spiral cloud system. The low-level cloud height at the spiral center is between 4 km and 5 km, while the low-level clouds in the outer spiral cloud system are distributed at heights between 2 km and 4 km. Since there are fewer usable pixels in low-level clouds, the mean window algorithm can provide an estimated wind field in cloudless areas by combining the movement of surrounding cloud systems. The inversion results show that the spiral center of the low-level cloud-guided winds is quite close to the optimal path typhoon eye location, indicating that the double-layer cloud-guided winds have significant application value in improving typhoon path and intensity forecasts.

[0039] Compared with the radiosonde observations, the cloud-guided wind speed deviation was 4.52 m / s, while the ERA5 wind speed error was slightly smaller at 3.48 m / s. Both ERA5 and cloud-guided winds underestimated the wind speed. In terms of wind direction, the cloud-guided winds matched the radiosonde winds better. Comparing the wind fields on the 8th and 9th, the ERA5 winds showed a significant deviation from the radiosonde winds, with an average deviation of 29.9 degrees and a maximum exceeding 90 degrees. In contrast, the average deviation between the cloud-guided winds and the radiosonde winds was only 12.3 degrees.

[0040] In summary, for typhoon systems with a two-layer cloud structure, the two-layer cloud-guided wind algorithm in this invention can effectively invert its wind field and determine the spiral center at different heights, accurately inverting the wind field of the lower typhoon core region. This is of great significance for improving the accuracy of typhoon path and intensity forecasts and promoting the application of cloud-guided wind data in numerical forecast assimilation.

[0041] The technical solution of the present invention is not limited to the specific embodiments described above; it can be applied to cloud systems with a two-layer cloud structure.

Claims

1. A two-layer cloud wind inversion algorithm based on two-layer cloud top height information, characterized in that, The specific steps are as follows: (1) Image preprocessing: The cloud top height images at different times are processed into a unified gray value matrix according to the same threshold range. Specifically, the maximum cloud top height of the upper layer is fixed at 18 km, and the maximum cloud top height of the lower double-layer water cloud is fixed at 8 km. The lowest cloud top height is set to 0.5 km. Pixels with cloud height less than 0.5 km are set to the default value, and pixels with cloud height greater than 18 km in the upper layer and greater than 8 km in the lower layer are all uniformly set to 18 km. After this preprocessing step, the numerical range of cloud height in the upper layer is between 0.5 and 18, and the numerical range of cloud height in the lower layer is between 0.5 and 8. In the upper and lower cloud layers, the cloud height is mapped to the gray value range of 0-255 by the maximum and minimum thresholds of the cloud height in each layer, as shown in formula (1): (2) Calculate the quadratic polynomial expansion coefficients for each cloud pixel; For each center pixel and the grayscale matrix of its surrounding neighborhood window PEW, it can be considered as a function f(x,y) of the spatial location of the independent variable, with a neighborhood window size of (2n+1)×(2n+1), where n is an integer; the function is expanded into a quadratic polynomial using formula (2): Where C = r1, r1 is the constant term of the polynomial expansion, and r2, r3, r4, r5, and r6 are the coefficients of each term in the polynomial expansion. The expansion coefficients A, B, and C are estimated using the weighted least squares method. (3) Calculate the displacement of each pixel; Based on the grayscale images of two consecutive time points, assuming the brightness of the center pixel at time T1 is f1(Xd) and the brightness after moving a distance d at time T2 is f2(X), if the assumption of constant brightness is satisfied, then we have formula (3): f1(X-d)=C1+B1 T (X-d)+A1(X-d) T (X-d)=f2(X)=C2+B2 T X+A2X T X,(3) Where A1, B1, C1, and A2, B2, C2 are the expansion coefficients of f1(Xd) and f2(X), respectively; By rearranging the coefficients according to formulas (2) and (3), we get: A1=A2,B2=B1-2A1d,C2=A1d T d-B1 T d+c1,(4) To solve for the displacement d while meeting the above conditions, an intermediate value is used. Substituting the estimated values ​​of A1 and A2 into the second constraint condition in formula (4), we get: To avoid noise in the results, a minimum variance objective function is constructed within a local optimization window LOW, specifically as shown in formula (7): in, I indicates the element's position within the LOW window, w I d(X) represents the weight of each element in LOW with respect to the center element. d(X) is obtained by finding the minimum value of the objective function. (4) For the optimal displacement d of the objective function, the idea of ​​pyramid layering is introduced to achieve iterative calculation from coarse to fine; firstly, the displacement d is calculated starting from the coarsest scale. i To obtain the first-order displacement estimate, and then according to formula (8), the scaling factor s and the displacement d estimated in the previous layer are adjusted. i The residual L is passed to the next layer at a finer scale to obtain a more accurate estimate; d i+1 <s×d i +L, (8) (5) Finally, based on the temporal and spatial resolution of the satellite image, the pixel displacement in the image coordinates is converted into velocity components (u, v) in the geographic coordinates; Let the two directional components of the estimated displacement d be Δx and Δy, which actually represent the change in the position of the pixel in the matrix. They are converted into displacements in geographic coordinates as Δx×Δlon and Δy×Δlon, where Δlon and Δlat represent the actual geospatial resolution corresponding to a pixel in the satellite image. After converting to actual geographic displacement, we divide by the time difference between the two images, which is the time resolution of the satellite image, to obtain the velocity components (u, v) in geographic coordinates.

Images