CN109738137B

CN109738137B - Real-time earth-rock dam leakage monitoring and rapid diagnosis method based on image contrast

Info

Publication number: CN109738137B
Application number: CN201910001572.4A
Authority: CN
Inventors: 王日升; 李居铜; 章传涛; 赵之仲
Original assignee: Shandong Jiaotong University
Current assignee: Shandong Jiaotong University
Priority date: 2019-01-02
Filing date: 2019-01-02
Publication date: 2021-05-07
Anticipated expiration: 2039-01-02
Also published as: CN109738137A

Abstract

The invention relates to an earth-rock dam leakage monitoring and diagnosing method based on image contrast, which is a nondestructive, dynamic, real-time monitoring and rapid monitoring method. Firstly, fixing an electrode on a dam body of an earth and rockfill dam, acquiring a three-dimensional resistivity chromatographic image in real time through an electrical method instrument, screening images based on convolutional neural network learning by adopting real-time transmission to acquire a mass of images, performing two-dimensional resistivity profile processing on a sample image, and judging through image color space conversion/color space separation/resistivity image comparison: detecting whether the mass center of the hidden danger abnormal area in the specimen image changes through an image mass center detection algorithm: or comparing the processed image with a damage threshold value through a Canny edge detection algorithm to judge that the measuring point is in a safe state. The method can monitor the leakage supervision of the earth-rock dam through image comparison, is consistent with the conclusion that no obvious leakage abnormity is found by a conventional detection means, and is rapid and intelligent.

Description

Real-time earth-rock dam leakage monitoring and rapid diagnosis method based on image contrast

Technical Field

The invention relates to the technical field related to earth and rockfill dam leakage monitoring and diagnosis.

Background

An invention application with the application number of 201610068583.0 was filed in 2016 by Sichuan university, and the document of the invention discloses a high-precision distributed optical fiber sensing monitoring technical scheme and system. The optical demodulation instrument adopts a high-precision prepulse Brillouin optical time domain analyzer PPP-BOTDA or Brillouin-Rayleigh synthesis system, and the sensing optical cable adopts an optical-electrical composite cable with sufficient electrothermal energy to ensure a mutation signal of electrothermal temperature rise amplitude at an optical fiber-infiltration line intersection so as to be reliably identified; the signal-to-noise ratio of the system is up to more than 10-20, so that the efficient and reliable observation of the infiltration line and the dam foundation underground water level can be realized. It can be used for the large-range-space-time double-covering on-line remote measurement of the earth and rockfill dam infiltration line and the dam foundation seepage-proofing curtain underground water level. The multifunctional tester has the advantages of obvious 'one machine with multiple functions', and can measure dam centralized leakage, reservoir silt and the like. The specific implementation mode is provided, and the specific implementation mode comprises the arrangement method of the sensing optical cables for monitoring the saturation line and the underground water level, namely the wave pattern and the technical details of optical cable laying, so that the survival rate of the optical fibers can be practically guaranteed, and the interference to dam construction is reduced to the maximum extent.

The application number of the Nanjing Water conservancy science research institute of the department of Water conservancy is 2017100742819, the Nanjing Water conservancy science research institute of the department of Water conservancy is in 2017, the disclosed device comprises a resistivity measurement system, a permeation pressure measurement system, a water level measurement system, a leakage amount measurement system, rainfall observation equipment, a cable and a detection box, wherein the resistivity measurement system and the leakage amount measurement system are related to the invention, the resistivity measurement system consists of a plurality of measurement electrodes, and the measurement electrodes can be distributed at corresponding parts of an earth-rock dam, such as longitudinal sections, and are used for collecting related data. And processing the resistivity image is realized by the following method, and effective conclusion is expected to be obtained: in the third step, the apparent resistivity is obtained through the resistivity measuring system, and the state value v is obtained through analysis and processing of the apparent resistivity, for example, the measured value of the apparent resistivity is multiplied by the area of the corresponding Thiessen polygon and accumulated, and a set of formula is disclosed for calculating the water level data, the leakage data and the rainfall data to be monitored to obtain a monitoring model.

The above is related art that the applicant can understand.

Disclosure of Invention

In order to solve the defects of the prior art, the invention provides a real-time monitoring and rapid diagnosis method for earth and rockfill dam leakage based on image contrast, and provides a nondestructive, dynamic, real-time monitoring and rapid monitoring method.

The technical scheme adopted by the invention for solving the technical problems is as follows:

the earth and rockfill dam leakage monitoring and diagnosing method based on image contrast comprises the following steps in sequence:

fixing an electrode on a dam body of the earth-rock dam, acquiring a three-dimensional resistivity tomography image in real time by an electrical method instrument, taking a section image in an initial stable state as an initial reference object for image comparison,

in the process, the quality, the size and the edge performance of the three-dimensional resistivity chromatographic images obtained at different time periods are completely consistent by strictly controlling the pre-buried position, the image size and the pixel size of the electrode;

step two, adopting real-time transmission, and transmitting the image data to a PC (personal computer) end through wireless transmission;

thirdly, screening the mass images obtained in the first step based on convolutional neural network learning, eliminating images without hidden dangers, deleting images overlapped in hidden danger overlapping areas, and carrying out scene labeling on the images with hidden dangers to obtain sample images;

step four, carrying out two-dimensional resistivity section processing on the specimen image obtained in the step three, carrying out gray processing on the obtained two-dimensional resistivity section image to obtain boundary images of various color images,

fifthly, image color space conversion is carried out;

sixthly, separating image color space; then, resistivity image comparison is carried out: the resistivity image comparison is carried out in the following mode or two modes according to the sequence:

the first method is as follows: detecting whether the mass center of the hidden danger abnormal area in the specimen image changes through an image mass center detection algorithm:

if so, alarming and processing the data,

if not, judging that the measuring point is in a safe state;

and secondly, comparing the processed image with a damage threshold value through a Canny edge detection algorithm, and if:

if the A/sample image change exceeds a damage threshold value, alarming and processing;

b/the specimen image changes and does not exceed the damage threshold value, the measurement point is judged to be in a safe state,

further, in the first step, three-dimensional forward and backward imaging or three-dimensional forward imaging processing is carried out on the data acquired by the field electrical method instrument.

Further, in the first step, the three-dimensional resistivity tomography image is one of a forward image or an inverse image.

Further, in the second step, in the wireless transmission, data transmission is performed between the electrical method instrument and the tablet computer through bluetooth, and 4G transmission is used between the tablet computer and the PC terminal.

Further, in the third step, the Gray value Gray is calculated according to a sensory weight value-specific gravity method, and then the images are subjected to in-situ one-to-one correspondence to complete the image graying processing.

Further, the destruction threshold is a series of images in a two-dimensional image interval from the time of initial stabilization to the time of destruction.

Further, the image color space conversion means that the RBG image represented by the color pattern is converted into an HSV map represented by the color characteristics.

The invention has the beneficial effects that:

based on Canny edge detection theoretical process and algorithm, the method can roughly obtain the positions of disease bodies of various hidden danger types of the earth and rockfill dam after three-dimensional forward and backward modeling processing is carried out on data collected by an on-site electrical method instrument, the hidden danger types are split along a certain position to obtain a two-dimensional image of a certain section hidden danger type, the two-dimensional image during initial stability is used as a reference object for comparing the image with an initial point, the image of the same position obtained in each time period is compared and analyzed with a threshold value, and diagnosis alarm prompt is carried out when the image threshold value exceeds the maximum limit of the earth and rockfill damage threshold value determined by test and numerical simulation. The method can monitor the leakage supervision of the earth-rock dam through image comparison, is consistent with the conclusion that no obvious leakage abnormity is found by a conventional detection means, is rapid and intelligent, and can know the specific effect in the embodiment part.

Drawings

FIG. 1 is a flow chart of an image contrast process obtained by an electrical method.

Fig. 2 is a neural network computation graph.

Fig. 3 is a resistivity profile.

Fig. 4 is a resistivity gradation processing chart.

FIG. 5 is a resistivity spatial color separation plot.

Fig. 6 is a flow chart of Canny edge detection.

FIG. 7 is a noise reduction plot of the resistivity profile at a standard deviation of 0.8.

FIG. 8 is a schematic diagram of gradient strength and direction calculation.

FIG. 9 is a schematic diagram of gradient strength and direction calculation.

Fig. 10 is a diagram illustrating the strong and weak boundaries.

Fig. 11 is a schematic diagram of color boundaries after being processed by the Canny edge detection algorithm.

Fig. 12 is a schematic diagram of resistivity profile processed by hough line detection algorithm.

Fig. 13 is before image binarization processing.

Fig. 14 shows the image after the binarization processing.

Fig. 15 is a centroid calculation area diagram (front) after etching treatment.

Fig. 16 is a centroid calculation region diagram (post) after the etching process.

FIG. 17 is an image centroid display diagram.

FIG. 18 is a two-dimensional resistivity profile of a potential volume location.

Fig. 19 is a schematic diagram of the resistivity after the graying processing and the Canny edge detection algorithm processing.

Fig. 20 is a hough line detection processing diagram.

Fig. 21 is an image calculation region map.

Fig. 22 is a diagram of image color space conversion processing.

Fig. 23 is a diagram of image color space separation processing.

FIG. 24 is a schematic of the subtraction of the resistivity at the same location at different times.

FIG. 25 is a diagram of the separation of the centroid detection method image transformations.

FIG. 26 is a comparison graph of image centroid at different time instances (time point one) at the same position.

Fig. 27 is a comparison graph of image centroid at different time points (time point two) at the same position.

Fig. 28 is a comparison graph of image centroid at different time points in the same position (time point three).

Fig. 29 is a comparison graph of image centroid at different time points in the same position (time point four).

FIG. 30 is a sectional view of section 205.

FIG. 31 is a sectional view of section 220.

FIG. 32 is a cut-away view of section 228.

Fig. 33 is a 205 sectional image subtraction result diagram.

Fig. 34 is a 220 cross-sectional image subtraction result chart.

Fig. 35 is a 228 sectional image subtraction result diagram.

Detailed Description

The method for monitoring and quickly diagnosing the leakage of the earth-rock dam in real time can also be called as a diagnosis system and a diagnosis model.

1. Implementation concept of diagnosis system

Based on images formed by a three-dimensional resistivity tomography method (related discussion about the three-dimensional resistivity tomography method, reference can be made to the Master graduate paper of Wangshengl, one of the inventors of the present application) and a test method, different hidden danger type damage thresholds of different earth-rock ratio dam body materials under different water environment effects are determined, and corresponding earth-rock dam leakage diagnosis systems are established by different image comparison algorithms under different forms.

2 diagnostic system establishment basis

The invention aims at the seepage of the earth-rock dam in different time periods to carry out the real-time acquisition, identification and diagnosis research of the resistivity image. In the process, in the contrast algorithm solving process, the quality, size, edge and other factors of the image affect the result of the diagnosis directly, the imaging data are subjected to different gray scale contrasts and edge conditions due to different electrode positions and imaging software in the acquisition process, and the different conditions bring difficulties for subsequent image processing, edge extraction and contrast calculation, so that in order to obtain a contrast result under the same conditions, the following imaging measures are adopted for process control in the embodiment:

and 2.1, fixing the electrode distance position. Because the difference of electrode ground connection condition can bring the influence to follow-up data collection to directly influence image quality, consequently in order to reach the data collection under the unified condition, all electrodes all adopt fixed position pre-buried, be about to the electrode pre-buried in the preset position that earth and rockfill dam needs monitoring, mark pre-buried electrode, will gather the instrument from taking the electrode to be connected with pre-buried electrode again, the kneck adopts insulating processing, avoids the inconsistent problem of ground connection condition that the electrode disturbance brought.

2.2, the same position image size is the same. In order to achieve an accurate one-to-one correspondence effect in the image comparison operation process, the sizes of two images to be compared before and after are required to be completely consistent, so that the sizes of the images are completely controlled in a fixed coordinate frame by adopting a coordinate control method in the imaging process, and the sizes of the images to be compared before and after at the same position can be ensured to be completely consistent.

And 2.3, image pixels are uniform. Since the image contrast is finally classified as pixel contrast, in order to ensure the consistency of the imaging pixels, the same control condition is applied to the imaging in advance, namely the pixel change condition meets the critical value of the resistivity, so that the change rate of each image pixel is stabilized in a controllable interval.

The diagnostic system is not only required to be accurate, but also efficient and rapid, which is one of the important consideration indexes, so that the rapid purpose can be achieved only by accelerating the process from each link of data acquisition, data transmission and image processing

And 2.4, wireless real-time data transmission. The electric method instrument for collecting data is improved, the collecting panel transmits the data to a tablet personal computer convenient to carry through Bluetooth, a mobile cos card is installed in the tablet personal computer, and the collected data is transmitted to a processed PC terminal through the flow of the card.

And 2.5, carrying out automatic image comparison processing. In order to achieve the purpose of rapid diagnosis and timely early warning of earth-rock dam leakage, massive image screening is performed on the basis of convolutional neural network learning, and scene labeling is performed on images with hidden danger bodies; and setting a threshold value (also called a damage threshold value) of the image variation range obtained by experiments based on a computer aided system, and providing an alarm for the variation range exceeding the threshold value. The threshold, also called a critical value, refers to the lowest value or the highest value that an effect can generate, and in this embodiment, the threshold is a series of images in a two-dimensional image interval from the first stabilization to the destruction.

3 Overall implementation of diagnostic Process

The earth and rockfill dam leakage detection diagnosis process is realized according to a three-dimensional high-quality resistivity chromatographic image formed in the detection process, the image can be a forward image or an inversion image according to actual needs, the forward image or the inversion image is subjected to two-dimensional subdivision, and the obtained two-dimensional subdivision image can be used as a specimen image for image comparison.

The method comprises the following steps of firstly adopting convolutional neural network learning to automatically screen dangerous region images of hidden patients for scene marking for images generated by a large number of measuring lines at different moments, or eliminating images without hidden danger bodies through screening, or deleting images in hidden danger overlapping regions, so that the number of compared images can be greatly reduced through convolutional neural network operation; and secondly, carrying out gray processing on the image with the screened saphenous body to reduce the calculation amount in the contrast process. The image subjected to gray processing needs to be subjected to image space conversion firstly, an RBG image represented by a color mode is converted into an HSV image represented by color characteristics, the HSV image subjected to space conversion is subjected to space separation, a channel which is only helpful for image comparison is separated, other color channels are ignored, the image color space subjected to gray processing is simple, the operation workload of a computer-aided system for image processing can be reduced, and a foundation is laid for rapid image comparison processing; and finally, performing image contrast analysis on the processed image according to different hidden danger types by adopting a Canny edge detection algorithm and a centroid movement algorithm, wherein the Canny edge detection algorithm can diagnose that the image threshold exceeds the damage threshold of the earth-rock-earth dam obtained by means of an object model and a digital model due to the expansion of the range of the hidden danger body, the centroid movement can warn that the position of the hidden danger body shifts although the range of the hidden danger body is not expanded, namely the earth-rock dam side slope caused by the hidden danger body slips, and the development condition of the hidden danger body is judged by the two different algorithms according to the result of the image contrast analysis.

The flow chart of the implementation described herein is shown in fig. 1/2.

Hereinafter, it is useful to understand the process and meaning of the present invention.

Regarding the massive image screening process based on convolutional neural network learning in the third step, in the scene labeling strategy based on convolutional neural network, firstly, the acquired image X, that is, the three-dimensional resistivity tomographic image obtained in the first step, is input, and may be regarded as being composed of multi-scale superpixels, that is, X ═ X₁,x₂,…x_n]And the input image X is a sample of convolution training. Convolutional neural network algorithm can provide a non-linear hypothesis model H_W,b(x)It has trainable parameters W, b, where W is the weight of the convolution kernel and b represents the convolution layer offset, whereby the convolution layer output in the neural network is:

wherein f is a nonlinear activation function, and a common activation function is a sigmoid function: (x) ═ 1+ exp^-x)^-1And hyperbalic tangent function: (x) ═ e^x-e^-x)/(e^x+e^-x)^-1。

To illustrate the neural network computation process in detail, take the example shown in FIG. 2, where x₁,x₂,x₃Is an entry, "+ 1" is an offset, and the rightmost side h_w,b(x) For the output item, the neural network calculation process is as follows:

a₁ ⁽²⁾＝f(W₁₁ ⁽¹⁾(x₁)+W₁₂ ⁽¹⁾(x₂)+W₁₃ ⁽¹⁾(x₃)+b₁ ⁽¹⁾)

a₂ ⁽²⁾＝f(W₂₁ ⁽¹⁾(x₁)+W₂₂ ⁽¹⁾(x₂)+W₂₃ ⁽¹⁾(x₃)+b₂ ⁽¹⁾)

a₃ ⁽²⁾＝f(W₃₁ ⁽¹⁾(x₁)+W₃₂ ⁽¹⁾(x₂)+W₃₃ ⁽¹⁾(x₃)+b₃ ⁽¹⁾)

h_W,b(x)＝a₁ ⁽³⁾＝f(W₁₁ ⁽²⁾(a₁ ⁽²⁾)+W₁₂ ⁽²⁾(a₂ ⁽²⁾)+W₁₃ ⁽²⁾(a₃ ⁽²⁾)+b₁ ⁽²⁾)

wherein: w_ij ^lAs associated parameters of the l +1 th layer neighbor cell, b_i ^(l)Is the bias term for the i cell in layer l +1, a_i ^lIs the activation value of the ith unit of the ith layer.

If z is used_i ^(l)Representing the ith cell input weighted sum of the ith layer and taking the activation function as a component representation, the above equation can be simply expressed as: z is a radical of⁽²⁾＝W⁽¹⁾x+b⁽¹⁾

a⁽²⁾＝f(z⁽²⁾)

z⁽³⁾＝W⁽²⁾a⁽²⁾+b⁽²⁾

h_W,b(x)＝a⁽³⁾＝f(z⁽³⁾)

Once the activation value a of the l-th layer is determined^lThen the activation value a of the l +1 th layer^l+1It can be obtained by the following formula:

z^(l+1)＝W^(l)a^(l)+b^(l)

a^(l+1)＝f(z^(l+1)) 5.3.2

in the calculation process of the convolutional neural network, if the t-th layer is the pooling layer, it can be expressed as:

H_t＝pool(H_t-1)

at this time, the output image is the input image H₀Class probability distribution Y of_(m)Can be expressed as^[145]：

Y_(m)＝P(L＝l_m/H₀；(W,b)

In the formula: m is an index tag, l_mAre category labels, both used as indexes.

In the convolutional neural network training process, the primary problem is to determine a training target, and the cost function is as follows:

to prevent overfitting, a weight decay term can be added to the cost function, so that the overall cost function of a fixed sample set containing m samples can be represented as:

the first term in the above formula is a mean square error term, the second term is a weight attenuation term, and lambda is a hyper-parameter, and the function of the hyper-parameter is to control the fitting degree of the whole network model.

For the above sample training, the objective of minimizing the cost function is achieved by changing the convolution kernel weight and the offset of the corresponding level by using a gradient descent method, that is:

w, b are updated with equation 5.3.3 for each iteration of the process, where α is the learning control, whose essential purpose is to control the decline of the velocity gradient during the learning process. When the loss rate of the trained neural network model is constant, the training is completed, and then the training parameters (W, b) can be determined accordingly, and the specific training steps are as follows:

1. l is calculated from the formula 5.3.2₂、L₃… … up to the output layer L_nlAn activation value of;

2. calculating the output residual error of each unit i of the output layer;

3. calculating a calculation residual error of any node i in any layer;

4. the final required partial derivative value is calculated.

The training steps aim to convert the spatial class probability of pixels of an image X into Y, label a scene L through the spatial distribution of different pixel points after conversion, and mark the pixel points by the common probability maximization principle in the labeling process^[145]Namely:

l_i＝max Y(i,j)(j∈[1,N]) 5.3.4

in the formula 5.3.4, i and j are respectively a pixel index and an identification index, and N is the total scene label number.

The convolutional neural network learns scene labeling, namely the class distribution probability of each pixel point of each image in the image space is obtained through training, so that an image labeling area can be expressed as a discrete random field:

G＝(V,E) 5.3.5

in the formula: v is the sum of all pixel points of each image, and E is the distribution of pixel points in different areas. The energy function of the discretized random field by equation 5.3.5 can be expressed as:

in formula 5.3.6, (i, i)^*) Respectively representing the combination of different pixel points; l_iRepresenting the identification category of the ith pixel point in the identification image; s represents the identification of different segmentation areas in the image; i is_cAll the pixel points in the segmentation region c are represented.

Thus, the whole function first order energy potential function can be defined:

in the formula, i and j have the same meanings as above, and sigma is a hyperparameter and can be freely set.

Defining a second order energy potential function as:

in the formula: x is the number of_iFor identifying pixel point corresponding values, x, in an image_i ^*Mu is a hyper-parameter and can be freely set for the RGB value of the corresponding point of the mapping area.

Defining a high-order energy potential function as:

in the formula: n (I)_c) Indicating the number of abnormal pixels, epsilon, in the partition_maxAnd expressing the quantity of all pixel points in the segmentation region, and Q expresses a truncation parameter.

Therefore, massive pictures can be screened through convolutional neural network learning, and the pictures screened out of the areas with hidden danger are subjected to scene marking.

In the fourth step, the processing is carried out by adopting a sensory weight value-specific gravity method in the image gray scale processing process, which has the advantages that the processing effect is more in line with the visual feeling of people, and the image processing meets the aim of simplicity and rapidness, therefore,

taking an abnormal body of a monitoring point of a certain earth and rockfill dam in the embodiment as an example, the data acquired by 32 electrodes of the abnormal body is three-dimensionally imaged and then is divided along a certain direction, so that a two-dimensional resistivity profile 3 (a depth resistivity change map) of the model in the certain direction is obtained.

And (3) calculating the Gray value Gray according to a sensory weight-to-weight method, then carrying out in-situ one-to-one correspondence on the image to finish the graying treatment of the image, wherein the treated image is shown in figure 4.

In the fifth step, in order to not change the color contrast effect of the color original image, in-situ assignment is performed on the image subjected to the graying processing so as to display the original color difference. Because the original color image generated by data is an RGB three-color channel image, the original image processing operation amount is large, and in order to convert the original RGB color image into an HSV color image represented by hue (H), saturation (S) and lightness (V) during simplification, the space conversion of the color image is realized, and the conversion formula for realizing the process is as follows:

V←max(R,G,B)

IfH is more than 0the H ← H +360.On output 0 is more than or equal to V and less than or equal to 1, S is more than or equal to 0 and less than or equal to 1, and H is more than or equal to 0 and less than or equal to 360. The brightness V is equal to the maximum value of three channels of R, G and B in each pixel, and the value range is between 0 and 1. The saturation S is given based on the lightness V, also in the range 0-1. And the hue H has different calculation formulas according to different values of the lightness V, and the value of the hue H is in the range of 0-360.

In the sixth step, the HSV color image mentioned in the fifth step is composed of three channels, and the color space is separated, namely pixel values of each channel are respectively extracted to form three gray level images. In the contrast process of the data generation image, only the color change amplitude is concerned, and the change of saturation and lightness is not required to be considered, so that the subsequent operation amount can be reduced by extracting only useful contrast information after the color space is separated, and the operation speed is greatly improved.

For the HSV color space, the hue H channel represents color information, and the variation range of the hue H channel is 0-360, but the variation range of the hue H channel in OpenCV is 0-180, and the hue H channel is linearly extended to 0-255, so that the value distribution interval is the same as the gray pixel interval. In the separation operation extraction process, three-channel transformation is firstly carried out according to the color image after space transformation, S, V channels are eliminated and shielded according to the transformed result, only the hue H is separated independently, and the color information is used as the quantization value of the three-dimensional resistivity profile so as to complete the space separation of the image color. The results after spatial color separation are shown in fig. 5.

In step seven, the Canny edge detection algorithm reduces the data size of the image operation as much as possible without changing the attribute of the original image, which is generally based on the gray image, and the implementation flow process thereof is shown in fig. 6.

The gaussian filtering denoising algorithm is a method for eliminating noise points generated by a computer in the process of generating an image due to difference in a calculation process in the process of outputting the image.

The main parameter in the process of generating the Gaussian filter template is the standard deviation sigma, and due to the data discrete degree, if sigma is smaller in processing, the central coefficient of the generated template is larger, but the coefficient around the center is smaller, so that the image smoothing processing effect is not obvious; conversely, if σ is large, the image smoothing processing effect is good. In general, σ may be 0.8, and the result after denoising is shown in fig. 7.

And after noise points are removed through Gaussian filtering and noise reduction, gradient calculation of the identification image can be carried out.

The general idea of the Canny algorithm is to obtain the position of the coordinate point with the maximum gray level intensity change in each recognition image. The location in the image where the change is greatest can be defined as the direction of the gradient of the change. The calculation process is as follows:

1. firstly, acquiring a gradient map of an identification image, namely respectively calculating the horizontal (X direction) derivative and the vertical (Y direction) derivative of the image after gray processing by adopting a Sobel operator, and taking the first derivative as the gradient map.

G_x(i,j)＝(f[i+1,j-1]-f[i,j-1]+f[i+1,j]-f[i,j]+f[i+1,j+1]-f[i,j+1])/2

G_y(i,j＝)(f[i-1,j]-f[i-1,j+1]+f[i,j]-f[i,j1]+f[i+1,j]-f[i+1,j+1])/2

2. From the obtained two gradient maps Gx and Gy, the gradient and direction of the image boundary are determined by equations 5.4.1 and 5.4.2.

After the gradient size and the gradient direction of the resistivity identification image are obtained, corresponding image contrast processing measures can be adopted to remove the interference points with larger influence. The removal method comprises the following steps: firstly, gradient calculation is carried out on pixels of each point in the resistivity image, whether the gradient of the pixels is the maximum gradient value in adjacent points in the same gradient direction or not is judged according to the calculation result, front and back comparison is carried out if the gradient direction is horizontal, and up and down comparison is carried out if the gradient direction is vertical.

As shown in fig. 8, the number represents the gradient strength of the pixel point, the arrow direction represents the gradient direction, and the gradient direction is vertical, so the comparison is performed in the vertical direction. Taking the third pixel point in the second row as an example, since the gradient direction is upward, the intensity 7 of the third pixel point is compared with the

intensities

5 and 4 of the upper and lower two pixel points, since the intensity of the position of the 7 is the largest in the intensity comparison of the upper and lower adjacent three points, the 7 is retained, the adjacent 5 and 4 are removed, and the like, the maximum intensity value along the gradient direction in each comparison graph can be obtained.

If the gradient direction is in either direction, as shown in FIG. 9, the intensity values M (C) in the horizontal and vertical directions are obtained₁)、M(C₂)、M(C₃)、M(C₄) It is known that non-maximum suppression can be performed by linear interpolation,the specific method comprises the following steps:

wherein: n-dis distance (C)₁₂,C₂)/dis tance(C₁,C₂)

The gradient value in any gradient direction can be interpolated by the formula 5.4.3, so that the calculated result can be compared in any gradient direction for non-maximum suppression.

Many noise points still exist in the image after the non-maximum suppression processing. Therefore, in order to obtain better image contrast result, the Canny algorithm removes noise by setting the upper and lower threshold (usually specified by human in opencv) techniques, and therefore also becomes dual-threshold control, that is, a control range of an image is defined by setting one threshold.

The double-threshold selection is considered that the image is identified to be a strong boundary when the pixel point value exceeds the upper boundary of the set threshold, the image is identified to be a weak boundary when the pixel point value is smaller than the lower boundary of the set threshold, and the weak boundary is located between the upper boundary and the lower boundary and belongs to the boundary processing category. As shown in FIG. 10, point A is greater than the set upper threshold bound and therefore a strong boundary, belonging to a true boundary, and point C is less than the upper threshold bound but greater than the lower threshold bound, which can be considered to belong to a boundary point because its boundary is connected to A. While point B, although lower than the upper threshold maxVal and higher than the lower threshold minVal, should be discarded because it is not connected to a true boundary point. The ideal boundary of each color contrast of the image can be obtained through double threshold selection.

The color boundaries corresponding to the noise reduction process of fig. 7 can be obtained by calculating through the processes of the Canny edge detection algorithm, as shown in fig. 11, a color boundary schematic diagram after the Canny edge detection algorithm is processed.

The hough line detection is to detect a line in an image by using hough transform, which is a feature extraction technology in image processing, and in this embodiment, the following procedure is followed:

1. firstly, acquiring edge information of each contrast recognition resistivity image based on a Canny edge detection algorithm;

2. drawing a straight line in k-b space for each point in the edge image;

3. accumulating the points on each straight line by adopting a 'voting' method, namely if one straight line passes through the point, adding 1 to the value of the point until the accumulation of the straight line passing through the point is finished;

4. and traversing the k-b space to find local maximum value points, wherein the coordinates (k, b) of the points are the slope and the intercept of a possible straight line in the original image.

The resistivity profile after the hough line detection algorithm is shown in fig. 12, where the detection line is marked in blue.

Therefore, the image contrast frame after Hough line detection is completely consistent, and Hough line interception can be carried out on the in-situ image.

In the step eight, an image centroid position shift algorithm is adopted in the implementation process of the other algorithm of the image contrast, and the method is suitable for the position shift of the saphenoid in the two contrast images.

The process applies the following technology, image binarization processing,

The image binarization processing is to make the image which participates in the contrast only show two different effects of non-black, namely white and black, namely, to turn the contrast image pixels to 0 or 255.

There are also different algorithms for the selection of the threshold. In this example, for anomalous resistivity cross-sections, because they have the same size and background, we chose the Otsu algorithm to manually select the threshold in order to achieve the desired effect. After many comparison and debugging, the result shows that the threshold value is selected to obtain the best effect at 233, and the image can be converted into a black and white effect picture with the pixel value of 255 when being larger than 233 and 0 when being smaller than 233 according to the result. The original image before processing and the effect map are shown in fig. 13/14.

In this embodiment, a morphology-erosion operation is used in the operation process of processing a binary image, after the image in fig. 14 is eroded, the coordinates and the basic noise of the binarized image may be removed, and then the area image involved in the centroid calculation is detected and calculated by using the contour detection function in opencv, and the area image extracted after the contour detection is as shown in fig. 15 and fig. 16.

The centroid in the image contrast refers to the average position of the pixel mass distribution in the particle system. If the contrast image is composed of n pixel points, the quality of the composition unit is m₁，m₂……m_nIf r is₁，r₂……r_nRepresenting the vector of each pixel point to any fixed point O, and being represented as r₀Namely:

and f (x, y) is more than or equal to 0 for a two-dimensional discrete image. Moment M of order p + q_pqAnd central moment mu_pqIs defined as:

formula 5.5.4 i_c，j_cIs a centroid coordinate and has:

i_c＝m₁₀/m₀₀ 5.5.5

j_c＝m₁₀/m₀₀ 5.5.6

as can be seen from equations 5.5.5 and 5.5.6, the centroid of the image after the erosion processing is the ratio of the 0th order moment to the 1 st order moment, so that the centroid of each image in the compared image can be obtained, and the centroid of the image can be labeled, as shown in fig. 17, the centroid of the image is the vertex of the lower right corner of the black rectangle.

The first embodiment is described by taking a certain actual project as an example:

1. first, a two-dimensional sectional view of a saphenous body at a certain position is taken, and the sectional view in the initial steady state is taken as an initial reference for image comparison, as shown in fig. 18.

2. The two-dimensional image obtained above is subjected to graying processing, and then image edge detection is performed by using a Canny algorithm, so as to obtain boundary images of color images, as shown in fig. 19.

3. And (4) carrying out Hough line detection on the original image, filtering out lines which are not concerned by people, obtaining a horizontal line of the section triangle, and drawing the line in blue. As shown in fig. 20.

4. After the hough line detection, two lines on the upper side and the left side can determine a rectangle, and the image is captured by taking the rectangular area as a calculation area as shown in fig. 21.

5. And performing Gaussian filtering processing on the determined calculation area image, namely converting the image from an RGB color space to an HSV color space. The resulting image is shown in fig. 22.

6. Further carrying out image space color separation on the processed image, separating out image color channels, and changing the color sequence: red, orange, yellow, green, cyan, blue, violet, corresponding to a change in value of [0, 255], isolated image display is shown in fig. 23.

7. The same processing operation is performed on the images at the same position obtained in any time interval in the steps of 1 to 6 to obtain the spatial color separation images.

8. Two space color separation images at the same position at different time before and after are subjected to phase diagram position color amplitude subtraction, particularly potential hazard position subtraction, so as to observe the resistivity change condition of the potential hazard position, and the schematic diagram after image subtraction is shown in fig. 24.

9. And averaging the amplitudes of the image pixels after the subtraction of the images at the color edge space position to obtain the threshold value of the resistivity change of the front image and the rear image.

10. And (3) comparing the average threshold value of the resistivity change obtained by the step (9) with the resistivity change threshold value obtained by destroying different soil-rock ratio models obtained by the test under different water-containing conditions, and if the resistivity change threshold value is larger than the threshold value obtained by the test, carrying out diagnosis and early warning, otherwise, judging as safe.

The specific implementation process of the preprocessing technology based on the resistivity two-dimensional profile and the algorithm for solving the image centroid is as follows:

1. and (5) image preprocessing. First, an arbitrarily obtained two-dimensional resistivity profile in a certain direction (as shown in fig. 13) is subjected to image binarization processing, the grayscale image is used as input, image binarization processing is performed, a binarization threshold value is obtained by an algorithm of Otsu, and the image after the binarization processing is shown in fig. 14.

2. And carrying out corrosion operation on the image after the binarization processing. The purpose of the etching operation is to remove noise points and coordinate systems with large interference in the image, and the effect after the etching treatment is as shown in fig. 15. The centroid calculation area of the image is extracted by using the opencv image contour detection function, and the result is shown in fig. 16.

3. Image color space conversion and separation. The preprocessed image is used for color space conversion, the RGB image is converted into an HSV color space image, the converted image is subjected to color separation, H-channel color information is extracted, and the converted and separated image is displayed as shown in fig. 25.

4. And (4) solving the centroid of the image. The centroid diagrams of the front and rear images at the same position are obtained by the equations 5.5.1-5.5.6, as shown in fig. 26/27/28/29, it is obvious that the centroid positions of the images at four different moments are shifted to the lower left.

In fact, the preprocessing technology based on the resistivity two-dimensional profile and the algorithm for finding the image centroid and the process and algorithm based on the Canny edge detection theory are parallel processes or alternative processes.

The advantages of the invention, a typical section resistivity imaging diagnostic (example), are illustrated by the following example.

And selecting three typical slice sections, amplifying the three typical slice sections, and then carrying out image monitoring contrast at the same position for use. In order to clearly show the resistivity change in the dam body, the monitoring time is selected to be 2017.12.2 and 2018.5.9, the section cutting surface is shown as 30/31/32, the three sections are respectively displayed in the elevation ranges of 205, 220 and 228, low-resistance abnormal areas exist, and the section after treatment shows that the radiation range of the low-resistance abnormal areas is not expanded, namely, the development trend of forming through channels does not exist.

Since no sliding phenomenon of the position of the abnormal body is found in all the slice images, the Canny edge detection early warning algorithm implementation process is adopted to process each comparison image respectively, and since the specific processing process is described in detail above, the detailed description is omitted here, and a subtraction result image after each group of images is processed is directly given, as shown in fig. 33/34/35.

It can be known from calculation of fig. 33/34/35 that, by subtracting the same position maps measured in the two previous and subsequent times to obtain values (103.077, 102.937, 104.111, respectively) with color average values of 14.736, 13.518, 11.773, which are obtained after color space transformation and separation from the three images measured in the first time, ratio processing is performed, so that the change rates of the three cross-sectional images monitored in the two previous and subsequent times are 14.3%, 13.1%, and 11.3%, respectively, and the three cross-sectional images do not reach the alarm condition of the diagnosis threshold when the earth-rock ratio is 7:3, so that the image comparison conclusion shows that the leakage is safe, which is consistent with the abnormal conclusion obtained by the conventional detection means that no obvious leakage is found.

It can be further known from the analysis of fig. 33/34/35 that there is an abnormal high-resistance region at the edge of the dam crest, on one hand, because the dam crest is far away from the water-facing region, the water content of the particles in the dam body is relatively low, and on the other hand, a local high-resistance region is formed after the cracks at the dam crest are filled with air, but it is clear from the comparison image that the apparent cracks are longitudinally developed but basically unchanged in the depth direction, and no through-connection occurs. The color change rate of the entire contrast image is also reduced if the average color value changes due to high resistance at the top surface of the dam is planed off.

In fact, whether the centroid of the hidden danger abnormal area in the specimen image changes is detected through an image centroid detection algorithm: and comparing the processed image with a damage threshold value through a Canny edge detection algorithm, wherein the two judgment processes can be used independently for judgment and can also be used in a combined way for judgment, and the judgment is within the protection scope of the invention.

The above-described embodiments are merely illustrative of the preferred embodiments of the present invention and do not limit the scope of the present invention, and various modifications and improvements of the present invention by those skilled in the art without departing from the spirit of the present invention are intended to fall within the scope of the present invention defined by the claims.

Claims

1. The real-time earth-rock dam leakage monitoring and rapid diagnosis method based on image contrast comprises the following steps in sequence:

fixing an electrode on a dam body of an earth-rock dam, acquiring a three-dimensional resistivity tomography image in real time through an electrical method instrument, taking a section image in an initial stable state as an initial reference object for image comparison, and performing three-dimensional forward and backward imaging or three-dimensional forward imaging processing on data acquired by the field electrical method instrument;

in the process, the pre-buried position, the image size and the pixel size of the electrode are controlled, so that the quality, the size and the edge performance of the three-dimensional resistivity tomography images acquired at different time periods are completely consistent;

step two, real-time transmission, namely transmitting the image data to a PC (personal computer) end through wireless transmission;

step three, screening the massive images obtained in the step one based on convolutional neural network learning, and screening the massive images based on convolutional neural network learning, wherein in a scene labeling strategy based on convolutional neural network, firstly, an acquired image X is input, namely, the three-dimensional resistivity tomography image obtained in the step one is regarded as being composed of multi-scale super pixels:

i.e. X ═ X₁,x₂,...x_n]，

In the formula: x is the number of₁Representing the first acquired image, x₂Representing the second acquired image, x_nRepresenting the image acquired the nth time,

the collected image X is a sample of convolution training, and a convolution neural network algorithm can provide a nonlinear hypothesis model H_W,b(x)It has trainable parameters W, b, where W is the weight of the convolution kernel and b represents the convolution layer offset, whereby the convolution layer output in the neural network is:

wherein f is a nonlinear activation function, Wi is the ith convolution kernel weight, i represents the acquisition times, x represents the acquired image,

the activation function is a sigmoid function: (x) ═ 1+ exp^-x)^-1

And hyperbalic tangent function: (x) ═ e^x-e^-x)/(e^x+e^-x)^-1；

In the formula, x represents the acquired image, exp represents a function with e as the base,

H_t＝pool(H_t-1)

in the formula: t is the layer of the pool formation,

at this time, the output image is the input image H₀Class probability distribution Y of_(m1)Can be expressed as:

Y_(m1)＝P(L＝l_m1/H₀；(W,b)

in the formula: m1 is an index tag, l_mThe label is a category label, the category label and the category label are used as indexes, P is a probability function, L is a loss function, W is the weight of a convolution kernel, and b represents the offset of the convolution layer;

to prevent overfitting, a weight decay term is added to the cost function, so that the overall cost function of a fixed sample set containing m samples can be represented as:

in the above formula, the first term is a mean square error term, the second term is a weight attenuation term, λ is a hyper-parameter, i and J are respectively category index labels, J is a cost function, x is an acquired image, y is a fixed sample number, h is a mean square error, m is a sample number, and sl is a sample attenuation number, and the function of the method is to control the fitting degree of the whole network model;

aiming at the sample training, a gradient descent method is adopted, and the goal of minimizing the cost function is achieved by changing the convolution kernel weight and the offset of the corresponding level, namely:

in the formula, l represents the number of layers, i and j are respectively pixel point class label indexes, and alpha is learning control

Updating W and b in each iteration in the process, wherein alpha is learning control and controls the reduction of the speed gradient in the learning process; when the loss rate of the trained neural network model is constant, the training is completed, and then the training parameters (W, b) can be determined accordingly, and the specific training steps are as follows:

1. calculating to obtain L₂、L₃… … up to the output layer L_nlAn activation value of;

2. calculating the output residual error of each unit i of the output layer;

3. calculating a calculation residual error of any node i in any layer;

4. calculating a final required partial derivative value;

in the formula, delta_iFor training the output residual, f is the training parameter, l is the number of training layers,

the purpose of the training steps is to perform space classification probability conversion Y on the pixels of the collected image X, perform scene L labeling on the converted pixels through the space distribution of different pixels, and identify the pixels by using a probability maximization principle in the labeling process, namely:

l_i＝max Y(i,j)(j∈[1,N])

in the formula, i and j are respectively a pixel index and an identification index, and N is a total scene labeling number;

G＝(V,E)

in the formula: v is the sum of all pixel points of each image, and E is the distribution of different areas of the pixel points;

the energy function of the discretized random field can be expressed as:

i，i^*respectively representing the combination of different pixel points; l_iRepresenting the identification category of the ith pixel point in the identification image; s represents the identification of different segmentation areas in the image; i is_cAll pixel points in the segmentation region c are represented; phi, phi,

Gamma is respectively a corresponding calculation function;

thus, the whole function first order energy potential function can be defined:

in the formula, i and j have the same meanings as above, and sigma is a first-order hyperparameter which can be freely set;

defining a second order energy potential function as:

in the formula: x is the number of_iFor identifying pixel point corresponding values, x, in an image_i ^*The RGB value of the corresponding point of the mapping area, mu is a second-order hyper-parameter which can be freely set, and exp is a function with e as a base;

defining a high-order energy potential function as:

in the formula: n (I)_c) Indicating the number of abnormal pixels, epsilon, in the partition_maxExpressing the number of all pixel points in the segmentation region, and Q expressing a truncation parameter;

screening a large number of pictures through convolutional neural network learning, eliminating images without hidden dangers, deleting images with overlapped hidden dangers, and carrying out scene labeling on the images with hidden dangers to obtain sample images; calculating a Gray value Gray according to a sensory weight value-to-weight method, and then carrying out in-situ one-to-one correspondence on the image to finish the graying processing of the image;

fifthly, image color space conversion is carried out;

the image color space conversion process is as follows: carrying out Hough line detection on the original image, filtering out unappreciated straight lines to obtain a horizontal straight line with a triangular cross section, and drawing the straight line in blue;

determining a rectangle by two straight lines on the upper side and the left side after Hough straight line detection, and intercepting an image by taking the rectangular area as a calculation area;

performing Gaussian filtering processing on the determined image of the calculation area to obtain an image;

sixthly, separating image color space; and then carrying out resistivity image comparison, wherein the resistivity image comparison is carried out in a mode of selecting the following mode or two modes and combining the modes according to the sequence:

detecting whether the centroid of the hidden danger abnormal area in the sample image changes through an image centroid detection algorithm, wherein the centroid in the image contrast refers to the average position of pixel mass distribution in a mass point system;

if so, alarming and processing the data,

if not, judging that the measuring point is in a safe state;

comparing the processed image with a damage threshold value through a Canny edge detection algorithm, wherein the damage threshold value is a series of images in a two-dimensional image interval from the initial stabilization to the destruction;

if:

A. if the change of the specimen image exceeds a damage threshold value, alarming and processing;

B. and if the change of the specimen image does not exceed the damage threshold value, judging that the measuring point is in a safe state.

2. The image-contrast-based real-time earth and rockfill dam leakage monitoring and rapid diagnosis method according to claim 1, wherein in the first step, the three-dimensional resistivity tomography image is one of a forward image or an inverted image.

3. The image contrast-based real-time earth and rockfill dam leakage monitoring and rapid diagnosis method according to claim 1, wherein in the second step, in the wireless transmission, data transmission is performed between an electrical method instrument and a tablet computer through Bluetooth, and 4G transmission is used between the tablet computer and a PC terminal.