Underground multi-target identification method based on slope distribution estimation and stability test

Abstract

The invention discloses a novel ground penetrating radar underground pipeline and underground cavity identification technology. According to the technology, target identification is carried out by using slope characteristics of underground pipelines and cavities. Firstly, ground penetrating radar B scanning data is preprocessed, direct coupling waves and multiple reflection waves are removed, self-adaptive threshold binarization is used, and an image is divided into a plurality of sub-images; carrying out non-zero sequence sparse coding on the binarized image to convert the binarized image into a sparse matrix; then, column offset coding is used for extracting the slope distribution condition of the sparse matrix, slope information and original position information are reserved, and irrelevant items are filtered out; clustering is carried out according to the slope information and the original position information; and finally, the slope distribution characteristics of each class are judged. Experiments show that the method has a good recognition effect.

Description

Technical Field

[0001] The present invention relates to the technical field of radar image processing, and particularly relates to a method for identifying multiple underground targets based on slope distribution estimation and stationarity test. Background Art

[0002] Road pipelines are regarded as the lifelines of cities, and the main pipelines include drainage pipelines, communication cables, etc. Nowadays, the underground pipeline network is becoming increasingly intricate. If the construction is improper, it is very likely to cause large-scale water and power outages, seriously affecting the quality of residents' daily life.

[0003] As a non-destructive detection method, ground penetrating radar has advantages such as high efficiency and low cost, making it have important applications in aspects such as road surfaces, bridges, walls, tunnels, etc. In the daily maintenance of urban infrastructure, ground penetrating radar can be used to detect urban underground pipelines, cables, etc., and can accurately locate the cables when maintenance is needed; in road construction, ground penetrating radar can ensure the detection of the underground state before road construction, such as the position of tree roots and cables; in cultural relic protection, it is necessary to use ground penetrating radar to detect the state of underground cultural relics in order to avoid them for protection; in military applications, using airborne ground penetrating radar to detect landmines has the advantage of strong safety; in hazard warning, ground penetrating radar can discover these potential risks at the initial stage of the development of hazards such as wall cracks and underground cavities, improving the ability to control risks.

[0004] The interpretation of ground penetrating radar data often adopts manual interpretation, which, however, requires the experience of operators. With the increase in the amount of data, the economic and time costs relying on manual interpretation have greatly increased. In contrast, intelligent recognition methods have gradually shown great advantages. Today's intelligent recognition methods are mainly divided into two categories: traditional methods and recognition methods based on machine learning. Traditional methods include methods such as Hough transform, pattern recognition, least squares fitting, and curvelet transform; machine learning is divided into traditional machine learning (such as SVM, random forest) and deep learning, and deep learning is the mainstream research direction in recent years. Common deep learning models include SSD, Fast-RCNN, YOLO, etc. However, today's intelligent recognition methods cannot completely solve the problem of automatic recognition of underground targets by ground penetrating radar. The reasons can be attributed to: 1. There is a lot of clutter interference in ground penetrating radar data, which greatly affects the recognition accuracy of traditional methods. At the same time, the processing speed of traditional methods for relatively complex images is slow, making it difficult to meet real-time requirements. 2. The number of training data is limited, which is manifested in two aspects: on the one hand, it is impossible to find a large amount of data available for training on the network; on the other hand, the time cost of collecting underground data by ground penetrating radar is high, and the search for underground targets requires the experience of professionals. This has certain limitations on the use of machine learning methods. Although simulated data can be used as a substitute, there are differences between the simulated scenario and the actual scenario, and the recognition accuracy of the model trained with simulated data is not very satisfactory. Summary of the invention

[0005] The present invention proposes a multi-target recognition method for ground penetrating radar B-scan images. While identifying underground hyperbolic targets, underground non-hyperbolic targets are marked to achieve the effect of simultaneously identifying hyperbolic and underground cavity diseases.

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

[0007] A method for underground multi-target recognition based on slope distribution estimation and stationarity test includes the following steps:

[0008] Step 1, preprocess the image (mean cancellation, normalization, adaptive threshold segmentation) to convert the ground penetrating radar B-scan image into a binary image, and divide the binary image into a series of sub-images by row;

[0009] Step 2, using non-zero interval sparse coding (NISE) to perform non-zero interval sparse coding on the preprocessed image to convert the binary image into a sparse image;

[0010] Step 3, using column offset encoding (COE) to perform column offset encoding operation on the sparse matrix and filter out monomials that do not contain useful information;

[0011] Step 4, filter out useless information in the polynomial; cluster parts with similar slopes according to the distribution of offset and non-zero point position information;

[0012] Step 5: Design the criteria for distinguishing the targets represented by each category in each part.

[0013] Furthermore, the pre-processing process in step 1 includes the following steps:

[0014] Step 1.1, perform mean cancellation on the image to remove the direct wave;

[0015] Step 1.2, normalize the image after mean cancellation;

[0016] Step 1.3, perform adaptive threshold processing on the normalized image.

[0017] Furthermore, the mean cancellation method in step 1 is: for each pixel, the mean value of the row is subtracted from the value of the pixel. Let the original image be I o , the size is M×N, and its expression is as follows:

[0018]

[0019] Output I m Represents the image after mean cancellation, with a size of M×N.

[0020] Furthermore, the normalization method in step 1 is:

[0021] Normalization is to convert the values ​​of all pixels in the image to between [0,1]. The process can be expressed as:

[0022]

[0023] Output I n is an M×N matrix with values ​​ranging from [0,1], max(I m ) and min(I m ) represent the maximum and minimum values ​​of the whole image respectively.

[0024] Furthermore, the adaptive threshold algorithm in step 1 is:

[0025] The adaptive threshold algorithm consists of three parts:

[0026] 1. Use median filtering to remove image I n Noise, get image I b1

[0027] 2. For each pixel, if the value of the point is less than the average value of the global image, then the value of the point is set to be equal to m, and Mean{} is used to represent the matrix average value, which can be expressed mathematically as:

[0028]

[0029] 3. Using the adaptive threshold algorithm based on longitudinal gradient information, the gradient region set is defined as:

[0030]

[0031] Where g is the minimum gradient threshold, and the adaptive threshold T is defined as:

[0032] T=Mean(G(i,j))

[0033] After segmentation with threshold T, a binary image I is obtained. b :

[0034]

[0035] Furthermore, the method for equally dividing the image in step 2 is:

[0036] For image I b (i, j) is divided into m equal parts by row; if the number of rows of the image cannot be divided by m, and the quotient is a and the remainder is b, then the first b sub-images are a+1 rows, and the remaining sub-images are a rows; the kth sub-image is represented by C k .

[0037] Furthermore, the non-zero interval sparse coding (NISE) method in step 2 is:

[0038] For each sub-image, divide it into a set of column vectors. Suppose the k-th sub-image C k The size is R×N, then:

[0039] C k ={c 1 ,c 2 ,…,c j ,…,c N}

[0040] c j (1≤j≤N) is a column vector of length R, for c j Each element of is searched from top to bottom. If there is a non-zero element, the segment is searched for the number of consecutive non-zero elements. The end position of the segment is the position of the non-zero value in the sparse matrix, and the length of the segment is the size of the corresponding non-zero value in the sparse matrix.

[0041] Furthermore, the column offset encoding method in step 3 is:

[0042] S k Partition into a collection of column vectors:

[0043] S k ={s 1 ,s 2 ,…,s j ,…,s N}

[0044] Let the column vector s j In the example, the non-zero position set is Y j ={y 1,j ,y 2,j ,…,y i,j ,…|Value(y i,j ,j)≠0},Value(y i,j ,j) represents the yth i,j The output of the column offset encoding is a 1×N symbol vector, which is encoded as follows: k Each column in is cyclically shifted, using cir(s j ,r) means s j After the vector is circularly shifted by r positions, if at this time, cir(s j ,r)·s j+1 ≠0 (i.e. Then output a polynomial, let PubY j,r Yes j After circular shift by r positions, add j+1The set of non-zero points with overlapping positions is expressed as:

[0045]

[0046] The output polynomial is of the form:

[0047]

[0048] In this encoding method, each non-zero point that coincides with the next column in the cyclic shift is represented by a monomial; the frequency of each monomial represents the number of bits r of the cyclic shift at this time; the amplitude represents the sum of the i-th non-zero values ​​that coincide with the position of the next column after the previous column is cyclically shifted by r positions; the phase represents the position of the non-zero point after the cyclic shift.

[0049] Furthermore, the method for filtering out the unconnected items in step 3 is:

[0050] Assume that in the circular shift r p When s j With s j+1 The set of non-zero position overlaps is: Pick like

[0051]

[0052] This means that the two points in the original image are not connected, so the monomial is filtered out.

[0053] Further, step 4 includes the following steps:

[0054] 4.1. Calculate the Fourier series for the j-th column (1≤j≤N) COE code, extract all frequency terms of the polynomial, and let the set of all frequencies in the j-th column be fre j = {r 1,j ,r 2,j ,…,r i,j ,…}, traverse fre and find the offset that has not been clustered If fre j If all elements in are clustered, they go to the next column;

[0055] 4.2. Extract all the Similar offset, autonomously set bandwidth T, all offset components extracted are Construct a class diagram E with a size of R×N and an initial value of 0; add these elements to the class diagram E, and for all make

[0056] 4.3. For freT jFor each offset in, extract all the corresponding original image positions to form a set With these cen j All elements in are taken as the center point, and T is the bandwidth to construct a bandpass filter F;

[0057] 4.4. Extract all frequencies in the next column to form a set fre j+1 , and the corresponding original image position cen j+1 ; For each offset r in nfre i,j+1 , so that the bandpass filter F range is shifted cyclically and the offset is equal to r i,j+1 ; will cen j+1 Middle and r i,j+1 All the corresponding original positions are input into the bandpass filter to observe whether the output is empty; assuming the offset is r i1,j+1 If the output value is not empty, Put it in this class and make

[0058] Furthermore, in step 5, the centroid position S of each column is calculated according to the position of the non-zero point and the size of the non-zero point, and its differential sequence is represented by d(S); three judgment criteria are selected to judge the type of target: the trend of S, the stability of d(S), and the variance of S.

[0059] Beneficial effects:

[0060] 1. This invention is applied to the field of underground target recognition by ground penetrating radar. It is a non-training method and solves the problem that data in the field of ground penetrating radar is difficult to obtain.

[0061] 2. The present invention has high recognition accuracy and robustness.

[0062] 3. The present invention can simultaneously identify underground pipelines and hollow targets, and has a multi-target recognition effect that existing classical non-learning methods do not have. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 , pretreatment process and effect;

[0064] Figure 2 , non-zero interval sparse coding diagram;

[0065] Figure 3 , column offset coding and filtering of unconnected items;

[0066] Figure 4 , indication of the criteria for judgment;

[0067] Figure 5 , Experimental results of underground pipeline identification;

[0068] Figure 6, underground voids, voids and other disease experimental results. DETAILED DESCRIPTION

[0069] The specific technical solutions of the present invention are as follows:

[0070] Step 1: Preprocess the image (mean cancellation, normalization, adaptive threshold segmentation) to convert the GPR B-scan image into a binary image, and divide the binary image into a series of sub-images by row.

[0071] The preprocessing process includes the following steps:

[0072] Step 1.1, perform mean cancellation on the image to remove the direct wave;

[0073] Step 1.2, normalize the image after mean cancellation;

[0074] Step 1.3, perform adaptive threshold processing on the normalized image.

[0075] The specific method of mean cancellation is: for each pixel, subtract the average value of the row from the value of the pixel, assuming that the original image is I o , the size is M×N, and its expression is as follows:

[0076]

[0077] Output I m Represents the image after mean cancellation, with a size of M×N.

[0078] Normalization:

[0079] Normalization is to convert the values ​​of all pixels in the image to between [0,1]. The process can be expressed as:

[0080]

[0081] Output I n is an M×N matrix with values ​​ranging from [0,1], max(I m ) and min(I m ) represent the maximum and minimum values ​​of the whole image respectively.

[0082] Adaptive Threshold Algorithm:

[0083] The adaptive threshold algorithm consists of three parts:

[0084] 1. Use median filtering to remove image I n Noise, get image I b1

[0085] 2 For each pixel, if the value of the point is less than the average value of the global image, then let the value of the point be equal to m, and use Mean{} to represent the matrix average value, which can be expressed mathematically as:

[0086]

[0087] 3. Using the adaptive threshold algorithm based on longitudinal gradient information, the gradient region set is defined as:

[0088]

[0089] Where g is the minimum gradient threshold, and the adaptive threshold T is defined as:

[0090] T=Mean(G(i,j))

[0091] After segmentation with threshold T, a binary image I is obtained. b :

[0092]

[0093] The effect after pretreatment is as follows Figure 1 shown.

[0094] Step 2: Use non-zero interval sparse coding (NISE) to perform non-zero interval sparse coding on the preprocessed image to convert the binary image into a sparse image.

[0095] Image division:

[0096] For image I b (i,j) is divided into m equal parts by row. If the number of rows in the image cannot be divided by m, and the quotient is a and the remainder is b, then the first b sub-images are a+1 rows, and the remaining sub-images are a rows. The kth sub-image is represented by C k .

[0097] Non-zero interval sparse coding:

[0098] For each sub-image, divide it into a set of column vectors. Suppose the k-th sub-image C k The size is R×N, then:

[0099] C k ={c 1 ,c 2 ,…,c j ,…,c N}

[0100] c j (1≤j≤N) is a column vector of length R, for c jEach element of is searched from top to bottom. If there is a non-zero element, the segment is searched for the number of consecutive non-zero elements. The end position of the segment is the position of the non-zero value in the sparse matrix, and the length of the segment is the size of the corresponding non-zero value in the sparse matrix. The encoding method is as follows: Figure 2 shown.

[0101] Use S k Represents the sparse matrix of the K-th sub-image after NISE.

[0102] Step 3: Use column offset encoding (COE) to perform column offset encoding on the sparse matrix and filter out monomials that do not contain useful information.

[0103] Column offset encoding:

[0104] S k Partition into a collection of column vectors:

[0105] S k ={s 1 ,s 2 ,…,s j ,…,s N}

[0106] Let the column vector s j In the example, the non-zero position set is Y j ={y 1,j ,y 2,j ,…,y i,j ,…|Value(y i,j ,j)≠0},Value(y i,j ,j) represents the yth i,j row, value of column j.

[0107] The output of column offset coding is a 1×N symbol vector, which is encoded as follows: k Each column in is cyclically shifted, using cir(s j ,r) means s j After the vector is circularly shifted by r positions, if at this time, cir(s j ,r)·s j+1 ≠0 (i.e. Then output a polynomial, let PubY j,r Yes j After circular shift by r positions, add j+1 The set of non-zero points with overlapping positions is expressed as:

[0108]

[0109] The output polynomial is of the form:

[0110]

[0111] In this encoding method, each non-zero point that coincides with the next column in the cyclic shift is represented by a monomial. The frequency of each monomial represents the number of bits r of the cyclic shift at this time; the amplitude represents the sum of the i-th non-zero values ​​that coincide with the next column position after the previous column is cyclically shifted by r bits; the phase represents the position of the non-zero point after the cyclic shift; the column offset encoding process is as follows Figure 3 shown.

[0112] Filter out unconnected items:

[0113] In many cases, the point segments represented by the non-zero values ​​of the previous column are not connected to the next column in the original image, but in the cyclic shift, these unconnected point segments will overlap with the next column, and are thus represented by a monomial and added to the COE code; this information is not needed, so we need to filter out these items.

[0114] Assume that in the circular shift r p When s j With s j+1 The set of non-zero position overlaps is: Pick like

[0115]

[0116] This means that the two points in the original image are not connected, so the monomial is filtered out. Figure 3 shown.

[0117] Step 4, filter out useless information in the polynomial; cluster parts with similar slopes according to the distribution of offset and non-zero point position information;

[0118] The result of column offset encoding can reflect the slope information of the sub-image, and records the slope in the sparse matrix S k In the clustering process, if a non-zero position y in the previous column i,j The sum of its offset r is equal to a non-zero point position y in the next column k,j+1 , then the two non-zero points are in one category, and the corresponding offset r will also be recorded in this category.

[0119] The specific steps are as follows

[0120] Step 4.1. Calculate the Fourier series for the j-th column (1≤j≤N) COE code, extract all frequency terms of the polynomial, and let the set of all frequencies in the j-th column be fre j = {r 1,j ,r 2,j ,…,r i,j,…}, traverse fre and find the offset that has not been clustered If fre j If all elements in are clustered, they go to the next column.

[0121] Step 4.2. Extract all the Similar offset, autonomously set bandwidth T, all offset components extracted are Construct a class diagram E with a size of R×N and an initial value of 0. Add these elements to the class diagram E, and for all make

[0122] Step 4.3. freT j For each offset in, extract all the corresponding original image positions to form a set With these cen j All elements in are taken as the center points and T is the bandwidth to construct a bandpass filter F.

[0123] Step 4.4. Extract all frequencies in the next column to form a set fre j+1 , and the corresponding original image position cen j+1 For each offset r in nfre i,j+1 , so that the bandpass filter F range is shifted cyclically and the offset is equal to r i,j+1 . j+1 Middle and r i,j+1 All the corresponding original positions are input into the bandpass filter to observe whether the output is empty. Assume that the offset is If the output value is not empty, Put it in this class and make

[0124] Example:

[0125] Assume that the size of a sub-image is 100×400, the bandwidth is set to 5, and the bandwidth is . Assume that the COE polynomial of the jth column is

[0126] COE(j)=2e (x+11)i +4e (x+33)i +5e (3x+5)i +8e (6x+75)i +8e (97x+52)i

[0127] The COE polynomial of the j+1th column is

[0128] COE(j+1)=2e (2x+12)i +5e (2x+34)i +4e (5x+8)i +8e (5x+81)i +7e (95x+65)i

[0129] +8e (99x+49)i

[0130] The set of all frequencies in the jth column is

[0131] fre j ={1,3,6,97}

[0132] Traversing fre j , we find that offset 1 is not clustered, and extract all offsets whose distance from 1 is within the set bandwidth of 5:

[0133] freT j ={1,3,97}

[0134] Let the jth column, 1st, 3rd, 97th row in E be equal to 6, 5, 8 respectively. Extract the original position set corresponding to the offset 1, 3, 97 (that is, the set of all points in the original sparse image whose offset from the next row is 1 or 3)

[0135] cen j ={11,33,5,52}

[0136] cen j The elements in construct a bandpass filter for the center frequency:

[0137] F=[0,10][6,16][28,38][47,57]

[0138] The set of all frequencies in the j+1th column is

[0139] fre j+1 ={2,5,95,99}

[0140] For an offset of 2, the corresponding filter is

[0141] F r＝2 =[2,12][8,18][30,40][49,59]

[0142] The phase (original image position) corresponding to offset 2 is {12, 34}, which can pass the filter, so offset 2 belongs to this category.

[0143] For an offset of 5, the corresponding filter is

[0144] F r＝2 =[5,15][11,21][33,43][52,62]

[0145] The phase (original image position) corresponding to offset 5 is {8,81}, which is not empty after passing through the filter, so offset 5 belongs to this category.

[0146] For an offset of 95, the corresponding filter is

[0147] F r＝95 =[-5,5][1,11][23,33][42,52]

[0148] The phase (original image position) corresponding to the offset 95 is {65}, which is empty after passing through the filter, so the offset 95 does not belong to this category.

[0149] For an offset of 99, the corresponding filter is

[0150] F r＝99 =[-1,9][5,15][27,37][46,56]

[0151] The phase (original image position) corresponding to the offset 99 is {49}, which can pass the filter, so the offset 99 belongs to this category.

[0152] Therefore freT j ={1,3,97} and the offset {2,5,99} of the next column are in the same category. Let the j+1th column and the {2,5,99}th row in E be equal to 7, 12, and 8 respectively. Continue the loop until the next column has no offset belonging to this category or j+1=N.

[0153] For fre j Repeat the above steps for other unclustered elements in . In this clustering method, the same offset can be allowed to be clustered into different clusters.

[0154] Step five: Design the criteria for judging, and for each category of each part, determine the target it represents.

[0155] In step 4, each class diagram E represents a target to be identified, and step 5 determines the target based on the class diagram.

[0156] Calculate the centroid position S of each column according to the position and size of the non-zero point, and use d(S) to represent its difference sequence. Three criteria are used to determine the type of target: the trend of S, the stability of d(S), and the variance of S.

[0157] The stability parameter of d(S) is stable:

[0158] The echo of the target object (pipeline, etc.) in the measured data does not satisfy a perfect hyperbolic equation. Obviously, the distribution of the random jitter generated by each Ascan in the ground penetrating radar has nothing to do with the starting time. Assuming that the jitter caused by the uneven ground and the operator's posture satisfies the Gaussian distribution g t , the hyperbolic equation that satisfies the original scene is H(t), and the target echo equation in the actual b-scan image can be approximately expressed as:

[0159] b(t)=H(t)+g t

[0160] After column offset encoding:

[0161] b c (t) = β(t) + Δg t

[0162] The jitter of the hyperbola in the actual B-scan satisfies the trend stable distribution after column offset encoding. β(t) is the trend term, which is a linear equation of t. Δg t g t The difference.

[0163] For b 1 (t) Take the difference and get the difference sequence b d (t), then:

[0164] b d (t) = β 0 +ΔΔg t

[0165] β 0 is the first-order coefficient of β(t), the sum and difference of Gaussian distribution is still Gaussian distribution, ΔΔg t is Δg t The difference of also satisfies the Gaussian distribution. Therefore, the sequence b d (t) is a stationary process.

[0166] If the echo is caused by a hole, then β(t) is a piecewise function, and finally the difference is obtained. d (t) is a non-stationary series.

[0167] Use the ADF test to test the stationarity of the sequence d(S). If the test result is stationary, set stable = 1, otherwise stable = 0.

[0168] S's changing trend:

[0169] The process of column offset coding is equivalent to the process of deriving the image. Since the echo of the pipeline in the image approximately satisfies a downward-opening hyperbola, the center of gravity S after coding approximately satisfies a straight line with a downward trend. Extract the approximate slope of S. When the slope is downward, set trend = 1, otherwise trend = 0. The hollow echo is approximately the superposition of several small hyperbolas, and the overall shape is still downward-opening. Therefore, after coding, its shape is a non-straight line with a downward trend.

[0170] The variance sigma of S is:

[0171] When the jitter is too severe, it is likely that it is a false target composed of various clutters. This situation needs to be eliminated. The variance of S represents the jitter situation. When the variance of S is too large, it means that this category is clutter. When D(S) <T sigma When , let sigma = 1, otherwise sigma = 0, where T sigma Self-set.

[0172] Pipeline inspection method:

[0173] If stable = 1, trend = 1, sigma = 1, the echo is a hyperbola and is considered to be a pipeline target.

[0174] Void inspection method:

[0175] If stable=0, trend=1, sigma=1, the echo is other shapes with a downward opening, and this type is considered to be a hollow target.

[0176] The identification part is shown as follows Figure 4 shown.

[0177] In order to verify the accuracy of the invention in identifying underground targets, this experiment uses simulation data and measured data to test. The simulation experiment only judges the recognition effect of pipelines, and selects 400M, 900M, 1.5G, and 2.6GHz antennas. The background medium is soil or concrete, and is divided into four groups. The simulation experiment parameters of the four groups are shown in Table 1

[0178] Table 1 Simulation experiment parameters

[0179]

[0180] The recognition results are shown in Table 2

[0181] Table 2 Simulation data recognition performance

[0182]

[0183] The measured data uses 8 pipeline images and 3 hole images as input to verify the effectiveness of the method. The measured data parameters are shown in Table 3. Figure 5 For pipeline identification effect, Figure 6 This is the hole recognition effect.

[0184] Table 3 Measured experimental parameters

[0185]

[0186]

[0187] In summary, the above are only preferred embodiments of the present invention and are not intended to limit the protection scope of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for underground multi-target recognition based on slope distribution estimation and stationarity test, characterized in that: The steps include: Step 1, preprocess the image (mean cancellation, normalization, adaptive threshold segmentation) to convert the ground penetrating radar B-scan image into a binary image, and divide the binary image into a series of sub-images by row; Step 2, using non-zero interval sparse coding (NISE) to perform non-zero interval sparse coding on the preprocessed image to convert the binary image into a sparse image; Step 3, using column offset encoding (COE) to perform column offset encoding operation on the sparse matrix and filter out monomials that do not contain useful information; Step 4, filter out useless information in the polynomial; cluster parts with similar slopes according to the distribution of offset and non-zero point position information; Step 5: Design the criteria for distinguishing the targets represented by each category in each part.

2. The method according to claim 1, characterized in that The pre-processing process described in step 1 includes the following steps: Step 1.1, perform mean cancellation on the image to remove the direct wave; Step 1.2, normalize the image after mean cancellation; Step 1.3, perform adaptive threshold processing on the normalized image.

3. The method according to claim 1, characterized in that The method of mean cancellation in step 1 is: for each pixel, subtract the average value of the row from the value of the pixel. Assuming the original image is Io, the size is M×N, and the expression is as follows: Output I m Represents the image after mean cancellation, with a size of M×N.

4. The method according to claim 1, characterized in that The normalization method described in step 1 is: Normalization is to convert the values ​​of all pixels in the image to between [0,1]. The process can be expressed as: Output I n is an M×N matrix with values ​​ranging from [0,1], max(I m ) and min(I m ) represent the maximum and minimum values ​​of the whole image respectively.

5. The method according to claim 1, characterized in that The adaptive threshold algorithm described in step 1 is: The adaptive threshold algorithm part consists of three parts: 1). Use median filtering to remove image I n Noise, get image I b1 2). For each pixel, if the value of the point is less than the average value of the global image, then the value of the point is set to be equal to m, and Mean{} is used to represent the matrix average value, which can be expressed mathematically as: 3). Using the adaptive threshold algorithm based on longitudinal gradient information, the gradient region set is defined as: Where g is the minimum gradient threshold, and the adaptive threshold T is defined as: T = Mean(G(i,j)) After segmentation with threshold T, a binary image Ib is obtained: 。 6. The method according to claim 1, characterized in that The method for equally dividing the image in step 2 is: For image I b (i, j) is divided into m equal parts by row; if the number of rows of the image cannot be divided by m, and the quotient is a and the remainder is b, then the first b sub-images are a+1 rows, and the remaining sub-images are a rows; the kth sub-image is represented by C k .

7. The method according to claim 1, characterized in that The non-zero interval sparse coding (NISE) method described in step 2 is: For each sub-image, divide it into a set of column vectors. Suppose the k-th sub-image C k The size is R×N, then: C k ={c1,c2,…,c j ,…,c N } cj (1≤j≤N) is a column vector of length R. Each element of cj is searched from top to bottom. If there is a non-zero element, the number of consecutive non-zero elements that make up the segment is searched. The end position of the segment is the position of the non-zero value in the sparse matrix, and the length of the segment is the size of the corresponding non-zero value in the sparse matrix.

8. The method according to claim 1, characterized in that The column offset encoding method in step 3 is: divide Sk into a set of column vectors: S k ={s1,s2,…,s j ,…,s N } Suppose that in the column vector sj, the set of non-zero positions is Y j ={y 1,j ,y 2,j ,…,y i,j ,…|Value(y i,j ,j)≠0},Value(y i,j ,j) represents the yth i,j The output of column offset coding is a 1×N symbol vector, which is encoded by cyclically shifting each column in Sk and using cir(s j ,r) means s j The vector after circular shift r positions, if at this time, Then output a polynomial, let PubY j,r Yes j After circular shift by r positions, add j+1 The set of non-zero points with overlapping positions is expressed as: The output polynomial is of the form: In this encoding method, each non-zero point that coincides with the next column in the cyclic shift is represented by a monomial; the frequency of each monomial represents the number of bits r of the cyclic shift at this time; the amplitude represents the sum of the i-th non-zero values ​​that coincide with the position of the next column after the previous column is cyclically shifted by r positions; the phase represents the position of the non-zero point after the cyclic shift.

9. The method according to claim 1, characterized in that The method for filtering out unconnected items in step 3 is: Assume that when the cyclic shift is rp, the set of non-zero positions where sj and sj+1 coincide is Pick like This means that the two points in the original image are not connected, so the monomial is filtered out.

10. The method according to claim 1, characterized in that Step 4 includes the following steps: Step 4.

1. Calculate the Fourier series for the j-th column (1≤j≤N) COE code, extract all frequency terms of the polynomial, and let the set of all frequencies in the j-th column be fre j = {r 1,j ,r 2,j ,…,r i,j ,…}, traverse fre and find the offset that has not been clustered If fre j If all elements in are clustered, they go to the next column; Step 4.

2. Extract all the Similar offset, autonomously set bandwidth T, all offset components extracted are Construct a class diagram E with a size of R×N and an initial value of 0; add these elements to the class diagram E, and for all make Step 4.

3. freT j For each offset in, extract all the corresponding original image positions to form a set With these cen j All elements in are taken as the center point, and T is the bandwidth to construct a bandpass filter F; Step 4.

4. Extract all frequencies in the next column to form a set fre j+1 , and the corresponding original image position cen j+1 ; For each offset r in nfre i,j+1 , so that the bandpass filter F range is shifted cyclically and the offset is equal to r i,j+1 ; will cen j+1 Middle and r i,j+1 All the corresponding original positions are input into the bandpass filter to observe whether the output is empty; assuming the offset is If the output value is not empty, Put it in this class and make 11. The method according to claim 1, characterized in that In step 5, the centroid position S of each column is calculated according to the position and size of the non-zero point, and its differential sequence is represented by d(S); Three criteria are used to determine the type of target: the trend of S, the stationarity of d(S), and the variance of S.

Images