CN109459753A - Weather radar data coordinate converts Fast Interpolation method - Google Patents

Abstract

The invention provides a fast interpolation method for coordinate conversion of weather radar data. The sinc interpolation kernel table is pre-calculated according to the preset number of interpolation kernel points and the quantized displacement; each network in the uniform three-dimensional grid of the geographic coordinate system to be converted is calculated. The grid point coordinate value of the point corresponds to the coordinates in the polar coordinate system centered on the current radar itself, and calculates its position in the uniform polar coordinate discrete grid; according to the integer part of the position value, from the given Extract a three-dimensional matrix data block from the weather radar volume scan data saved by discretization; query from the sinc interpolation kernel table according to the fractional part of the position value and obtain three rows of elements in the table to form column vectors respectively and Using 3D matrix data blocks and Get a two-dimensional matrix data block; use the two-dimensional matrix data block and Get a one-dimensional column vector; use the one-dimensional column vector and Get the interpolation result of the current grid value. The present invention improves the interpolation speed.

Description

A Fast Interpolation Method for Coordinate Transformation of Weather Radar Data

technical field

The invention relates to the field of weather radar signal and data processing, in particular to a fast interpolation method for coordinate transformation of weather radar data.

Background technique

The data products obtained by the weather radar are usually in polar coordinates centered on the radar itself, including the distance, azimuth and pitch coordinates of the meteorological target. For the application of networked observation of multiple weather radars, a unified coordinate system needs to be constructed. For example, the geographic coordinate system composed of longitude, latitude and height is used as the benchmark. The data products obtained by each radar are transformed from their own polar coordinate system to Under the new coordinate system, it is convenient for the networking puzzle of weather radar data. Since most meteorological objects, such as clouds, are spatially continuous and have many fine-scale structures, it is hoped that the interpolated scatter rate field should be spatially continuous, and the radar should be preserved as much as possible during the interpolation process. The original echo structural features present in the data. Common interpolation methods include: nearest neighbor method, linear interpolation method and sinc kernel interpolation method. Among them, although the nearest neighbor method is fast, its accuracy is too low, which easily leads to the spatial discontinuity of the interpolated radar scattering rate; the linear interpolation method introduces a large smoothing to the original data before interpolation, so it is easy to lose the original scattering rate Fine-scale structural features. The sinc kernel interpolation method is a method of interpolation according to the Nyquist sampling theorem, which has the highest accuracy, but usually has a high computational complexity.

SUMMARY OF THE INVENTION

In view of the above technical problems, the present invention provides a fast interpolation method for coordinate transformation of weather radar data.

The fast interpolation method for coordinate conversion of weather radar data in the present invention includes: step 1, pre-calculating a sinc interpolation kernel table according to a preset number of interpolation kernel points and quantized displacement; step 2, calculating a geographic coordinate system to be converted to through a mapping relationship The coordinates in the polar coordinate system with the current radar itself as the center corresponding to the grid point coordinate value of each network point in the uniform three-dimensional grid; step 3, for a specific distance, pitch calculated in step 2 and azimuth coordinates, calculate its position in the uniform polar coordinate discrete grid of the weather radar volume scan data; step 4, according to the integer part of the position value of the position obtained in step 3, save from the given discretization A three-dimensional matrix data block is extracted from the weather radar volume scanning data of ; Step 5, according to the fractional part of the position value of the position obtained in step 3, query from the sinc interpolation kernel table and obtain three rows of elements in the table to respectively make up a column vector and Step 6, use the three-dimensional matrix data block obtained in step 4 and the column vector obtained in step 5 Carry out weighted operation to obtain a two-dimensional matrix data block; Step 7, utilize the two-dimensional matrix data block obtained in step 6 and the column vector obtained in step 5 Perform a weighted operation to obtain a one-dimensional column vector; step 8, use the one-dimensional column vector obtained in step 7, and the column vector obtained in step 5 Perform a weighting operation to obtain the interpolation result of the current grid value.

Preferably, in the step 1, the typical value of the interpolation kernel number P is any even number between 6 and 16, and the typical value of L in the quantization displacement 1/L is an even number greater than or equal to 10;

The sinc interpolation kernel table is a numerical table of row L+1 and column P;

The value of the element w _i,j of the i-th row and the j-th column of the sinc interpolation kernel table is calculated by the following formula

Among them, i=1,2,...L+1, j=1,2,...P, sin c(x)=sin(πx)/(πx) represents the sinc function.

Preferably, in the step 2, the mapping from the grid point coordinate values (x _lat,m' ,y _lon,n' ,h _k' ) of the geographic coordinate system to the coordinates of the polar coordinate system centered on the current radar itself The expression for the relationship is:

Among them, (x _r , y _r , hr ₎ is the longitude, latitude and altitude coordinates of the weather radar itself, R represents the earth’s radius, s=R×arccos[sin(x _lat,m′ )sin(x _r )+cos(x _lat,m′ )cos(x _r )cos(y _lon,n′ -y _r )]x _lat,m′ ,y _lon,n′ ,h _{k′represent} the m′th latitude, The n'th longitude and the k'th altitude coordinate.

Preferably, in the step 3, a specific distance, pitch and azimuth coordinates Calculate its uniform polar discrete grid in weather radar volume scan data

The method for the position (x1,x2,x3) in the

Preferably, the step 4 includes:

Step 41, round up x1, x2, and x3 respectively to get in represents the rounding operator;

Step S42, the weather radar volume scan data saved in the given discretization , where m=1,2,...M, n=1,2,...N, k=1,2,...K, indexed by the first dimension from n ₁ -(P/ 2-1), n ₁ -P/2, n ₁ -P/2+1, ..., n ₁ +P/2, the second dimension is indexed from n ₂ -(P/2-1), n ₂ -P/2, n ₂ -P/2+1, ..., n ₂ +P/2, the third dimension is indexed from n ₃ -(P/2-1), n ₃ -P/2, n ₃ -P/2+1,...,n ₃ +P/2, extract a P×P×P dimensional three-dimensional matrix data block s(i,j,l), i,j,l=1,2,. ..P.

Preferably, the step 5 includes:

Calculate the fractional part of x1 Divide it by the value of the quantization shift 1/L and round to an integer in Indicates the rounding operator, queries the sinc interpolation table and selects the elements of the L+1-m1 row in the interpolation table as the weighted value to form a row vector round() indicates the rounding operator;

Calculate the fractional part of x2 Divide it by the value of the quantization shift 1/L and round to an integer Query the sinc interpolation table and select the elements of the L+1-m2 row in the interpolation table as the weighted value to form a row vector

Calculate the fractional part of x3 Divide it by the value of the quantization shift 1/L and round to an integer Query the sinc interpolation table and select the elements of the L+1-m3 row in the interpolation table as the weighted value to form a row vector

Preferably, the step 6 includes:

In the three-dimensional matrix data block s _P×P×P , for a fixed index pair (j,l), the P data elements s(1,j,l) of the three-dimensional matrix data block s _P×P×P , s(2,j,l),...,s(P,j,l) form a column vector Then find the vector and the column vector inner product of As a data element of the jth row and the lth column of a P×P two-dimensional matrix data block s′ _P×P×P , the superscript T represents the matrix or vector transposition. All index pairs (j,l), j,l=1,2,...P are traversed until a P×P two-dimensional matrix data block s′ _P×P×P is obtained by calculation.

Preferably, the expression of the one-dimensional column vector in the step 7 is:

in, is a one-dimensional column vector, s′ _P×P×P is a two-dimensional matrix data block, is a column vector, and T represents the matrix or vector transpose.

Preferably, in the step 8, the interpolation result of the current grid value is:

where T represents matrix or vector transpose.

As can be seen from the above technical solutions, the present invention has the following beneficial effects:

(1) In the technical solution of the present invention, the coordinate transformation of weather radar data uses sinc kernel interpolation, and the interpolation accuracy is high, which is helpful to retain the fine-scale structural characteristics of the scattering rate of meteorological targets;

(2) The implementation of the interpolation in the technical solution of the present invention is to obtain the weighted value by querying the pre-calculated sinc interpolation kernel table, which avoids the process of calculating the sinc function with high complexity in real time, and improves the interpolation speed.

Description of drawings

1 is a schematic diagram of volume scanning performed by a weather radar according to an embodiment of the present invention;

FIG. 2 is a schematic diagram of a sinc kernel interpolation principle used in an embodiment of the present invention.

FIG. 3 is a flowchart of a fast interpolation method for coordinate transformation of weather radar data according to an embodiment of the present invention.

Detailed ways

In order to make those skilled in the art better understand the technical solutions of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

In the technical solution of the present invention, sinc kernel interpolation is used for the coordinate transformation of weather radar data, and the interpolation accuracy is high, which is helpful for retaining the fine-scale structural features of the scattering rate of meteorological targets; the interpolation is realized by querying the pre-calculated sinc interpolation kernel table to obtain the weighting value, avoids the process of real-time calculation of the sinc function with high complexity, and improves the interpolation speed.

The embodiments of the present invention are carried out for the transformation of weather radar data from polar coordinate data to other coordinate systems such as geographic coordinate systems represented by latitude, longitude and altitude. Before introducing the details of the embodiments of the present invention in detail, the principle of realizing coordinate transformation by data interpolation is briefly introduced first.

Please refer to Figure 1. Usually, the weather radar performs volume scanning to obtain the distance r and azimuth angle of the meteorological target in the polar coordinate system. The basic data product of the radar scatter factor of the pitch angle θ coordinate, let it be expressed as Assuming that the data after the coordinate transformation of the data to the geographic coordinate system should be Z'(x _lat , y _lon , h), where x _lat represents the longitude coordinate, y _lon represents the dimensional coordinate, and h represents the height coordinate. For a meteorological target at a specific coordinate value (x _lat , y _lon , h), the coordinate-transformed data Z′ (x _lat , y _lon , h) is given by Obtained by coordinate mapping:

Among them, r(x _lat ,y _lon ,h),θ(x _lat ,y _lon ,h), It can be obtained through the theory of great circle geometry. If the latitude and longitude coordinates of the weather radar's own position are (x _r , y _r , hr ₎ , then

Among them, R is the radius of the earth, and the expression of s is

s=R×arccos[sin(x _lat )sin(x _r )+cos(x _lat )cos(x _r )cos(y _lat -y _r )]

Typically, data obtained from weather radar volume scans It is uniformly grid distributed in the polar coordinate system, and the data stored in discretization , m=1,2,...M, n=1,2,...N, k=1,2,...K. M, N, and K represent the data dimensions of the weather radar data in the range, pitch, and azimuth directions in the polar coordinate system before interpolation, respectively. Polar coordinates r(x _lat ,y _lon ,h),θ(x _{lat ,} y _{lon ,} h), It is not necessarily exactly on the integer grid point on the polar coordinate system, and the transformed data Z'(x _lat , y _lon , h) needs to be obtained by interpolation method using the data on the surrounding grid points.

Common interpolation methods include: nearest neighbor method, linear interpolation method and sinc kernel interpolation method. Among them, although the nearest neighbor method is fast, its accuracy is too low, which easily leads to the spatial discontinuity of the interpolated radar scattering rate; the linear interpolation method introduces a large smoothing to the original data before interpolation, so it is easy to lose the original scattering rate Fine-scale structural features. The sinc kernel interpolation method is a method of interpolation according to the Nyquist sampling theorem, which has the highest accuracy, but usually has a high computational complexity.

In the embodiment of the present invention, the interpolation is performed based on sinc kernel interpolation. Before introducing the details of the embodiments of the present invention in detail, it is necessary to briefly describe some concepts and principles of sinc kernel interpolation.

Generally, on the premise of satisfying the signal band limit and the sampling rate is greater than the Nyquist sampling frequency, the discrete sequence s[i], i=1,2,... can be obtained by sinc kernel interpolation to obtain any non-integer sampling point x The high-precision interpolation result s'(x) is expressed as follows:

Among them, sin c(x)=sin(πx)/(πx) is called the sinc interpolation kernel. The above interpolation formula indicates that the signal value at any non-integer sampling point x can be obtained by the weighted sum of the discrete sequence s[i], The weight of each discrete sequence sample point is the sinc interpolation kernel.

FIG. 2 is a schematic diagram of the sinc kernel interpolation principle used in the embodiment of the present invention. Generally, it is not necessary to calculate the sinc kernel weights for all sample points, but a discrete sequence s of 8 to 16 points around the non-integer sample point x is used. The sample value of [i] is involved in the interpolation operation. For example, in Figure 2, if the value at the non-integer sample point x needs to be interpolated, then four integer sample point values on the left and right of it are used (the four “o” symbols on the left and right of x represent the The sample value of ) and the corresponding weight of the sinc function (the sample value of the sinc function represented by "□") are weighted and summed. When calculating the sinc kernel interpolation formula, since the process of calculating the value of the interpolation kernel sinc(x) according to the input x value involves trigonometric function operations, the computational complexity is high. In order to overcome this shortcoming, the patent of the present invention will pre-calculate The interpolation kernels are stored in a table format to avoid real-time calculations, thereby increasing the interpolation speed.

FIG. 1 is a flowchart of a fast interpolation method for coordinate transformation of weather radar data according to an embodiment of the present invention. Referring to FIG. 1, in one embodiment of the present invention, a fast interpolation method for coordinate transformation of weather radar data is provided, and the method may include:

Step S1, pre-calculate the sinc interpolation kernel table according to the preset number of interpolation kernel points P and the quantization displacement 1/L, as shown in Table 1;

The typical value of the number of interpolation kernel points P is any even number between 6 and 16, and the typical value of L in the quantization displacement 1/L is an even number greater than or equal to 10;

The sinc interpolation kernel table is a numerical table of row L+1 and column P;

The value of the element w _i,j of the i-th row and the j-th column of the sinc interpolation kernel table is calculated by the following formula

Among them, i=1,2,...L+1, j=1,2,...P, sinc(x)=sin(πx)/(πx) represents the sinc function.

Table 1 is the interpolation kernel table of 13 rows and 8 columns when P=8, L=12

0 0 0 0 1 0 0 0 -0.021 0.028 -0.043 0.090 0.990 -0.076 0.040 -0.027 -0.042 0.057 -0.087 0.192 0.961 -0.137 0.074 -0.051 -0.061 0.083 -0.130 0.304 0.912 -0.182 0.101 -0.070 -0.077 0.105 -0.169 0.422 0.843 -0.211 0.120 -0.084 -0.088 0.122 -0.199 0.540 0.756 -0.222 0.130 -0.092 -0.093 0.131 -0.218 0.653 0.653 -0.218 0.131 -0.093 -0.092 0.130 -0.222 0.756 0.540 -0.199 0.122 -0.088 -0.084 0.120 -0.211 0.843 0.422 -0.169 0.105 -0.077 -0.0717-1059 0.101 -0.182 0.912 0.304 -0.130 0.083 -0.061 -0.051 0.074 -0.137 0.961 0.192 -0.087 0.057 -0.042 -0.027 0.040 -0.076 0.990 0.090 -0.043 0.028 -0.021 0 0 0 1 0 0 0 0

Step S2, set the uniform three-dimensional grid of the geographic coordinate system to be converted to Ω={(x _lat,m′ ,y _lon,n′ ,h _k′ )|m′=1,2,...M′ , n′=1,2,...N′, k′=1,2,...K′}, where x _lat,m′ ,y _lon,n′ ,h _{k′represent} the grid represented by The m'th latitude, n'th longitude and k'th height coordinates of , M, N and K respectively represent the size dimension of the grid, for each specific grid point coordinate value of all grid points (x _lat,m′ ,y _lon,n′ ,h _k′ ), first calculate the corresponding coordinates in the polar coordinate system centered on the current radar itself through the mapping relationship f Then, the following steps S3 to S8 are repeatedly executed, wherein r′, θ′, Represent the distance, pitch and azimuth coordinates in the polar coordinate system of the mapping;

The expression of the mapping relationship f from the grid point coordinate values of the geographic coordinate system (x _{lat, m′} , y _{lon, n′} , h _k′ ) to the coordinates of the polar coordinate system centered on the current radar itself is: :

Among them, (x _r , y _r , hr ₎ is the latitude, longitude and altitude coordinates of the weather radar itself, R is the radius of the earth, and the expression of s is

s=R×arccos[sin(x _lat,m′ )sin(x _r )+cos(x _lat,m′ )cos(x _r )cos(y _lon,n′ -y _r )]

Step S3, for a specific distance, pitch and azimuth coordinates calculated in step S2 Calculate its uniform polar discrete grid in weather radar volume scan data The position (x1, x2, x3) in , which indicates that the distance coordinate r' is located at the x1th sample point in the rm sequence, _m =1, 2,...M; the pitch coordinate θ' is located in the θ _n sequence The x2th sample point in , n=1,2,...N; azimuth coordinates lie in The x3th sample point in the sequence, k=1,2,...K; x1,x2,x3 are integers or non-integers not less than 1;

described by a specific distance, pitch and azimuth coordinates Calculate its uniform polar discrete grid in weather radar volume scan data The method for the position (x1, x2, x3) in the :

Step S4, according to the integer part of the x1, x2, x3 position values obtained in step S3, from the given discretized saved weather radar volume scan data m=1,2,...M, n=1,2,...N, k=1,2,...K, extract a P×P×P dimensional three-dimensional matrix data block s( i,j,l), i,j,l=1,2,...P, wherein the integer P is the preset number of interpolation kernel points P shown in the step s1;

Said weather radar volume scan data saved from a given discretization according to the integer parts of the x1, x2, x3 position values m=1,2,...M, n=1,2,...N, k=1,2,...K, extract a P×P×P dimensional three-dimensional matrix data block s( i,j,l), i,j,l=1,2,...P The steps can include:

Sub-step S41: respectively round x1, x2, and x3 to obtain in represents the rounding operator;

Sub-step S42: Weather radar volume scan data saved for a given discretization In m=1,2,...M, n=1,2,...N, k=1,2,...K, according to the index of the first dimension from n ₁ -(P/2- 1), n ₁ -P/2, n ₁ -P/2+1, ..., n ₁ +P/2, the second dimension is indexed from n ₂ -(P/2-1), n ₂ -P /2, n ₂ -P/2+1, ..., n ₂ +P/2, the third dimension is indexed from n ₃ -(P/2-1), n ₃ -P/2, n ₃ -P /2+1,...,n ₃ +P/2, extract a P×P×P dimensional three-dimensional matrix data block s(i,j,l), i,j,l=1,2,... P;

Step S5: According to the fractional parts of the x1, x2, and x3 position values obtained in step S3, query from the sinc interpolation kernel table and obtain three rows of elements in the table, which form column vectors respectively and

According to the fractional part of the position value of x1, x2, x3, query from the sinc interpolation kernel table and obtain three rows of elements in the table, which form column vectors respectively and The method is:

Calculate the fractional part of x1 Divide it by the value of the quantization shift 1/L and round to an integer in Indicates the rounding operator, queries the sinc interpolation table and selects the elements of the L+1-m1 row in the interpolation table as the weighted value to form a row vector round() indicates the rounding operator;

Calculate the fractional part of x2 Divide it by the value of the quantization shift 1/L and round to an integer Query the sinc interpolation table and select the elements of the L+1-m2 row in the interpolation table as the weighted value to form a row vector

Calculate the fractional part of x3 Divide it by the value of the quantization shift 1/L and round to an integer Query the sinc interpolation table and select the elements of the L+1-m3 row in the interpolation table as the weighted value to form a row vector

Step S6, use the three-dimensional matrix data block s _P×P×P obtained in step S4 and the column vector obtained in step S5 Perform a weighted operation to obtain a P×P two-dimensional matrix data block s′ _P×P×P ;

The use of three-dimensional matrix data blocks s _P×P×P and column vectors The process of performing the weighting operation is as follows:

In the three-dimensional matrix data block s _P×P×P , for a fixed index pair (j,l), the P data elements s(1,j,l) of the three-dimensional matrix data block s _P×P×P , s(2,j,l),...,s(P,j,l) form a column vector Then find the vector and the column vector inner product of As a data element of the jth row and the lth column of a P×P two-dimensional matrix data block s′ _P×P×P , the superscript T represents the matrix or vector transposition. Traverse all index pairs (j,l), j,l=1,2,...P until the P×P two-dimensional matrix data block s′ _P×P×P is obtained by calculation;

Step S7, use the two-dimensional matrix data block s′ _P×P×P obtained in step S6 and the column vector obtained in step S5 Perform a weighted operation to get a one-dimensional column vector containing P elements

The described use of two-dimensional matrix data blocks s′ _P×P×P and column vectors The weighting operation is performed using the multiplication of vectors and matrices, and its expression is:

get a 1D column vector with P elements Among them, the superscript T represents the matrix or vector transpose;

Step S8, use the one-dimensional vector containing P elements obtained in step S7 with the column vector obtained in step S5 Perform a weighted operation to get the interpolation result of the current grid value, the expression is

where the superscript T represents the matrix or vector transpose.

The above embodiments are only exemplary embodiments of the present invention, and are not intended to limit the present invention, and the protection scope of the present invention is defined by the claims. Those skilled in the art can make various modifications or equivalent replacements to the present invention within the spirit and protection scope of the present invention, and such modifications or equivalent replacements should also be regarded as falling within the protection scope of the present invention.

Claims (9)

1. A weather radar data coordinate conversion fast interpolation method is characterized by comprising the following steps:

step 1, pre-calculating a sinc interpolation kernel table according to a preset interpolation kernel number and a quantitative displacement;

step 2, calculating coordinates under a polar coordinate system taking the current radar as the center, which correspond to grid point coordinate values of each network point in a uniform three-dimensional grid of the geographic coordinate system to be converted through the mapping relation;

step 3, calculating the position of a certain specific distance, pitch and azimuth coordinates obtained by calculation in the step 2 in a uniform polar coordinate discrete grid of the volume scanning data of the weather radar;

step 4, extracting a three-dimensional matrix data block from the given discretized stored weather radar volume scanning data according to the integer part of the position value of the position obtained in the step 3;

step 5, respectively inquiring and obtaining three rows of elements in the table from the sinc interpolation kernel table according to the decimal part of the position value of the position obtained in the step 3 to respectively form column vectors And

step 6, utilizing the three-dimensional matrix data block obtained in the step 4 and the column vector obtained in the step 5Performing weighting operation to obtain a two-dimensional matrix data block;

step 7, using the two-dimensional matrix data block obtained in step 6 and the column vector obtained in step 5Performing weighting operation to obtain a one-dimensional column vector;

step 8, using the one-dimensional column vector obtained in step 7 and the column vector obtained in step 5And performing weighting operation to obtain an interpolation result of the current grid value.

2. The method according to claim 1, wherein in the step 1, a typical value of the number P of the interpolation kernels is any even number between 6 and 16, and a typical value of L in quantization displacement 1/L is an even number greater than or equal to 10;

the sinc interpolation kernel table is a numerical table with L +1 rows and P columns;

element w of ith row and jth column of the sinc interpolation kernel table_i,jIs calculated by the following formula

Where i ═ 1,2,. L +1, j ═ 1,2,. P, sinc (x) ═ sin (pi x)/(pi x) denotes a sinc function.

3. Method according to claim 2, characterized in that in step 2 the grid point coordinate values (x) of the geographical coordinate system are determined from the grid point coordinate values (x) of the geographical coordinate system_lat,m′,y_lon,n′,h_k′) The expression of the mapping relationship to the polar coordinate system coordinate centered on the current radar itself is:

wherein (x)_r,y_r,h_r) Is the longitude and latitude height coordinate of the weather radar self position, R represents the earth radius, s is R multiplied by arccos [ sin (x)_lat,m′)sin(x_r)+cos(x_lat,m′)cos(x_r)cos(y_lon,n′-y_r)]x_lat,m′,yl_on,n′,h_k′Respectively, the mth 'latitude, the nth' longitude and the kth altitude coordinate represented by the grid.

4. The method of claim 3, wherein in step 3, the distance, pitch and azimuth are specified by a certain valueCoordinates of the objectComputing its uniform polar discrete grid in weather radar volume scan dataThe method of the position (x1, x2, x3) is

5. The method of claim 4, wherein the step 4 comprises:

step 41, rounding x1, x2 and x3 respectively to obtain WhereinRepresenting a rounding operator;

step S42, storing weather radar volume scan data at given discretizationWhere M is 1,2,. M, N is 1,2,. N, K is 1,2,. K, K is indexed in a first dimension from N₁-(P/2-1)、n₁-P/2、n₁-P/2+1、…、n₁+ P/2, index of the second dimension from n₂-(P/2-1)、n₂-P/2、n₂-P/2+1、…、n₂+ P/2, index of the third dimension from n₃-(P/2-1)、n₃-P/2、n₃-P/2+1、…、n₃+ P/2, a three-dimensional matrix data block s (i, j, l) of dimensions P × P is extracted, i, j, l being 1, 2.

6. The method of claim 5, wherein the step 5 comprises:

calculate the fractional part of x1Divide it by the value of the quantized displacement 1/L and round it down to an integerWhereinExpressing an operator, inquiring a sinc interpolation table, selecting the elements of the L +1-m1 th lines in the interpolation table as weighted values to form a line vectorround () represents the rounding operator;

calculate the fractional part of x2Divide it by the value of the quantized displacement 1/L and round it down to an integerInquiring the sinc interpolation table and selecting the elements of the L +1-m2 th line in the interpolation table as the weighted value to form a line vector

Calculate the fractional part of x3Divide it by the value of the quantized displacement 1/L and round it down to an integerInquiring the sinc interpolation table and selecting the elements of the L +1-m3 th line in the interpolation table as the weighted value to form a line vector

7. The method of claim 6, wherein the step 6 comprises:

in a three-dimensional matrix data block s_P×P×PFor a fixed index pair (j, l), the three-dimensional matrix data block s_P×P×PP data elements s (1, j, l), s (2, j, l),.. s (P, j, l) form a column vectorThen, the vector and the column vector are obtainedInner product of (2)As a P × P two-dimensional matrix data block s'_P×P×PThe jth row of (a), the ith column of data elements, where superscript T represents a matrix or vector transpose; p is traversed for all index pairs (j, l), j, l ═ 1, 2.. until a P × P two-dimensional matrix data block s 'is computed'_P×P×P。

8. The method of claim 7, wherein the expression of the one-dimensional column vector in step 7 is:

wherein,is one dimensional column vector, s'_P×P×PIn the form of a two-dimensional matrix of data blocks,for a column vector, T represents a matrix or vector transpose.

9. The method according to claim 8, wherein in step 8, the interpolation result of the current grid value is:

where T represents a matrix or vector transpose.