A typhoon wind field inversion method and system

Abstract

The present invention discloses a typhoon wind field inversion method and system. The method obtains the corresponding sector grids of each sector grid in the typhoon search area that have correlation at adjacent times by analyzing the change in echo intensity at adjacent times, thereby inverting the advection wind vector field and iteratively smoothing and correcting it, eliminating some points with incomplete velocity deblurring in the process of radar-based data quality control, so that the wind field has better integrity and continuity, and improves the feasibility and stability of business operation.

Description

Technical Field

[0001] The invention belongs to the technical field of atmospheric remote sensing monitoring and analysis, and relates to a typhoon wind field inversion method and system. Background Art

[0002] Doppler radar uses the Doppler effect to detect the position of the target, and calculates the target's speed relative to the radar by measuring the frequency changes of the electromagnetic waves emitted and returned by the target. Radar data covers a wide area, has high temporal and spatial resolution, and can provide reflectivity factors and radial velocity information under all-weather conditions. The new generation of Doppler radar networks deployed in my country's coastal areas can better detect the three-dimensional fine structure of typhoons approaching landfall, which greatly compensates for the lack of conventional observation data on the open ocean, and is of great significance for the monitoring, analysis and forecasting of typhoons. However, the radial velocity only includes the wind speed component in the radial direction of the radar, not the true wind speed vector. Therefore, scholars at home and abroad have studied and developed various radar wind field inversion methods.

[0003] At present, there are two main types of radar wind field inversion technologies suitable for typhoon circulation. One is the GBVTD method, which is based on radial velocity and uses geometric relationships and Fourier series transformation to invert the axisymmetric circulation and asymmetric tangential circulation structure in the inner core of the typhoon. However, it needs to meet the geometric restrictions of the algorithm and has a limited inversion range.

[0004] The other type is the TREC method, which estimates the three-dimensional horizontal flow field through the movement of radar echoes based on the reflectivity factor. It can obtain circulation information with a larger range than the radial velocity inversion wind field, and can be applied before the typhoon enters the observation range of the radar radial velocity. However, since the echoes in the typhoon eyewall area are usually strong and uniformly distributed, and in the strong and deep convection in the eyewall area, precipitation particles do not move with the horizontal wind, but are controlled by strong low-level convergence or high-level outflow. The wind field inverted by the TREC method is significantly underestimated in the typhoon eyewall area, resulting in weak circulation.

[0005] On this basis, the T-TREC method introduces radial velocity as an additional inversion condition, objectively selects the tangential search range and establishes a wind field correlation matrix, which effectively solves the problem of underestimated intensity caused by the relatively uniform echo in the typhoon eyewall area and reduces the inversion error. However, the inversion accuracy of the T-TREC method strongly depends on the quality of radar data, is extremely sensitive to velocity ambiguity, and has low feasibility in business operation. In addition, the method directly uses the observed radial velocity to correct the weak inverted radial velocity, which easily causes discontinuity in the wind field, resulting in discontinuity of various meteorological elements and unbalanced wind pressure relationships during analysis, affecting the analysis and application results. At the same time, the T-TREC method still underestimates the wind field, especially when the typhoon approaches land, the inversion accuracy is significantly reduced, and the velocity underestimated phenomenon is serious. The reason for this phenomenon is, on the one hand, the influence of ground object echoes, and on the other hand, the interaction between the circulation and terrain of the typhoon during landing continuously stimulates new convection. Summary of the invention

[0006] In order to solve the problems existing in the prior art, the present invention proposes a typhoon wind field inversion method and system.

[0007] The technical solution of the present invention is as follows:

[0008] A typhoon wind field inversion method, comprising:

[0009] Preprocessing the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids with a number of data points distributed therein;

[0010] Based on the echo intensity changes at adjacent moments, the corresponding sector grids of each sector grid in the typhoon search area with correlation at adjacent moments are obtained, and then the advection wind vector field is inverted;

[0011] The inverted advection wind vector field is iteratively smoothed and corrected, data points with incomplete velocity deblurring are eliminated, and the typhoon wind field is output.

[0012] Furthermore, the specific steps of preprocessing the radar-based data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and having a plurality of data points distributed therein include:

[0013] S11. Iterative positioning of the typhoon center:

[0014] S111. Determine the initial estimate of the typhoon center coordinates (X 0 ,Y 0 ), the initial radius R0 of the typhoon eye area, the radius variation range ΔR of the typhoon eye area and the weak echo threshold Z0, and set the number of iterations m=1, set the maximum number of iterations threshold β, the first convergence standard threshold α, and the second convergence standard threshold γ;

[0015] S112, obtain the typhoon center in this iteration as (X m ,Y m ),

[0016] in,

[0017] In the formula, (X, Y) represents any data point in the radar base data that meets the following two criteria; among them, criterion one is that the reflectivity of the data point (X, Y) is less than Z0; criterion two is (XX m-1 ) 2 +(YY m-1 ) 2 ≤R 2 , R is the actual radius of the typhoon eye area, R = R0-ΔR; A is the total area enclosed by the data points in the set of data points (X, Y);

[0018] S113. If any of the following convergence conditions is met, the current typhoon center (X m ,Y m ), otherwise, the number of iterations is increased by 1, the actual radius R is increased by 1, and the process returns to step S112;

[0019] Convergence condition 1: R = R0 + ΔR;

[0020] Convergence condition 2: The current number of iterations is greater than the maximum number of iterations threshold β;

[0021] Convergence condition three: (X m ,Y m ) and (X m-1 ,Y m-1 ) is less than the first convergence standard threshold α and the ratio of excessive points to valid points in the typhoon eye area boundary data points is greater than the second convergence standard threshold γ; the excessive points are the data points in the typhoon eye area boundary data points whose reflectivity values ​​are greater than the weak echo threshold Z0; the valid points are the data points in the typhoon eye area boundary data points whose reflectivity values ​​are less than or equal to the weak echo threshold Z0;

[0022] S12. Coordinate transformation:

[0023] The radar base data with the radar center as the origin is now converted into rectangular coordinate data with the typhoon center as the origin, and then further converted into polar coordinate data with a number of data points distributed with the typhoon center as the origin using the nearest interpolation method;

[0024] S13. Sector grid division:

[0025] According to a predetermined target, the arc length and radial length of the sector grid are determined, thereby dividing the polar coordinate data into sector grids.

[0026] Furthermore, the typhoon search area is determined by the following method:

[0027] Obtain the average tangential velocity of each data point in the radial distance library in the polar coordinate data. For any data point with a radial distance R in the radial distance library, if R is less than or equal to the distance between the radar and the typhoon center, the average tangential velocity of the data point is Among them, V rmax (R) is the maximum radial velocity of the data point at all azimuths; V rmin (R) is the minimum radial velocity of the data point at all azimuths; if R is greater than the distance between the radar and the typhoon center, the average tangential velocity v of the data point is T (R) = v T (R RMW )(R RMW / R) κ ; R RMW is the maximum wind speed radius of the tangential wind, K is an empirical constant;

[0028] The radial profile of the average tangential velocity is completed by using the observed maximum wind speed and the maximum wind speed radius, thereby obtaining the average tangential velocity of each data point on the radial profile;

[0029] Further determine the tangential reference search distance DA within the typhoon search area ref and radial reference search distance DR ref , thereby determining the search area of ​​the inverted wind vector end point; wherein the tangential reference search distance DA in the typhoon search area ref =V T (R)·Δt; radial reference search distance is determined as DR ref , DR ref =α·V T (R)·Δt, Δt is the time interval between adjacent moments of radar-based data, and α is the weight factor.

[0030] Furthermore, the specific method of obtaining the corresponding sector grids with correlation at adjacent times of each sector grid in the typhoon search area based on the echo intensity change at adjacent times includes:

[0031] S21, calculate the average echo intensity of each grid unit in the search area at adjacent times, and the calculation formula is:

[0032]

[0033] In the formula, E ave(r,θ,t) is the average echo intensity of the sector grid at position (r,θ) at time t, θ is the azimuth, and r is the distance between the sector grid and the typhoon center; Z(r,θ,i,t) is the echo intensity of the i-th data point in the grid unit at position (r,θ) at time t; n is the number of non-missing data points in the grid unit;

[0034] S22, when the difference in the average echo intensity of the grid cells in the search area at adjacent moments is less than the difference threshold E th When , the echo threshold is used to determine the corresponding sector grids with correlation at adjacent times in each sector grid in the typhoon search area, which is specifically expressed as:

[0035] E diff (r2,θ2,t+1)=|E ave (r2,θ2,t+1)-E ave (r1,θ1,t)| <E th

[0036] In the formula, E diff (r2,θ2,t+1) is the difference in the average echo intensity between the sector grid at any position (r2,θ2) and the sector grid at the starting position (r1,θ1) in the search area at adjacent moments. E ave (r2,θ2,t+1) is the average echo intensity of the sector grid at the position (r2,θ2) at the time t+1; E ave (r1,θ1,t) is the average echo intensity of the sector grid at the position (r1,θ1) at the tth moment.

[0037] Furthermore, the specific steps of inverting the advection wind vector field include:

[0038] S23, calculate the cross-correlation coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time according to the following formula, thereby forming a cross-correlation matrix; wherein, for the cross-correlation coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time z The calculation scheme is as follows:

[0039]

[0040] Where Z1(k) and Z2(k) are the reflectivity factor data of the kth data point in any sector grid and its corresponding sector grid with correlation at adjacent moments, respectively, and N is the number of data points contained in the sector grid;

[0041] S24, calculate the radial wind weight coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time, thereby forming the radial wind weight coefficient; for the radial wind weight coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time V The calculation scheme is as follows:

[0042]

[0043] Among them, DA ref is the tangential reference search distance; β is the weight factor; DA is the distance between the sector grid and its corresponding sector grid with correlation at adjacent moments.

[0044] S25. According to the cross-correlation matrix and the radial wind weight coefficient matrix, a total correlation coefficient matrix is ​​obtained, and then the end point of the inversion wind vector is determined according to the position of the maximum value of the total correlation coefficient matrix, and the advection wind speed vector value is calculated through the time interval between adjacent moments of the radar base data and the change of the sector grid position, so as to complete the inversion of the advection wind vector field.

[0045] Furthermore, the specific steps of iteratively smoothing and correcting the inverted advection wind vector field, removing data points with incomplete velocity deblurring, and outputting the typhoon wind field include:

[0046] S31. Perform quality control on the inverted advection wind vector field, and fill in the wind field difference and missing values ​​in the rectangular coordinate system;

[0047] S32, the inverted advection wind vector field (u 0 ,v 0 ) is projected onto the radar radial direction to obtain the inverted radial wind component

[0048]

[0049] In the formula, dx and dy represent (u 0 ,v 0 ) is located at a distance from the center of the radar in the x-axis and y-axis directions, and r is (u 0 ,v 0 ) is the distance between the grid point and the radar center.

[0050] S33, the inverted radial wind component Radar observation of radial wind v r For comparison, if Then iterative correction is performed; otherwise, the output (u 0 ,v 0 ); The specific steps of iterative smoothing correction include:

[0051] S331, if The wind direction Dir of the inverted advection wind vector is correct, and the wind speed Spd 0 Weak, set the number of iterations m = 1;

[0052] S332, update the wind speed with a fixed step size Δs, and keep the wind direction Dir unchanged. The updated wind speed is recorded as Spd m , namely Spd m =Spd 0 +m×Δs,

[0053] Calculate the updated wind speed component u m 、v m , the calculation formula is as follows:

[0054]

[0055] S333、(u m ,v m ) projected onto the radar radial direction to obtain Recalculate the feature quantity m is the number of iterations,

[0056] S334, when the feature quantity F<ε or m>β, the iteration is stopped, ε is the convergence threshold, β is the iteration number threshold; if the iteration is stopped when m>β, the u that minimizes F is selected m 、v m Replace u 0 、v 0 , and the final output of the typhoon three-dimensional wind field is obtained; if the iteration is stopped when F<ε, the u at the time of iteration stop is output m 、v m Replace u 0 、v 0 , and obtain the final output typhoon three-dimensional wind field; otherwise, set m=m+1 and return to step S332.

[0057] Furthermore, the specific steps of performing quality control on the inverted advection wind vector field and filling in the wind field difference and missing values ​​in the rectangular coordinate system include:

[0058] Inversion wind quality control: Determine the rectangular coordinate grid with the typhoon center as the origin for the final output of the inverted advection wind vector, and average the inverted advection wind speed within the search neighborhood radius of each grid point. If the absolute value of the difference between a certain inverted advection wind speed and the average value within the range exceeds the threshold, the inverted advection wind vector is discarded;

[0059] Rectangular coordinate wind field interpolation: For each grid point, if the number of valid advection wind vectors within its search neighborhood radius is not less than the preset minimum value and is distributed in at least two quadrants, the advection wind vector of the grid point is interpolated using the least squares method.

[0060] Filling in missing values: Traverse the rectangular coordinate grid points and screen out the missing grid points. For each missing grid point, determine its search neighborhood radius. When the number of valid advection wind vectors within the search range is not less than the preset minimum value, use least squares interpolation to fill in the missing values.

[0061] A typhoon wind field inversion system comprises a preprocessing module, an inversion module and a smoothing correction module;

[0062] The preprocessing module is used to preprocess the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and distributed with a number of data points;

[0063] The inversion module is used to obtain the corresponding sector grids with correlation between the sector grids in the typhoon search area at adjacent times based on the echo intensity changes at adjacent times, and then perform the inversion of the advection wind vector field;

[0064] The smoothing correction module is used to perform iterative smoothing correction on the inverted advection wind vector field, eliminate data points with incomplete velocity deblurring, and output the typhoon wind field.

[0065] An electronic device comprises a memory and a processor, wherein the memory stores a computer program, and the processor is used to call and run the computer program stored in the memory to execute any of the methods described above.

[0066] A computer-readable storage medium stores a computer program, wherein the computer program implements the steps of any of the above methods when executed by a processor.

[0067] Compared with the prior art, the present invention has the following beneficial effects:

[0068] The present invention proposes a typhoon wind field inversion method and system. The method obtains the corresponding sector grids of each sector grid in the typhoon search area that have correlation at adjacent times by analyzing the change in echo intensity at adjacent times, thereby inverting the advection wind vector field and iteratively smoothing and correcting it, eliminating some points with incomplete velocity deblurring in the process of radar-based data quality control, so that the wind field has better integrity and continuity, and improves the feasibility and stability of business operation.

[0069] The typhoon wind field inversion method of the present invention projects the advection wind vector field in the radial direction of the radar to obtain the radial wind component, and then iteratively smoothes and corrects the radial wind component, taking into account the influence of the radial wind on the wind field, thereby reducing the error of wind field inversion and making the inverted typhoon three-dimensional wind field smoother and more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 Schematic diagram of the process of typhoon wind field inversion method in the embodiment. DETAILED DESCRIPTION

[0071] Embodiment 1:

[0072] A typhoon wind field inversion method of the present invention is as follows: Figure 1 As shown, including:

[0073] S1. Preprocess the radar-based data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids with a number of data points distributed therein; wherein, the purpose of preprocessing is to perform quality control on the radar-based data, and the preprocessing includes processing such as de-velocity blurring, eliminating ground clutter, and removing isolated noise points. In this example, the radar-based data is radar PPI data. The preprocessed radar PPI data is interpolated onto the contour surface to obtain CAPPI data with the radar center as the origin. The typhoon center can be objectively located by identifying the geometric center of the radar echo eye area. The CAPPI data with the radar center as the origin can be converted into polar coordinate data with the typhoon center as the origin. Grid division is then performed to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids with a number of data points distributed therein.

[0074] S2. Based on the change in echo intensity at adjacent moments, the corresponding sector grids of each sector grid in the typhoon search area with correlation at adjacent moments are obtained, and then the advection wind vector field is inverted. The present invention performs echo threshold judgment on the grid units in the search area at adjacent moments, and then calculates the cross-correlation coefficient matrix and radial wind weight coefficient matrix between the corresponding sector grids with correlation at adjacent moments. Then, the total correlation coefficient matrix is ​​obtained by multiplying the cross-correlation coefficient matrix and the radial wind weight coefficient matrix, and the end point of the inverted wind vector is determined according to the position of the maximum value of the total correlation coefficient matrix. The advection wind speed is calculated by the time difference between adjacent moments and the change in the position of the sector grid, and the advection wind vector field is inverted.

[0075] S3. Iteratively smooth and correct the inverted advection wind vector field, remove data points with incomplete velocity deblurring, and output the typhoon wind field.

[0076] Embodiment 2:

[0077] This embodiment is further designed on the basis of the first embodiment in that, in this embodiment, the radar base data is preprocessed to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and having a plurality of data points distributed therein, including:

[0078] S11. Iterative positioning of the typhoon center:

[0079] S111, assuming that the typhoon eye is circular, determine the initial estimate of the typhoon center coordinates (X 0 ,Y 0 ), the initial radius R0 of the typhoon eye area, the radius variation range ΔR of the typhoon eye area and the weak echo threshold Z0, and set the number of iterations m=1, set the maximum number of iterations threshold β, the first convergence standard threshold α, and the second convergence standard threshold γ;

[0080] S112, obtain the typhoon center in this iteration as (X m ,Y m ),

[0081] in,

[0082] In the formula, (X, Y) represents any data point in the radar base data that meets the following two criteria; among them, criterion one is that the reflectivity of the data point (X, Y) is less than Z0; criterion two is (XX m-1 ) 2 +(YY m-1 ) 2 ≤R 2 , R is the actual radius of the typhoon eye area, R = R0-ΔR; A is the total area enclosed by the data points in the set of data points (X, Y);

[0083] S113. If any of the following convergence conditions is met, the current typhoon center (X m ,Y m ), otherwise, the number of iterations is increased by 1 (i.e., m=m+1), the actual radius R is increased by 1 (i.e., R=R+1), and the process returns to step S112;

[0084] Convergence condition 1: R = R0 + ΔR;

[0085] Convergence condition 2: The current number of iterations is greater than the maximum number of iterations threshold β;

[0086] Convergence condition three: (X m ,Y m ) and (X m-1 ,Y m-1) is less than the first convergence standard threshold α and the ratio of excessive points to valid points in the typhoon eye area boundary data points is greater than the second convergence standard threshold γ; the excessive points are the data points in the typhoon eye area boundary (outermost layer) whose reflectivity values ​​are greater than the weak echo threshold Z0; the valid points are the data points in the typhoon eye area boundary (outermost layer) whose reflectivity values ​​are less than or equal to the weak echo threshold Z0;

[0087] S12. Coordinate transformation:

[0088] The radar base data with the radar center as the origin is now converted into rectangular coordinate data with the typhoon center as the origin, and then further converted into polar coordinate data with a number of data points distributed with the typhoon center as the origin using the nearest interpolation method;

[0089] S13. Sector grid division:

[0090] According to the predetermined target, the arc length and radial length of the sector grid are determined, thereby dividing the polar coordinate data into sector grids.

[0091] Embodiment three:

[0092] This embodiment is further designed on the basis of the second embodiment in that the typhoon search area in this embodiment is determined by the following method:

[0093] Get the average tangential velocity of each data point in the radial distance library in the polar coordinate data. For any data point with a radial distance R in the radial distance library, if R is less than or equal to the distance between the radar and the typhoon center, the average tangential velocity of the data point is Among them, V rmax (R) is the maximum radial velocity of the data point at all azimuths; V rmin (R) is the minimum radial velocity of the data point at all azimuths; if R is greater than the distance between the radar and the typhoon center, the average tangential velocity of the data point (, the Rankine vortex model is used for fitting, and the fitting formula is) V T (R) = V T (R RMW )(R RMW / R) κ ; R RMW is the maximum wind speed radius of the tangential wind, K is an empirical constant;

[0094] The observed maximum wind speed and the maximum wind speed radius are used to complete the radial profile of the average tangential velocity, thereby obtaining the average tangential velocity of each data point on the radial profile; the observed maximum wind speed and the maximum wind speed radius can be obtained from the typhoon dataset.

[0095] Further determine the tangential reference search distance DA within the typhoon search area refand radial reference search distance DR ref , thereby determining the search area of ​​the inverted wind vector end point; among them, the tangential reference search distance DA in the typhoon search area ref =V T (R)·Δt; radial reference search distance is determined as DR ref , DR ref =α·V T (R)·Δt, Δt is the time interval between adjacent moments of radar-based data, and α is the weight factor.

[0096] Embodiment 4:

[0097] This embodiment is further designed on the basis of the third embodiment in that, in this embodiment, based on the echo intensity change at adjacent moments, a specific method for obtaining corresponding sector grids having correlation at adjacent moments of each sector grid in the typhoon search area includes:

[0098] S21, calculate the average echo intensity of each grid unit in the search area at adjacent times, and the calculation formula is:

[0099]

[0100] In the formula, E ave (r,θ,t) is the average echo intensity of the sector grid at position (r,θ) at time t, θ is the azimuth, and r is the distance between the sector grid and the typhoon center; Z(r,θ,i,t) is the echo intensity of the i-th data point in the grid unit at position (r,θ) at time t; n is the number of non-missing data points in the grid unit;

[0101] S22, when the difference in the average echo intensity of the grid cells in the search area at adjacent moments is less than the difference threshold E th When , the echo threshold is used to determine the corresponding sector grids with correlation at adjacent times in each sector grid in the typhoon search area, which is specifically expressed as:

[0102] E diff (r2,θ2,t+1)=|E ave (r2,θ2,t+1)-E ave (r1,θ1,t)| <E th

[0103] In the formula, E diff (r2,θ2,t+1) is the difference in the average echo intensity between the sector grid at any position (r2,θ2) and the sector grid at the starting position (r1,θ1) in the search area at adjacent moments. E ave(r2,θ2,t+1) is the average echo intensity of the sector grid at the position (r2,θ2) at the time t+1; E ave (r1,θ1,t) is the average echo intensity of the sector grid at the position (r1,θ1) at the tth moment.

[0104] Embodiment five:

[0105] This embodiment is further designed on the basis of the first embodiment in that the specific steps of inverting the advection wind vector field in this embodiment include:

[0106] S23, calculate the cross-correlation coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time according to the following formula, thereby forming a cross-correlation matrix; wherein, for the cross-correlation coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time z The calculation scheme is as follows:

[0107]

[0108] Where Z1(k) and Z2(k) are the reflectivity factor data of the kth data point in any sector grid and its corresponding sector grid with correlation at adjacent moments, respectively, and N is the number of data points contained in the sector grid;

[0109] S24, calculate the radial wind weight coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time, thereby forming the radial wind weight coefficient; the actual tangential speed is V T (R) fluctuates in the neighborhood, so the search area is set to be within the tangential reference search distance DA from the previous grid cell. ref The next grid unit in the adjacent domain has a weight of 1, and the weight decreases with the distance to the sides, that is, the radial wind weight coefficient ρ between any sector grid and its corresponding sector grid with correlation at the adjacent time V The calculation scheme is as follows:

[0110]

[0111] Among them, DA ref is the tangential reference search distance; β is the weight factor; DA is the distance between the sector grid and its corresponding sector grid with correlation at adjacent moments.

[0112] S25. According to the cross-correlation matrix and the radial wind weight coefficient matrix, the total correlation coefficient matrix is ​​obtained, and then the end point of the inversion wind vector is determined according to the position of the maximum value of the total correlation coefficient matrix. The advection wind speed vector value is calculated through the time interval between adjacent moments of the radar base data and the change of the sector grid position, and the inversion of the advection wind vector field is completed.

[0113] Embodiment six:

[0114] This embodiment is further designed on the basis of the first embodiment in that, in this embodiment, the inverse advection wind vector field is iteratively smoothed and corrected, and data points with incomplete velocity deblurring are eliminated. The specific steps of outputting the typhoon wind field include:

[0115] S31. Perform quality control on the inverted advection wind vector field, and fill in the wind field difference and missing values ​​in the rectangular coordinate system;

[0116] S32, the inverted advection wind vector field (u 0 ,v 0 ) is projected onto the radar radial direction to obtain the inverted radial wind component

[0117]

[0118] In the formula, dx and dy represent (u 0 ,v 0 ) is located at a distance from the center of the radar in the x-axis and y-axis directions, and r is (u 0 ,v 0 ) is the distance between the grid point and the radar center.

[0119] S33, the inverted radial wind component Radar observation of radial wind V r For comparison, if Then iterative correction is performed; otherwise, the output (u 0 ,v 0 ); The specific steps of iterative smoothing correction include:

[0120] S331, if The wind direction Dir of the inverted advection wind vector is correct, and the wind speed Spd 0 Weak, set the number of iterations m = 1;

[0121] S332. In each iteration, the wind speed is updated with a fixed step size Δs, and the wind direction Dir is kept unchanged. The updated wind speed is recorded as Spd m , namely Spd m =Spd 0 +m×Δs,

[0122] Calculate the updated wind speed component u m 、v m , the calculation formula is as follows:

[0123]

[0124] S333、(um ,v m ) projected onto the radar radial direction to obtain Recalculate the feature quantity m is the number of iterations,

[0125] S334, when the feature quantity F<ε or m>β, the iteration is stopped, ε is the convergence threshold, β is the iteration number threshold; if the iteration is stopped when m>β, the u that minimizes F is selected m 、v m Replace u 0 、v 0 , and the final output of the typhoon three-dimensional wind field is obtained; if the iteration is stopped when F<ε, the u at the time of iteration stop is output m 、v m Replace u 0 、v 0 , and obtain the final output typhoon three-dimensional wind field; otherwise, set m=m+1 and return to step S332.

[0126] Embodiment seven:

[0127] This embodiment is further designed on the basis of the sixth embodiment in that the specific steps of quality control of the inverted advection wind vector field and filling in the wind field difference and missing values ​​in the rectangular coordinate system in this embodiment include:

[0128] Inversion wind quality control: Determine the rectangular coordinate grid with the typhoon center as the origin of the final output of the inverted advection wind vector, and average the inverted advection wind speed within the search neighborhood radius of each grid point. If the absolute value of the difference between a certain inverted advection wind speed and the average value within this range exceeds the threshold, the inverted advection wind vector is eliminated; this design is to ensure the credibility and accuracy of the data.

[0129] Rectangular coordinate wind field interpolation: For each grid point, if the number of valid advection wind vectors within its search neighborhood radius is not less than the preset minimum value and is distributed in at least two quadrants, the advection wind vector of the grid point is interpolated using the least squares method.

[0130] Filling in missing values: Traverse the rectangular coordinate grid points and screen out the missing grid points. For each missing grid point, determine its search neighborhood radius. When the number of valid advection wind vectors within the search range is not less than the preset minimum value, use least squares interpolation to fill in the missing values.

[0131] Embodiment eight:

[0132] A typhoon wind field inversion system of the present invention comprises a preprocessing module, an inversion module and a smoothing correction module;

[0133] A preprocessing module is used to preprocess the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and distributed with a number of data points;

[0134] An inversion module is used to obtain the corresponding sector grids with correlation at adjacent times of each sector grid in the typhoon search area based on the echo intensity changes at adjacent times, and then invert the advection wind vector field;

[0135] The smoothing correction module is used to iteratively smooth the inverted advection wind vector field, eliminate data points with incomplete velocity deblurring, and output the typhoon wind field.

[0136] Embodiment nine:

[0137] An electronic device of the present invention includes a memory and a processor, wherein the memory stores a computer program, and the processor is used to call and run the computer program stored in the memory to execute the method described in any of the above embodiments.

[0138] A computer-readable storage medium of the present invention stores a computer program, and when the computer program is executed by a processor, the steps of the method described in any of the above embodiments are implemented.

[0139] Application examples:

[0140] This example simulates and implements the typhoon wind field inversion method of the present invention. The simulation object of this example is the typhoon "Lotus" at 18:00 (UTC) on July 8, 2015. The specific steps of the method are as follows:

[0141] S1: Quality control of radar-based data is performed, including velocity blur removal, ground clutter removal, and isolated noise point removal. The radar PPI data after quality control are interpolated onto a contour surface with a horizontal and vertical grid resolution of 1 km in a Cartesian rectangular coordinate system to obtain CAPPI data within the range of -300 km to 300 km centered on the radar at an altitude of 1 to 8 km.

[0142] S2: Objectively locate the typhoon center by identifying the geometric center of the radar echo eye area, and convert the CAPPI data with the radar center as the origin into polar coordinate data with the typhoon center as the origin. The specific contents include:

[0143] S21: Typhoon center positioning. Assuming that the typhoon eye is circular, the initial estimate of the typhoon center (X 0 ,Y 0 ), eye radius R0, eye radius variation range ΔR and weak echo threshold Z0 = 15dBZ. In each iteration, the typhoon center is (X m ,Y m),in m is the number of iterations, (X,Y) represents any point whose reflectivity is less than Z0 and satisfies (XX m-1 ) 2 +(YY m-1 ) 2 ≤R 2 , is the total area of ​​all data points that meet the above two criteria. In the first iteration, R = R0-ΔR, and R increases by 1 in subsequent iterations. When the number of iterations is greater than the threshold or both convergence conditions are met at the same time, the iteration is stopped when the iteration convergence condition is met to obtain the final typhoon center position.

[0144] Convergence condition 1: R = R0 + ΔR;

[0145] Convergence condition 2: The current number of iterations is greater than the maximum number of iterations threshold β = 30;

[0146] Convergence condition three: (X m ,Y m ) and (X m-1 ,Y m-1 ) is less than the first convergence standard threshold α=0.5 and the ratio of excessive points to valid points in the typhoon eye boundary data points is greater than the second convergence standard threshold γ=0.8; the excessive points are the data points in the typhoon eye boundary (outermost layer) whose reflectivity values ​​are greater than the weak echo threshold Z0; the valid points are the data points in the typhoon eye boundary (outermost layer) whose reflectivity values ​​are less than or equal to the weak echo threshold Z0.

[0147] S22: Coordinate conversion. First, the CAPPI data with the radar center as the origin is converted into rectangular coordinate system data with the typhoon center as the origin, and then converted into polar coordinate system data with the typhoon center as the origin using the nearest neighbor interpolation method.

[0148] S3: Determine the parameters of the fan-shaped grid unit, use the radial velocity to calculate the radial profile of the average tangential velocity to determine the search area, perform echo threshold determination on the grid units in the search area at adjacent moments, and then calculate the cross-correlation coefficient matrix and radial wind weight coefficient matrix. The specific contents include:

[0149] S31: Determine the grid unit parameters. The centerline arc length and radial length of the fan-shaped grid unit are related to the strength and size of the typhoon. Stronger and smaller typhoons should use smaller grid units, and vice versa, larger grid units should be used, and the centerline arc length is greater than the radial length. The grid unit size selected in this example is 10km×25km (radial×tangential).

[0150] S32: Determine the search area. Assuming that the typhoon circulation meets the Rankine vortex condition, for each radial distance bin, calculate its maximum radial velocity V at all azimuths. rmaxand minimum radial velocity V rmin , then the azimuth average tangential velocity with the typhoon center as the origin and the radial distance R can be expressed as V T is the average tangential velocity, which varies with the radius R. When R exceeds the distance between the radar and the typhoon center, the average tangential velocity is fitted using the Rankine vortex model, V T (R) = V T (R RMW )(R RMW / R) κ , R RMW is the maximum wind speed radius of the tangential wind; finally, the radial profile of the average tangential speed is completed using the observed maximum wind speed and radius. The tangential reference search distance is determined as DA ref , D.A. ref =V T (R)·Δt, the radial reference search distance is determined as DR ref , DR ref =Δ·V T (R) · Δt, Δt = 6min is the time interval between adjacent moments of radar data, and α = 0.3 is the weight factor. The tangential reference search distance DA in the typhoon search area ref and radial reference search distance DR ref , thereby determining the search area for the end point of the inverted wind vector;

[0151] S33: Echo threshold determination. Calculate the average echo intensity E of the grid cells in the search area at adjacent times. ave , n is the number of non-missing data points in the grid unit. When the difference in the average echo intensity of the grid units in the search area at adjacent moments is less than a certain threshold, the echo threshold is used for judgment, that is,

[0152] E diff (r2,θ2,t+1)=|E ave (r2,θ2,t+1)-E ave (r1,θ1,t)| <E th , E diff It represents the average echo intensity E of any position (r2, θ2) and the starting position (r1, θ1) in the search area starting from position (r1, θ1) at adjacent moments. ave difference,

[0153] E th =5dBZ is the difference threshold.

[0154] S34: Calculation of cross-correlation coefficient. According to the formula Calculate the cross-correlation coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time z , Z1 and Z2 are the reflectivity factor data of the grid cells of the previous and next two times, respectively, N is the number of data points contained in the grid cell, and the corresponding cross-correlation matrix ρ can be obtained by calculating the grid cells of the previous time and the grid cells of the next time in the search area respectively. z .

[0155] S35: Calculation of radial wind weight coefficient. The actual tangential velocity is the average tangential velocity V T (R) fluctuates in the neighborhood, so the search area is set to be within DA ref The next grid unit in the adjacent domain has a weight of 1, and the weight decreases with the distance to the sides. The radial wind weight coefficient ρ between any sector grid and its corresponding sector grid with correlation at the adjacent time is calculated as follows: V :

[0156]

[0157] β=0.3 is a weight factor; DA is the distance between a sector grid and its corresponding sector grid with correlation at adjacent moments.

[0158] S4: The total correlation coefficient matrix is ​​obtained by multiplying the cross-correlation coefficient matrix and the radial wind weight coefficient matrix. The end point of the inversion wind vector is determined according to the position of the maximum value of the total correlation coefficient matrix, and the advection wind speed is calculated by the time difference between adjacent moments and the change in the position of the grid unit.

[0159] S5: After quality control of the inverted advection wind vector field, it is interpolated to the grid points to obtain the rectangular coordinate wind field with the typhoon center as the origin, and the missing values ​​are filled. The inverted wind is projected to the radial direction of the radar to obtain the inverted radial wind component. If the inverted radial wind component is weaker than the observed radial wind, it is corrected by iterative smoothing to obtain the final typhoon wind field. The specific contents include:

[0160] S51: Inversion wind quality control. Determine the rectangular coordinate grid with the typhoon center as the origin for the final output of the inverted advection wind vector, with a horizontal resolution of 5km, a vertical resolution of 1km, and a vertical range of 1 to 8km. For each grid point, search for the average of the inverted advection wind speed within a neighborhood radius of 10km. If the absolute value of the difference between a certain inverted advection wind speed and the average value within this range exceeds the threshold of 10m / s, the inverted advection wind vector is discarded to ensure the credibility and accuracy of the data.

[0161] S52: Cartesian coordinate system wind field interpolation. For each grid point, if the number of valid advection wind vectors within the search neighborhood radius of 10 km is not less than a preset minimum value of 3 and is distributed in at least two quadrants, the advection wind vector of the grid point is interpolated using the least squares method.

[0162] S53: Filling in missing values. Traverse the rectangular coordinate grid points to screen out missing grid points. For each missing grid point, determine its search neighborhood radius as 10km. When the number of valid advection wind vectors within the search range is not less than the preset minimum value of 8, use least squares interpolation to fill in the missing values.

[0163] S54: Iterative smoothing correction. Using the formula The inverted advective wind vector field (u 0 ,v 0 ) is projected onto the radar radial direction to obtain the inverted radial wind component Where dx and dy represent (u 0 ,v 0 ) is located at a distance from the center of the radar in the x-axis and y-axis directions, and r is (u 0 ,v 0 ) is located at the distance from the radar center. The inverted radial wind component Radar observation of radial wind V r For comparison, if Iterative correction is then performed. Assuming that the wind direction Dir of the inverted advection wind vector is correct, the wind speed Spd 0 Weaker, in each iteration, the wind speed is increased with a fixed step size of Δs = 0.5 m / s (keeping the wind direction Dir unchanged) to obtain Spd m , namely Spd m =Spd 0 +m×Δs, through Calculate the updated wind speed component u m 、v m , and projected onto the radar radial direction to obtain calculate m is the number of iterations. When F<ε or m>β, the iteration stops. ε=0.1m / s is the convergence threshold, β=20 is the iteration threshold.

[0164] When the iteration stops when m>β, the u that minimizes F is selected. m 、v m Replace u 0 、v 0 , and obtain the final output typhoon three-dimensional wind field.

Claims

1. A typhoon wind field inversion method, characterized in that: include: Preprocessing the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids with a number of data points distributed therein; Based on the echo intensity changes at adjacent moments, the corresponding sector grids of each sector grid in the typhoon search area with correlation at adjacent moments are obtained, and then the advection wind vector field is inverted; Iteratively smooth and correct the inverted advection wind vector field, remove data points with incomplete velocity deblurring, and output the typhoon wind field; The specific steps of iteratively smoothing and correcting the inverted advection wind vector field, eliminating data points with incomplete velocity deblurring, and outputting the typhoon wind field include: S31. Perform quality control on the inverted advection wind vector field, and fill in the wind field difference and missing values ​​in the rectangular coordinate system; S32, the inverted advection wind vector field (u 0 ,v 0 ) is projected onto the radar radial direction to obtain the inverted radial wind component In the formula, dx and dy represent (u 0 ,v 0 ) is located at a distance from the center of the radar in the x-axis and y-axis directions, and r is (u 0 ,v 0 ) is located at a grid point and the radar center; S33, the inverted radial wind component Radar observation of radial wind V r For comparison, if Then iterative correction is performed; otherwise, the output (u 0 ,v 0 ); The specific steps of iterative smoothing correction include: S331, if The wind direction Dir of the inverted advection wind vector is correct, and the wind speed Spd 0 Weak, set the number of iterations m = 1; S332, update the wind speed with a fixed step size Δs, and keep the wind direction Dir unchanged. The updated wind speed is recorded as Spd m , namely Spd m =Spd 0 +m×Δs, Calculate the updated wind speed component u m 、v m , the calculation formula is as follows: S333、(u m ,v m ) is projected onto the radar radial direction to obtain Recalculate the feature quantity m is the number of iterations, S334, when the feature quantity F<ε or m>β, the iteration is stopped, ε is the convergence threshold, β is the iteration number threshold; if the iteration is stopped when m>β, the u that minimizes F is selected m 、v m Replace u 0 、v 0 , and the final output of the typhoon three-dimensional wind field is obtained; if the iteration is stopped when F<ε, the u at the time of iteration stop is output m 、v m Replace u 0 、v 0 , and obtain the final output typhoon three-dimensional wind field; otherwise, set m=m+1 and return to step S332.

2. The typhoon wind field inversion method according to claim 1, characterized in that: The specific steps of preprocessing the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and distributed with a number of data points include: S11. Iterative positioning of the typhoon center: S111. Determine the initial estimate of the typhoon center coordinates (X 0 ,Y 0 ), the initial radius R0 of the typhoon eye area, the radius variation range ΔR of the typhoon eye area and the weak echo threshold Z0, and set the number of iterations m=1, set the maximum number of iterations threshold β, the first convergence standard threshold α, and the second convergence standard threshold γ; S112, obtain the typhoon center in this iteration as (X m ,Y m ), in, In the formula, (X, Y) represents any data point in the radar base data that meets the following two criteria; among them, criterion one is that the reflectivity of the data point (X, Y) is less than Z0; criterion two is (XX m-1 ) 2 +(YY m-1 ) 2 ≤R 2 , R is the actual radius of the typhoon eye area, R = R0-ΔR; A is the total area enclosed by the data points in the set of data points (X, Y); S113. If any of the following convergence conditions is met, the current typhoon center (X m ,Y m ), otherwise, the number of iterations is increased by 1, the actual radius R is increased by 1, and the process returns to step S112; Convergence condition 1: R = R0 + ΔR; Convergence condition 2: The current number of iterations is greater than the maximum number of iterations threshold β; Convergence condition three: (X m ,Y m ) and (X m-1 ,Y m-1 ) is less than the first convergence standard threshold α and the ratio of excessive points to valid points in the typhoon eye area boundary data points is greater than the second convergence standard threshold γ; the excessive points are the data points in the typhoon eye area boundary data points whose reflectivity values ​​are greater than the weak echo threshold Z0; the valid points are the data points in the typhoon eye area boundary data points whose reflectivity values ​​are less than or equal to the weak echo threshold Z0; S12. Coordinate transformation: The radar base data with the radar center as the origin is now converted into rectangular coordinate data with the typhoon center as the origin, and then further converted into polar coordinate data with a number of data points distributed with the typhoon center as the origin using the nearest interpolation method; S13. Sector grid division: According to a predetermined target, the arc length and radial length of the sector grid are determined, thereby dividing the polar coordinate data into sector grids.

3. The typhoon wind field inversion method according to claim 2, characterized in that: The typhoon search area is determined by the following method: Obtain the average tangential velocity of each data point in the radial distance library in the polar coordinate data. For any data point with a radial distance R in the radial distance library, if R is less than or equal to the distance between the radar and the typhoon center, the average tangential velocity of the data point is Among them, V rmax (R) is the maximum radial velocity of the data point at all azimuths; V rmin (R) is the minimum radial velocity of the data point at all azimuths; if R is greater than the distance between the radar and the typhoon center, the average tangential velocity V of the data point is T (R) = V T (R RMW )(R RMW / R) κ ; R RMW is the maximum wind speed radius of the tangential wind, K is an empirical constant; The radial profile of the average tangential velocity is completed by using the observed maximum wind speed and the maximum wind speed radius, thereby obtaining the average tangential velocity of each data point on the radial profile; Further determine the tangential reference search distance DA within the typhoon search area ref and radial reference search distance DR ref , thereby determining the search area of ​​the inverted wind vector end point; wherein the tangential reference search distance DA in the typhoon search area ref =V T (R)·Δt; radial reference search distance is determined as DR ref , DR ref =α·V T (R)·Δt, Δt is the time interval between adjacent moments of radar-based data, and α is the weight factor.

4. The typhoon wind field inversion method according to claim 3, characterized in that: The specific method of obtaining the corresponding sector grids with correlation at adjacent times of each sector grid in the typhoon search area based on the echo intensity change at adjacent times includes: S21, calculate the average echo intensity of each grid unit in the search area at adjacent times, and the calculation formula is: In the formula, E ave (r,θ,t) is the average echo intensity of the sector grid at position (r,θ) at time t, θ is the azimuth, and r is the distance between the sector grid and the typhoon center; Z(r,θ,i,t) is the echo intensity of the i-th data point in the grid unit at position (r,θ) at time t; n is the number of non-missing data points in the grid unit; S22, when the difference in the average echo intensity of the grid cells in the search area at adjacent moments is less than the difference threshold E th When , the echo threshold is used to determine the corresponding sector grids with correlation at adjacent times in each sector grid in the typhoon search area, which is specifically expressed as: E diff (r2,θ2,t+1)=|E ave (r2,θ2,t+1)-E ave (r1,θ1,t)|<E th In the formula, E diff (r2,θ2,t+1) is the difference in the average echo intensity between the sector grid at any position (r2,θ2) and the sector grid at the starting position (r1,θ1) in the search area at adjacent moments. E ave (r2,θ2,t+1) is the average echo intensity of the sector grid at the position (r2,θ2) at the time t+1; E ave (r1,θ1,t) is the average echo intensity of the sector grid at the position (r1,θ1) at the tth moment.

5. The typhoon wind field inversion method according to claim 1, characterized in that: The specific steps of inverting the advective wind vector field include: S23, calculate the cross-correlation coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time according to the following formula, thereby forming a cross-correlation matrix; wherein, for the cross-correlation coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time z The calculation scheme is as follows: Where Z1(k) and Z2(k) are the reflectivity factor data of the kth data point in any sector grid and its corresponding sector grid with correlation at adjacent moments, respectively, and N is the number of data points contained in the sector grid; S24, calculate the radial wind weight coefficient between each sector grid and its corresponding sector grid with correlation at adjacent time, thereby forming the radial wind weight coefficient; for the radial wind weight coefficient ρ between any sector grid and its corresponding sector grid with correlation at adjacent time V The calculation scheme is as follows: Among them, DA ref is the tangential reference search distance; β is the weight factor; DA is the distance between the sector grid and its corresponding sector grid with correlation at adjacent time; S25. According to the cross-correlation matrix and the radial wind weight coefficient matrix, a total correlation coefficient matrix is ​​obtained, and then the end point of the inversion wind vector is determined according to the position of the maximum value of the total correlation coefficient matrix, and the advection wind speed vector value is calculated through the time interval between adjacent moments of the radar base data and the change of the sector grid position, so as to complete the inversion of the advection wind vector field.

6. The typhoon wind field inversion method according to claim 5, characterized in that: The specific steps of performing quality control on the inverted advection wind vector field and filling in the wind field difference and missing values ​​in the rectangular coordinate system include: Inversion wind quality control: Determine the rectangular coordinate grid with the typhoon center as the origin for the final output of the inverted advection wind vector, average the inverted advection wind speed within the search neighborhood radius of each grid point, and if the absolute value of the difference between a certain inverted advection wind speed and the average value within the search neighborhood radius exceeds the threshold, the inverted advection wind vector is discarded; Cartesian coordinate system wind field interpolation: For each grid point, if the number of valid advection wind vectors within its search neighborhood radius is not less than the preset minimum value and is distributed in at least two quadrants, the advection wind vector of the grid point is interpolated using the least squares method; Filling in missing values: Traverse the rectangular coordinate grid points and screen out the missing grid points. For each missing grid point, determine its search neighborhood radius. When the number of valid advection wind vectors within the search range is not less than the preset minimum value, use least squares interpolation to fill in the missing values.

7. A typhoon wind field inversion system, characterized in that: It includes pre-processing module, inversion module and smoothing correction module; The preprocessing module is used to preprocess the radar base data to obtain polar coordinate data with the typhoon center as the origin and divided into multiple sector grids and distributed with a number of data points; The inversion module is used to obtain the corresponding sector grids with correlation between the sector grids in the typhoon search area at adjacent times based on the echo intensity changes at adjacent times, and then perform the inversion of the advection wind vector field; The smoothing correction module is used to perform iterative smoothing correction on the inverted advection wind vector field, remove data points with incomplete velocity deblurring, and output the typhoon wind field; The specific steps of iteratively smoothing and correcting the inverted advection wind vector field, eliminating data points with incomplete velocity deblurring, and outputting the typhoon wind field include: S31. Perform quality control on the inverted advection wind vector field, and fill in the wind field difference and missing values ​​in the rectangular coordinate system; S32, the inverted advection wind vector field (u 0 ,v 0 ) is projected onto the radar radial direction to obtain the inverted radial wind component In the formula, dx and dy represent (u 0 ,v 0 ) is located at a distance from the center of the radar in the x-axis and y-axis directions, and r is (u 0 ,v 0 ) is located at a grid point and the radar center; S33, the inverted radial wind component Radar observation of radial wind V r For comparison, if Then iterative correction is performed; otherwise, the output (u 0 ,v 0 ); The specific steps of iterative smoothing correction include: S331, if The wind direction Dir of the inverted advection wind vector is correct, and the wind speed Spd 0 Weak, set the number of iterations m = 1; S332, update the wind speed with a fixed step size Δs, and keep the wind direction Dir unchanged. The updated wind speed is recorded as Spd m , namely Spd m =Spd 0 +m×Δs, Calculate the updated wind speed component u m 、v m , the calculation formula is as follows: S333、(u m ,v m ) is projected onto the radar radial direction to obtain Recalculate the feature quantity m is the number of iterations, S334, when the feature quantity F<ε or m>β, the iteration is stopped, ε is the convergence threshold, β is the iteration number threshold; if the iteration is stopped when m>β, the u that minimizes F is selected m 、v m Replace u 0 、v 0 , and the final output of the typhoon three-dimensional wind field is obtained; if the iteration stops when f<ε, the u at the time of iteration stop is output m 、v m Replace u 0 、v 0 , and obtain the final output typhoon three-dimensional wind field; otherwise, set m=m+1 and return to step S332.

8. An electronic device, characterized in that: The electronic device comprises a memory and a processor, wherein the memory stores a computer program, and the processor is configured to call and run the computer program stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 6 are implemented.

Images