CN111060887B - Gm-APD laser radar low signal-to-noise ratio echo data signal extraction method based on concave-convex search - Google Patents

Abstract

The invention discloses a method for extracting low signal-to-noise ratio echo data signals of a Gm-APD laser radar based on concave-convex search. Step 1: performing preprocessing by using convolution of a Gaussian function and the Gaussian function so as to remove abnormal peaks and obtain a smooth distribution histogram; step 2: extracting the features of the target on the smooth distribution histogram, solving first and second derivatives of the smooth distribution histogram, and determining the distribution of maximum points under the current variance; and step 3: and combining the distance value results of the first-order search and the second-order search, and judging and reserving a correct target distance value by referring to the distance value of the cross neighborhood pixels. The invention is used for signal extraction of echo data with low peak signal-to-noise ratio, and can realize the work of extracting a long-distance target detection signal.

Description

Gm-APD laser radar low signal-to-noise ratio echo data signal extraction method based on concave-convex search

Technical Field

The invention belongs to the technical field; in particular to a method for extracting low signal-to-noise ratio echo data signals of a Gm-APD laser radar based on concave-convex search.

Background

At present, Gm-APD laser imaging radar signal extraction methods mainly comprise two types, namely nonparametric estimation and parametric estimation. The nonparametric estimation is used for extracting echo characteristics, wherein a peak value method considers that the trigger frequency of a target echo is higher than noise, the peak value of a trigger frequency histogram is extracted to be a target position for reconstruction, the method can better extract the position of the target echo under the conditions of low noise, high signal-to-noise ratio and multi-frame number, and the reconstruction effect on echo data with high noise data and few frames is poor. The time correlation method reconstructs signals in a pulse width according to the correlation, and the signals are considered as a target echo position by sequencing and judging whether the triggering interval of two times is smaller than the pulse width. In general, the reconstruction effect of the existing non-parameter estimation method on high-noise data is poor. The parameter estimation utilizes a Poisson distribution trigger model of Gm-APD to reconstruct, more frames are needed during reconstruction, wherein the maximum likelihood estimation utilizes the trigger model and a likelihood function to search in a parameter space to obtain a distance corresponding to the maximum value of the likelihood function to reconstruct, and when the actual trigger histogram is different from the trigger model, a larger error is caused. Bayesian estimation utilizes reversible jump Monte Carlo Markov to combine with trigger model to obtain the target position of multimodal echo data, depending on the trigger model. And a method for reversely pushing and correcting the trigger frequency histogram reduces the trigger frequency of the background noise at the opening end and highlights the peak value part of the target. In addition, the spatial correlation processing method can be used for improving the statistical frame number by combining with other methods, the reconstruction effect can be improved under the condition of less frame number or extremely low signal-to-noise ratio, and the target edge blurring phenomenon can be caused when the signal-to-noise ratio is slightly high.

Disclosure of Invention

The invention provides a method for extracting a Gm-APD laser radar low signal-to-noise ratio echo data signal based on concave-convex search, which combines a first-order search result and a second-order search result by utilizing the convex characteristic of an echo signal; and when the combination is carried out, the neighborhood information is used for carrying out judgment for many times, the target pixel point is reserved, and the noise pixel point is removed.

The invention is realized by the following technical scheme:

a method for extracting a Gm-APD laser radar low signal-to-noise ratio echo data signal based on concave-convex search comprises the following steps:

step 1: performing preprocessing by using convolution of a Gaussian function and the Gaussian function so as to remove abnormal peaks and obtain a smooth distribution histogram;

step 2: extracting the features of the target on the smooth distribution histogram, solving first and second derivatives of the smooth distribution histogram, and determining the distribution of maximum points under the current variance;

and step 3: and combining the distance value results of the first-order search and the second-order search, and judging and reserving a correct target distance value by referring to the distance value of the cross neighborhood pixels.

Further, the preprocessing expression of step 1 is as follows:

wherein v is a kernel density function, i.e., a gaussian smoothing function, h is a variance value, u is a trigger histogram, and w is a smoothed histogram.

Further, the step 2 first derivative expression is:

the second derivative expression is:

wherein v is₁Is a variance of h₁Kernel density function of v₂Is a variance of h₂Kernel density function of (1), w₁And w₂Corresponding to the first and second derivatives of the preprocessed data;

the first derivative target location is characterized by w₁In the two adjacent values, the left value is larger than zero, the right value is smaller than zero, the position closest to the middle zero point in the two values is a set of target possible values, and finally the position with the largest difference between the left nearest maximum value and the right nearest minimum value in the possible values is extracted; second derivative position characterized by w₂Position corresponding to the minimum value。

Further, the expression of extracting the signal in the first derivative and the second derivative is as follows:

θ＝(i,i+1)|w₁(i)＞0,w₁(i+1)＜0

pos_1stpro＝argmin|w₁(θ)|

θ₁＝sort(maximum(w₁),pos_1stpro)

θ₂＝sort(minimum(w₁),pos_1stpro)

pos_2nd＝argmin|w₂(i)|

where θ is the possible value of the first characteristic satisfying the first derivative, pos_1stproTo satisfy the first and second characteristic possible values of the first derivative, pos_1stFor first order search results, pos_2ndIs a second order search result.

Further, the specific process in step 3 is as follows:

step 3.1: keeping and judging the value of the target pixel point;

step 3.2: and removing the noise pixel points and supplementing the target pixel points.

Further, the step 3.1 specifically comprises:

step 3.1.1: calculating the difference value between the central pixel point and the cross neighborhood, and sequencing the difference values from small to large to obtain D through first-order search and second-order search respectively_1stAnd D_2stThe cross neighborhood difference sorting correspondence formula of the two sequences is as follows:

D_1st＝sort(D_{value_1st})＝d₁,d₂,d₃,d₄

D_2nd＝sort(D_{value_2nd})＝d₅,d₆,d₇,d₈

step 3.1.2: respectively judge the difference D_1st，D_2ndWhether the number of the medium valued points is more than or equal to 2 or not is determined as follows:

(isnan(D_1st||D_2nd)≤2)

if D is_1stThe number of the valued points is more than or equal to 2, S₁Searching for a difference D to a first order_1stAverage value of the first two values in the sequence, if not, S₁Searching for a difference D to a first order_1stA first value of the sequence; s₂The same process is carried out;

step 3.1.3: judgment of | P_1st(i,j)-P_2nd(i, j) | and threshold A relation judge S at the same time₁And S₂In relation to a threshold A;

when | P_1st(i,j)-P_2nd(i, j) | is less than threshold A or S₁And S₂Meanwhile, when the distance is smaller than the threshold A, the obtained distance is considered as a target value, the distance value is reserved, and the judgment expression is as follows:

[P_1st(i,j)-P_2nd(i,j)＜th||(S₁＜th&&S₂＜th)]

when S is₁Less than threshold A, while S₂When the first-order search result is greater than or equal to the threshold A, the first-order search result P is considered_1st(i, j) the target, second order search result P_2nd(i, j) obtaining the wrong distance value, assigning the first-order search result of the pixel point to a second-order search result, and judging the expression as follows:

(S₁＜th&&S₂≥th)

when S is₂Less than threshold A, while S₁When the second-order search result is greater than or equal to the threshold A, the second-order search result P is considered_2nd(i, j) obtaining the first-order search result P_1st(i, j) obtaining the wrong distance value, assigning the second-order search result of the pixel point to the first-order search result, and judging the expression as follows:

(S₁≥th&&S₂＜th)

when the above is to S₁And S₂When the relation with the threshold A is not met, namely the relation is greater than the threshold, the two search results are considered to be non-target values, and are discarded and set as a nonnumber;

step 3.1.4: the routine is ended.

Further, the step 3.2 specifically includes:

step 3.2.1: calculating the number N of the current pixel array non-number points₁The calculation expression is as follows:

N₁＝isnan(P(i,j))

step 3.2.2: d obtained according to step 3.1.1_1st，D_2stS obtained according to step 3.1.2₁，S₂Simultaneously, setting the pixel points of which the cross neighborhood points are all nonnumbers or S is greater than a threshold value as nonnumbers;

step 3.2.3: judging whether three or more valued pixel points exist in the cross neighborhood, wherein the judgment expression is as follows:

(isnan(D_1st‖D_2nd)≤1)

when present, calculate D_1stVariance ss of₁And D_2stVariance ss of₂The calculation formulas are respectively as follows:

ss₁＝var(D_1st(1:3))

ss₂＝var(D_2nd(1:3))

variance ss₁The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood; variance ss₂The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood;

step 3.2.4: judging the relation between the variance ss and the threshold B, wherein the judgment expression is as follows:

(ss＜th)

when the threshold value is smaller than the variance, setting the central pixel point as the mean value of the valued pixel points with the maximum deviation; when the threshold is greater than the variance, setting it to be a nonnumber;

step 3.2.5: repeating steps 3.1.1 to 3.1.4;

step 3.2.6: calculating the number N of the non-number points of the pixel array at the moment₂The calculation formula is as follows:

N₂＝isnan(P₁(i,j))

step 3.2.7: judgment of N₁And N₂The judgment formula is as follows:

(N₁＝N₂)

when N is present₁Is not equal to N₂Then N will be₂Value is given to N₁Returning to the step 3.2.2 for circulation; when N is present₁Is equal to N₂The routine is ended.

The invention has the beneficial effects that:

the target echo peak value lower than the noise trigger frequency can also be extracted, the method does not depend on a trigger model, and the image signal to noise ratio of the three-dimensional range profile can be improved by judging for many times.

Drawings

FIG. 1 is a flow chart of the present invention.

FIG. 2 is a graph of the high noise echo data and the pre-processing result of the present invention, FIG. 2- (a) the raw data of the high noise echo, and FIG. 2- (b) the pre-processing data of the high noise echo.

FIG. 3 first derivative w of the present invention₁And second derivative w₂FIG. 3- (a) first derivative w1 and FIG. 3- (b) second derivative w 2.

FIG. 4 is a flow chart of target selection according to the present invention, FIG. 4- (a) a flow chart of a process combining first-order search and second-order search, and FIG. 4- (b) a flow chart of target pixel point determination.

FIG. 5 is a diagram of an experimental scene and a multi-frame distance map for repairing, FIG. 5- (a) a target scene and FIG. 5- (b) a distance map.

FIG. 6 shows a result diagram of a reconstructed distance image extracted from 500 frames of targets according to the present invention, FIG. 6- (a) a peak method diagram, FIG. 6- (b) a maximum likelihood estimation diagram, and FIG. 6- (c) a concave-convex search combination method diagram.

FIG. 7 is a diagram showing the variation of objective evaluation indexes of a reconstructed three-dimensional image with different frame numbers, FIG. 7- (a) a diagram showing the variation of target reduction degree, FIG. 7- (b) a diagram showing the variation of average distance measurement error, and FIG. 7- (c) a diagram showing the variation of inverse signal-to-noise ratio of the image.

Detailed Description

The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

A method for extracting a Gm-APD laser radar low signal-to-noise ratio echo data signal based on concave-convex search comprises the following steps: as shown in figure 1 of the drawings, in which,

step 1: performing preprocessing by using convolution of a Gaussian function and the Gaussian function so as to remove abnormal peaks and obtain a smooth distribution histogram; as shown in figure 2 of the drawings, in which,

step 2: extracting the features of the target on the smooth distribution histogram, solving first and second derivatives of the smooth distribution histogram, and determining the distribution of maximum points under the current variance;

and step 3: and combining the distance value results of the first-order search and the second-order search, and judging and reserving a correct target distance value by referring to the distance value of the cross neighborhood pixels.

Further, the preprocessing expression of step 1 is as follows:

wherein v is a kernel density function, i.e., a gaussian smoothing function, h is a variance value, u is a trigger histogram, and w is a smoothed histogram.

Further, the step 2 first derivative expression is:

the second derivative expression is:

wherein v is₁Is a variance of h₁Kernel density function of v₂Is a variance of h₂Kernel density function of (1), w₁And w₂Corresponding to the first and second derivatives of the preprocessed data.

The first derivative target location is characterized by w₁In the two adjacent values, the left value is larger than zero, the right value is smaller than zero, the position closest to the middle zero point in the two values is a set of target possible values, and finally the position with the largest difference between the left nearest maximum value and the right nearest minimum value in the possible values is extracted;

second derivative position characterized by w₂The position corresponding to the minimum value.

Further, the expression of extracting the signal in the first derivative and the second derivative is as follows:

θ＝(i,i+1)|w₁(i)＞0,w₁(i+1)＜0

the upper formula is w₁The left value of the two adjacent points is larger than the right value and is less than zero;

pos_1stpro＝argmin|w₁(θ)|

the point where the first derivative of the two points is closer to 0 is a possible point;

θ₁＝sort(maximum(w₁),pos_1stpro)

the above formula is that all maximum value points and possible points in the first derivative are sorted from small to large;

θ₂＝sort(minimum(w₁),pos_1stpro)

the above formula is that all minimum value points and possible points in the first derivative are ordered from small to large;

the point with the maximum difference between the left-side nearest maximum and the right-side nearest minimum of all possible points is the finally obtained first derivative search result;

pos_2nd＝argmin|w₂(i)|

the corresponding point of the minimum value of the second derivative in the above formula is a second derivative search result;

where θ is the possible value of the first characteristic satisfying the first derivative, pos_1stproTo satisfy the first and second characteristic possible values of the first derivative, pos_1stFor first order search results, pos_2ndIs a second order search result.

Further, the process distribution combining the processing and determining the target pixel points in step 3 includes the following processes:

the process distribution combining the processing and the target pixel point judgment in the step 3 comprises the following processes:

step 3.1: keeping and judging the value of the target pixel point;

in FIG. 4a), P (i, j) is the jth pixel (pixel to be determined) in the ith row, i.e., the first-order search result is P_1st(i, j), the second order search result corresponds to P_2nd(i, j), the process judges the difference value of the two search results and the relation between the difference value and the respective cross neighborhood, and the value which is considered as the target pixel point is reserved and judged.

Step 3.1.1: calculating the difference value between the central pixel point and the cross neighborhood, and sequencing the difference values from small to large to obtain D through first-order search and second-order search respectively_1stAnd D_2stThe cross neighborhood difference sorting correspondence formula of the two sequences is as follows:

D_1st＝sort(D_{value_1st})＝d₁,d₂,d₃,d₄

D_2nd＝sort(D_{value_2nd})＝d₅,d₆,d₇,d₈

step 3.1.2: respectively judge the difference D_1st，D_2stWhether the number of the medium valued points is more than or equal to 2 or not is determined as follows:

(isnan(D_1st||D_2nd)≤2)

if D is_1stThe number of the valued points is more than or equal to 2, S₁Searching for a difference D to a first order_1stAverage value of the first two values in the sequence, if not, S₁Searching for a difference D to a first order_1stA first value of the sequence; s₂The same process is carried out;

step 3.1.3: judgment of | P_1st(i,j)-P_2nd(i, j) | and threshold A relation judge S at the same time₁And S₂In relation to a threshold A;

when | P_1st(i,j)-P_2nd(i, j) | is less than threshold A or S₁And S₂Meanwhile, when the distance is smaller than the threshold A, the obtained distance is considered as a target value, the distance value is reserved, and the judgment expression is as follows:

[P_1st(i,j)-P_2nd(i,j)＜th||(S₁＜th&&S₂＜th)]

when S is₁Less than threshold A, while S₂When the first-order search result is greater than or equal to the threshold A, the first-order search result P is considered_1st(i, j) obtaining a second order search result P targeting_2nd(i, j) obtaining the wrong distance value, assigning the first-order search result of the pixel point to a second-order search result, and judging the expression as follows:

(S₁＜th&&S₂≥th)

when S is₂Less than threshold A, while S₁When the second-order search result is greater than or equal to the threshold A, the second-order search result P is considered_2nd(i, j) obtaining the first-order search result P_1st(i, j) obtaining the wrong distance value, assigning the second-order search result of the pixel point to the first-order search result, and judging the expression as follows:

(S₁≥th&&S₂＜th)

when the above is to S₁And S₂When the relation with the threshold A is not met, namely the relation is greater than the threshold, the two search results are considered to be non-target values, and are discarded and set as a nonnumber;

step 3.1.4: the routine is ended.

Step 3.2: removing noise pixel points and supplementing target pixel points;

in FIG. 4b), P (i, j) is the jth pixel (pixel to be determined) in the ith row, i.e., the first-order search result is P_1st(i, j), the second order search result corresponds to P_2nd(i, j), the process judges the relation of the respective cross neighborhoods, removes noise pixel points and supplements the targetsMarking pixel points;

step 3.2.1: calculating the number N of the current pixel array non-number points₁The calculation expression is as follows:

N₁＝isnan(P(i,j))

step 3.2.2: d obtained according to step 3.1.1_1st，D_2stS obtained according to step 3.1.2₁，S₂Simultaneously, setting the pixel points of which the cross neighborhood points are all nonnumbers or S is greater than a threshold value as nonnumbers;

step 3.2.3: judging whether three or more valued pixel points exist in the cross neighborhood, wherein the judgment expression is as follows:

(isnan(D_1st||D_2nd)≤1)

when present, calculate D_1stVariance ss of₁And D_2stVariance ss of₂The calculation formulas are respectively as follows:

ss₁＝var(D_1st(1:3))

ss₂＝var(D_2nd(1:3))

variance ss₁The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood; variance ss₂The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood;

step 3.2.4: judging the relation between the variance ss and the threshold B, wherein the judgment expression is as follows:

(ss＜th)

when the threshold value is smaller than the variance, setting the central pixel point as the mean value of the valued pixel points with the maximum deviation; when the threshold is greater than the variance, setting it to be a nonnumber;

step 3.2.5: repeating steps 3.1.1 to 3.1.4;

step 3.2.6: calculating the number N of the non-number points of the pixel array at the moment₂The calculation formula is as follows:

N₂＝isnan(P₁(i,j))

step 3.2.7: judgment of N₁And N₂The judgment formula is as follows:

(N₁＝N₂)

when N is present₁Is not equal to N₂Then N will be₂Value is given to N₁Returning to the step 3.2.2 for circulation; when N is present₁Is equal to N₂The routine is ended.

The echo data in the real scene is used for carrying out comparative analysis on the method results, and the experimental parameters are shown in table 1. The experimental scene graph and the reference range profile obtained by manually restoring part of the pixel points by combining the multi-frame experimental data reconstruction result and the photo information are shown in FIG. 5,

TABLE 1 Experimental parameters

In the case of using 500 frames of echo data, the three methods of signal extraction and reconstruction of the range profile result as shown in fig. 6.

As shown in fig. 6, the peak method and the maximum likelihood estimation method are basically unable to identify the target contour, and have more noise points; the concave-convex searching combination method has prominent target, more complete target pixel points and less noise points. In order to more clearly compare the differences of the three methods, the three methods need to be evaluated by using objective evaluation indexes, namely, target reduction degrees:

wherein d is a reconstructed distance value; d_sIs a standard distance value; d_bIs an allowable error distance value; n is the number of target total pixels, and m is the number of pixels of the allowable error distance value. K₁The value represents the proportion of the number of pixels obtained by the reconstruction of the method in the pixels occupied by the target. Second, relative average range error is used:

wherein d is_iDetermining a reconstructed distance value for the target pixel; d_siThe distance value is judged to be the standard distance value of the target pixel. K₂The value represents the average ranging error value of all target pel distance values reconstructed by the method. And finally, using the pixel signal-to-noise ratio:

wherein p is the number of pixels in the background where noise is not filtered by the threshold and which are judged as non-target in the target. K₃The value represents the ratio of the number of target pixel elements to the number of noise pixel elements which are not filtered out, wherein the target pixel elements are obtained by reconstruction in the whole pixel array. Target degree of reduction K₁The closer to 1, the better, the average range error K₂The smaller the signal to noise ratio, the better the signal to noise ratio, and the corresponding K₃The smaller the better. D is used for judging objective evaluation index_bThe results are shown in Table 2 for the various 500 frames set at 10 time intervals.

TABLE 2 average objective evaluation results for multiple groups of 500 frames

The objective evaluation result under the fixed 500 frames can be seen, the concave-convex searching method is superior to a peak value method and a maximum likelihood estimation method in signal extraction target reconstruction, the target reduction degree can be improved by 15.10% compared with the peak value method under the 500 frames of the scene, and the image signal-to-noise ratio can be improved by 21 times; compared with maximum likelihood estimation, the target reduction degree can be improved by 3.83%, and the image signal-to-noise ratio can be improved by 15 times. For more comprehensive comparison, the objective evaluation index changes at different frame numbers are analyzed, and the reconstruction effects of various methods at different frame numbers are compared, and the result is shown in fig. 7.

FIG. 7 shows that as the number of frames increases, the target reduction degree tends to be stable after increasing; the average distance measurement error tends to be stable after being reduced; the image signal-to-noise ratio reciprocal peak value method and the maximum likelihood estimation tend to be stable after being reduced, the concave-convex search combination method is slightly increased, the analysis reason is that when the combination processing is carried out, the distance values of noise points are close along with the increase of the frame number, and parts are identified as targets to be reserved when the judgment is carried out. And the three methods are longitudinally compared, the obtained concave-convex search combination method has the best reconstruction result, and the signal extraction effect of the Gm-APD can be improved.

Claims (5)

1. A Gm-APD laser radar low signal-to-noise ratio echo data signal extraction method based on concave-convex search is characterized by comprising the following steps:

step 1: performing preprocessing by using convolution of a Gaussian function and the Gaussian function so as to remove abnormal peaks and obtain a smooth distribution histogram;

step 2: extracting the features of the target on the smooth distribution histogram, solving first and second derivatives of the smooth distribution histogram, and determining the distribution of maximum points under the current variance;

and step 3: combining the distance value results of the first-order search and the second-order search, and judging and reserving a correct target distance value by referring to the distance value of the cross neighborhood pixels;

the expression of the first derivative in the step 2 is as follows:

the second derivative expression is:

wherein v is₁Is a variance of h₁Kernel density function of v₂Is a variance of h₂Kernel density function of (1), w₁And w₂Corresponding to the first and second derivatives of the preprocessed data;

the first derivative target location is characterized by w₁In the two adjacent values, the left value is larger than zero, the right value is smaller than zero, the position closest to the middle zero point in the two values is a set of target possible values, and finally the position with the largest difference between the left nearest maximum value and the right nearest minimum value in the possible values is extracted; second derivative position characterized by w₂The position corresponding to the minimum value;

the expression of signal extraction in the first derivative and the second derivative is as follows:

θ＝(i,i+1)|w₁(i)＞0,w₁(i+1)＜0

pos_1stpro＝argmin|w₁(θ)|

θ₁＝sort(maximum(w₁),pos_1stpro)

θ₂＝sort(minimum(w₁),pos_1stpro)

pos_2nd＝argmin|w₂(i)|

where θ is the possible value of the first characteristic satisfying the first derivative, pos_1stproTo satisfy the first and second characteristic possible values of the first derivative, pos_1stFor first order search results, pos_2ndIs a second order search result.

2. The method of claim 1, wherein the step 1 preprocessing expression is:

wherein v is a kernel density function, i.e., a gaussian smoothing function, h is a variance value, u is a trigger histogram, and w is a smoothed histogram.

3. The method according to claim 1, wherein the step 3 comprises the following specific steps:

step 3.1: keeping and judging the value of the target pixel point;

step 3.2: and removing the noise pixel points and supplementing the target pixel points.

4. The method according to claim 3, characterized in that said step 3.1 is in particular:

step 3.1.1: calculating the difference value between the central pixel point and the cross neighborhood, and sequencing the difference values from small to large to obtain D through first-order search and second-order search respectively_1stAnd D_2stThe cross neighborhood difference sorting correspondence formula of the two sequences is as follows:

D_1st＝sort(D_{value_1st})＝d₁,d₂,d₃,d₄

D_2nd＝sort(D_{value_2nd})＝d₅,d₆,d₇,d₈

step 3.1.2: respectively judge the difference D_1st，D_2ndWhether the number of the medium valued points is more than or equal to 2 or not is determined as follows:

(isnan(D_1st||D_2nd)≤2)

if D is_1stThe number of the valued points is more than or equal to 2, S₁Searching for a difference D to a first order_1stAverage value of the first two values in the sequence, if not, S₁Searching for a difference D to a first order_1stA first value of the sequence; s₂The same process is carried out;

step 3.1.3: judgment of | P_1st(i,j)-P_2nd(i, j) | and threshold A relation judge S at the same time₁And S₂In relation to a threshold A;

when | P_1st(i,j)-P_2nd(i, j) | is less than threshold A or S₁And S₂Meanwhile, when the distance is smaller than the threshold A, the obtained distance is considered as a target value, the distance value is reserved, and the judgment expression is as follows:

[P_1st(i,j)-P_2nd(i,j)＜th||(S₁＜th&&S₂＜th)]

when S is₁Less than threshold A, while S₂When the first-order search result is greater than or equal to the threshold A, the first-order search result P is considered_1st(i, j) the target, second order search result P_2nd(i, j) obtaining the wrong distance value, assigning the first-order search result of the pixel point to a second-order search result, and judging the expression as follows:

(S₁＜th&&S₂≥th)

when S is₂Less than threshold A, while S₁When the second-order search result is greater than or equal to the threshold A, the second-order search result P is considered_2nd(i, j) obtaining the first-order search result P_1st(i, j) obtaining the wrong distance value, assigning the second-order search result of the pixel point to the first-order search result, and judging the expression as follows:

(S₁≥th&&S₂＜th)

when the above is to S₁And S₂When the relation with the threshold A is not met, namely the relation is greater than the threshold, the two search results are considered to be non-target values, and are discarded and set as a nonnumber;

step 3.1.4: the routine is ended.

5. The method according to claim 3, characterized in that said step 3.2 is in particular:

step 3.2.1: calculating the number N of the current pixel array non-number points₁The calculation expression is as follows:

N₁＝isnan(P(i,j))

step 3.2.2: d obtained according to step 3.1.1_1st，D_2stS from step 3.1.2₁，S₂Simultaneously, setting the pixel points of which the cross neighborhood points are all nonnumbers or S is greater than a threshold value as nonnumbers;

step 3.2.3: judging whether three or more valued pixel points exist in the cross neighborhood, wherein the judgment expression is as follows:

(isnan(D_1st||D_2nd)≤1)

when present, calculate D_1stVariance ss of₁And D_2stVariance ss of₂The calculation formulas are respectively as follows:

ss₁＝var(D_1st(1:3))

ss₂＝var(D_2nd(1:3))

variance ss₁The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood; variance ss₂The meaning of (1) is that the variance value of three values is left after the maximum deviation value is removed from the four values in the cross neighborhood;

step 3.2.4: judging the relation between the variance ss and the threshold B, wherein the judgment expression is as follows:

(ss＜th)

when the threshold value is smaller than the variance, setting the central pixel point as the mean value of the valued pixel points with the maximum deviation; when the threshold is greater than the variance, setting it to be a nonnumber;

step 3.2.5: repeating steps 3.1.1 to 3.1.4;

step 3.2.6: calculating the number N of the non-number points of the pixel array at the moment₂The calculation formula is as follows:

N₂＝isnan(P₁(i,j))

step 3.2.7: judgment of N₁And N₂The judgment formula is as follows:

(N₁＝N₂)

when N is present₁Is not equal to N₂Then N will be₂Value is given to N₁Returning to the step 3.2.2 for circulation; when N is present₁Is equal to N₂The routine is ended.

Images