Overwater pollution event situation analysis inversion method

Abstract

The invention discloses a water pollution event situation analysis inversion method, and relates to the field of remote sensing technology and environment model technology, the method comprises the following steps: carrying out feature fusion on an obtained infrared image and an ultraviolet fluorescence image, building an oil film thickness calculation model based on the fused comprehensive features, and obtaining the observation concentration of a detection area by using the model; establishing a convective diffusion model of the target water area, introducing a loss item related to time and flow velocity, estimating the predicted concentration of the detection area by using the convective diffusion model, establishing an inversion model based on the observed concentration and the predicted concentration, and speculating the leakage position and leakage time in the target water area through an optimization iteration mode; according to the technical scheme, the infrared image and the ultraviolet fluorescence image are subjected to feature fusion, the oil concentration of the detection area can be accurately calculated, the loss item related to the time and the flow velocity is introduced into the convection diffusion model, and the inversion precision of the leakage position and the leakage time is improved.

Description

Technical Field

[0001] The present application relates to the fields of remote sensing technology and environmental model technology, and in particular to a method for inversion of situation analysis of water pollution incidents. Background Art

[0002] Due to the limitations of transportation costs and carrying capacity, water transportation has become the main form of oil transportation. With the prosperity and development of the social economy, the demand for oil transportation is strong, and with it comes the frequent occurrence of oil spills and water pollution incidents. Oil water transportation uses tankers or oil barges as carriers, and the transportation range covers waters such as oceans, rivers and lakes. The vast majority of my country's imported crude oil is transported to the country by tankers, and major refineries are also built along the coast or around navigable rivers and lakes. Most of the refined oil is also completed by coastal and inland shipping. Therefore, the frequent occurrence of water pollution incidents has become the focus of social attention. By developing and applying advanced water pollution incident situation analysis inversion methods, we can effectively identify and track pollution sources, predict pollutant diffusion paths, provide strong support for environmental management and emergency response, and thus protect the water ecological environment and public health.

[0003] In the prior art, a single detection method or field measurement is usually used to detect the oil film thickness of the selected target waters, which is easily affected by external background light or other waves, such as infrared remote sensing, ultrasonic detection and other methods. The accuracy of the oil film thickness detected by a single method is low and cannot meet the requirements of subsequent inversion calculations. The film thickness detection process using field measurement is cumbersome and difficult to implement. The detection cycle of this detection method is long and data cannot be obtained in the same time period. In addition, in the prior art, during the subsequent inversion analysis process, the amount of oil lost due to external factors is usually ignored when establishing the model, resulting in the oil spill concentration ultimately calculated by the model being not accurate enough, or linear or nonlinear reaction terms are used when establishing the model to describe the lost part of the oil volume, and the coefficient of the reaction term has nothing to do with the flow rate of the river itself and the length of time, and cannot be well applied to different waters, and the calculation accuracy is relatively low. Summary of the invention

[0004] The purpose of this application is to solve the problem in the prior art that the oil film thickness detection is not accurate enough and the description of the amount of oil lost during the water flow is not accurate enough, which leads to the inability of subsequent inversion calculation to accurately infer the oil leakage point and leakage time.

[0005] The technical solution of the present application is to provide a method for inversion of situation analysis of water pollution incidents, the method comprising:

[0006] Step 1, obtaining an infrared image and an ultraviolet fluorescence image of the detection area, and performing image preprocessing respectively, and registering the preprocessed infrared image and ultraviolet fluorescence image so that the pixel positions of the two registered images are aligned;

[0007] Step 2, extracting features and contrast information of corresponding pixels in the infrared image and the ultraviolet fluorescence image, assigning weights according to the contrast information of corresponding pixels in the two images, and performing feature fusion on corresponding pixels in the two images based on the weights;

[0008] Step 3: Establish an oil film thickness calculation model based on the fused features, use the oil film thickness calculation model to calculate the oil film thickness corresponding to each pixel in the fused image, and calculate the oil spill concentration corresponding to each pixel according to the oil film thickness as the observed concentration;

[0009] Step 4, establish a convection diffusion model for the target waters, introduce loss terms related to time and flow rate into the convection diffusion model, preset the initial oil spill location and initial oil spill time, divide the target waters based on the spatial coordinate system, discretize the convection diffusion model, and calculate the predicted concentration of each detection point in the spatial coordinate system through the discretized model;

[0010] Step 5: Establish an inversion model based on the observed concentration and the predicted concentration calculated by the convection-diffusion model, iteratively optimize the parameters of the inversion model, and obtain the optimal values ​​of the initial oil spill position and initial time.

[0011] Furthermore, step 2 specifically includes:

[0012] Extract the radiation temperature T of each pixel in the infrared image IR (i, j), texture feature F IR (i, j) and thermal contrast C IR (i, j), extract the fluorescence intensity I of each pixel in the UV fluorescence image UV (i, j), texture feature F UV (i,j) and fluorescence contrast C UV (i,j), according to the extracted thermal contrast C IR (i,j) and fluorescence contrast C UV (i, j) assigns weights to the two corresponding pixels in the infrared image and the ultraviolet fluorescence image:

[0013]

[0014] In the formula, ω IR (i,j) is the weight of the pixel in the infrared image, ω UV (i, j) is the weight of the pixel in the UV fluorescence image, (i, j) is the coordinate label of the pixel, ω IR(i,j) and ω UV The sum of (i, j) is 1, and the registered infrared image and ultraviolet fluorescence image have the same reference coordinate system;

[0015] Based on the weight distribution of different pixel positions, the pixels in the two images are fused in turn to obtain a fused image. Each pixel in the fused image corresponds to a comprehensive feature vector P fusion (i,j):

[0016] P fusion (i,j)=[ω IR ·T IR (i,j),ω UV I UV (i,j),ω IR ·F IR (i,j),ω UV ·F UV (i,j)] T

[0017] In the formula, [·] T Represents the transpose of a vector.

[0018] Furthermore, in step 3, establishing an oil film thickness calculation model based on the fused features specifically includes:

[0019] Based on the comprehensive feature vector P fusion (i,j) Establish the oil film thickness calculation model:

[0020]

[0021] Where h is the oil film thickness, β0 is the intercept term, and β n is the regression coefficient, P fusion,n is the nth component of the comprehensive feature vector, ∈ is the error term, n is the label of the pixel in the fused image, and N is the total number of pixels in the fused image;

[0022] Optimize the oil film thickness calculation model and calculate the parameters β0 and β n The optimal value of , converts the oil film thickness calculation model into:

[0023] H=Xβ+∈

[0024] Where H is the actual observation value matrix of oil film thickness, X is the input feature matrix, and β is the intercept term β0 and the regression coefficient β n Assuming that a total of M samples are measured, H is the M×1 observation value vector, X is the M×(N+1) input feature matrix, β is the (N+1)×1 parameter vector, and the matrix H, matrix X and matrix β are as follows:

[0025]

[0026] Minimize the sum of squared error terms:

[0027]

[0028] In the formula, m is the label of the sample, ∈ T is the transpose of ∈, and the sum of squares of the error terms min β Find the partial derivative of parameter β and set the derivative to 0, and we get:

[0029]

[0030] Solving the equation yields:

[0031] β=(X T X) -1 X T h

[0032] Finally, the intercept term β0 and regression coefficients β1, β2, …, β n ; The comprehensive feature vector P of different pixels fusion (i,j) and regression coefficient β n Substitute the oil film thickness calculation model to calculate the oil film thickness h at different pixel points n .

[0033] Furthermore, in step 3, the oil spill concentration corresponding to each pixel is calculated according to the oil film thickness, which specifically includes:

[0034] According to the area S of the detection area corresponding to the fused image and the total number of pixels N in the fused image, the specific area s corresponding to each pixel in the detection area is calculated, s = S / N. The volume of water corresponding to each pixel is calculated based on the area s: According to the volume of water corresponding to the pixel and oil film thickness h n , calculate the oil spill concentration C corresponding to the pixel point obs,n :

[0035]

[0036] In the formula, C obs,n is the observed concentration corresponding to the nth pixel of the fused image.

[0037] Furthermore, in step 4, a convection-diffusion model of the target waters is established, and loss terms related to time and flow velocity are introduced into the convection-diffusion model, specifically including:

[0038] A convection-diffusion model is established, and the loss term R(C) is introduced:

[0039]

[0040] Where C is the oil spill concentration, is the velocity vector of the water body, D is the diffusion coefficient, t is the time, is the gradient operator, represents the rate of change of concentration C with time t, Represents the concentration C along the velocity vector The rate of change of direction, represents the diffusion effect of concentration C, represents the sum of the second-order derivatives of concentration in all directions;

[0041] The loss term R(C) is related to time t and flow velocity vector The formula is as follows:

[0042]

[0043] In the formula, α is the loss proportional constant, ε is the influence coefficient of flow velocity and time on the loss proportion, is the model of water flow velocity, t0 is the initial time of oil leakage, and (t-t0) is the time difference between the current time and the initial time; oil pollution and clean water are sampled in the target water area, and the simulated water flow velocity is used in the laboratory instead of Measure multiple sets of data including concentration C, observed time difference (t-t0) and final loss value, and put these data into the equation to calculate the proportional constant α and coefficient ε;

[0044] Establish a spatial coordinate system for the target water area, take any point in the target water area as the origin, take the water flow direction as the x-axis, take the direction from the center of the water body toward the shore as the y-axis, and take the vertical upward direction as the z-axis, and measure the flow velocity vector in the target water area. Among them, u x 、u y and u z are the actual measured flow velocities in the x-axis, y-axis and z-axis directions respectively. The convection-diffusion model is expressed as:

[0045]

[0046] Where D x , D y , D z are the diffusion coefficients in the x-, y-, and z-axis directions, respectively.

[0047] Furthermore, in step 4, the target water area is divided based on the spatial coordinate system, the convection-diffusion model is discretized, and the predicted concentration of each detection point in the spatial coordinate system is calculated through the discretized model, specifically including:

[0048] Assume that the initial oil spill position is (x0, y0, z0), the initial time is t0, the initial concentration of the current grid point is C0, the oil floats on the surface of the water, and the flow velocity in the z-axis direction is ignored, and we get:

[0049] C(x,y,t0)=C0δ(x-x0)δ(y-y0)

[0050] In the formula, δ represents the instantaneous release of oil pollution at the point (x0, y0) at time t0, and the initial concentration of oil pollution at this point is C0;

[0051] Discrete convection-diffusion model divides the space where the water body is located into several grid points according to the coordinates, and defines the discrete concentration at the grid points as Where (x, y) are the coordinates of the grid points, τ is the time step label, time is discretized, the time step is set to Δt, space is discretized, and the space step is set to Δx and Δy. The discretized convection-diffusion model is converted to:

[0052]

[0053] In the formula, is the concentration of the grid point (x, y) at time step τ, is the concentration of the grid point (x,y) at time step τ+1, is the discretization of the convection term, is the discretization of the diffusion term, It is the discretization of the loss term. Based on the initial oil spill position (x0, y0, z0) and the initial time t0, the coordinates of each detection point in the target waters spatial coordinate system are brought in, and the predicted concentration C of each detection point at different time steps is gradually calculated. mod .

[0054] Furthermore, step 5 specifically includes:

[0055] The predicted concentration C of each detection point calculated according to the convection-diffusion model mod , and the observed concentration C obtained in step 3 obs,n Establish the inversion model:

[0056]

[0057] In the formula, C obs is the observed concentration, C mod is the model predicted concentration, θ is the parameter that needs to be inverted, parameter θ∈{t0,x0,y0,z0}, θ * is the optimized parameter, parameter k is the label of the detection point in the target water area spatial coordinate system, K is the total number of detection points, where K = N, and all detection points are distributed in the detection area. The detection points correspond to the pixels in the fused image one by one. The gradient descent rule is used to calculate the gradient of Φ to the parameters (t0, x0, y0, z0) to update the parameters (t0, x0, y0, z0), as follows:

[0058]

[0059] in:

[0060]

[0061] In the formula, η is the learning rate. The iteration is stopped when the inversion model reaches the predetermined number of iterations. After the iteration is completed, the initial oil spill position is output. and initial time

[0062] Furthermore, step 1 specifically includes:

[0063] Use a drone equipped with an infrared sensor, an ultraviolet light source, and a fluorescent imaging device to collect infrared images and ultraviolet fluorescent images of the detection area at night or without obvious sunlight interference, wherein the infrared sensor is used to collect infrared images, which captures and records the infrared radiation intensity of each pixel point to generate an infrared image; the ultraviolet light source and the fluorescent imaging device are used to collect ultraviolet fluorescent images, aim the ultraviolet light source at the detection area, and the oil film in the detection area generates a fluorescent signal under the excitation of ultraviolet light, and the fluorescent imaging device captures and records the fluorescent signal to generate an ultraviolet fluorescent image;

[0064] The collected infrared image and ultraviolet fluorescence image are subjected to image denoising respectively, and the denoised infrared image and ultraviolet fluorescence image are subjected to image registration so that the pixels in the infrared image and ultraviolet fluorescence image are aligned in spatial position. After registration, the two images have the same spatial reference coordinate system.

[0065] The beneficial effects of this application are:

[0066] The technical solution in this application collects infrared images and ultraviolet fluorescence images of the target water area detection area in the same time period, and performs feature fusion on the infrared images and ultraviolet fluorescence images, and uses the fused comprehensive features to calculate the oil film thickness in the detection area. The multi-dimensional data collection and fusion method is used to calculate the oil film thickness, which can obtain more information from different images, thereby improving the accuracy and robustness of the oil film thickness in the detection area. Compared with the method of single-dimensional oil film thickness in the required technology, the technical solution in this application is more accurate, which is conducive to subsequent inversion analysis and reduces errors.

[0067] The technical solution in the present application introduces loss terms related to time and flow rate into the convection-diffusion model of the water body, which can more accurately describe the amount of oil lost due to reactions such as evaporation, degradation, sedimentation, and dissolution caused by the advancement of time and the movement of the water body, improve the accuracy of the convection-diffusion model, further reduce the error of the inversion analysis, and make the parameters of the oil leakage point and leakage time finally obtained more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] The advantages of the above and / or additional aspects of the present application will become apparent and easily understood in the description of the embodiments in conjunction with the following drawings, in which:

[0069] Figure 1 It is a schematic flowchart of a method for analyzing and inverting a water pollution incident situation according to an embodiment of the present application;

[0070] Figure 2 is a schematic block diagram of image fusion according to an embodiment of the present application. DETAILED DESCRIPTION

[0071] In order to more clearly understand the above-mentioned purposes, features and advantages of the present application, the present application is further described in detail below in conjunction with the accompanying drawings and specific implementation methods. It should be noted that the embodiments of the present application and the features in the embodiments can be combined with each other without conflict.

[0072] In the following description, many specific details are elaborated to facilitate a full understanding of the present application. However, the present application may also be implemented in other ways different from those described herein. Therefore, the protection scope of the present application is not limited to the specific embodiments disclosed below.

[0073] like Figure 1 As shown, this embodiment provides a method for inversion of situation analysis of water pollution incidents, the method comprising:

[0074] Step 1: Obtain an infrared image and an ultraviolet fluorescence image of the detection area, perform image preprocessing on them respectively, and register the preprocessed infrared image and ultraviolet fluorescence image so that the pixel positions of the two registered images are aligned.

[0075] Use a drone equipped with an infrared sensor, an ultraviolet light source, and a fluorescent imaging device to collect infrared images and ultraviolet fluorescent images of the detection area (the area selected in the target waters for the measurement experiment) in the target waters at night or without obvious sunlight interference. Among them, the infrared sensor is used to collect infrared images, which can capture and record the infrared radiation intensity of each pixel point to generate an infrared image; the ultraviolet light source and the fluorescent imaging device are used to collect ultraviolet fluorescent images. The ultraviolet light source is aimed at the detection area, and the oil film in the detection area generates a fluorescent signal under ultraviolet light excitation. The fluorescent imaging device captures and records the fluorescent signal to generate an ultraviolet fluorescent image.

[0076] In this embodiment, collecting images at night or without obvious interference from sunlight can reduce interference and make the collected images clearer. For infrared sensors, the background radiation is low and the ambient temperature is more uniform, which reduces interference from sunlight and other light sources. Usually, clearer infrared images can be obtained. For fluorescent imaging equipment, it is easier to distinguish the fluorescent signal generated by the ultraviolet light source, thereby ensuring the accuracy and reliability of detection.

[0077] The collected infrared images and ultraviolet fluorescence images are subjected to image denoising respectively to enhance the recognition accuracy. The infrared images can be denoised by Gaussian filtering, wavelet denoising and other methods, and the ultraviolet fluorescence images can be denoised by histogram equalization and beamforming algorithm. The denoised infrared images and ultraviolet fluorescence images are subjected to image registration to ensure that the infrared images and ultraviolet fluorescence images are aligned in space. After registration, the two images have the same spatial reference coordinate system. Registration can align the pixel positions of the infrared image and the ultraviolet fluorescence image, that is, the coordinates of the pixel points in the two images correspond one to one, which is convenient for subsequent fusion.

[0078] In this embodiment, feature point matching or other image registration algorithms can be used to implement it, such as using a feature point matching method for registration. Specifically, the SIFT algorithm is used to extract feature points of the two images, the nearest neighbor algorithm is used to match the feature point pairs of the two images, and the RANSAC algorithm is used to remove mismatches to obtain accurate feature point pairs. The least squares method or other optimization methods are used to calculate the affine transformation matrix between the two images, and the affine transformation matrix is ​​applied to one of the images and transformed to the coordinate system of the other image to complete the registration.

[0079] Step 2: extract the features and contrast information of the corresponding pixels in the infrared image and the ultraviolet fluorescence image, assign weights according to the contrast information of the corresponding pixels in the two images, and perform feature fusion on the corresponding pixels in the two images based on the weights.

[0080] Extract the radiation temperature T of each pixel in the infrared image IR (i, j), texture feature FIR (i, j) and thermal contrast C IR (i,j); extract the fluorescence intensity I of each pixel in the UV fluorescence image UV (i, j), texture feature F UV (i,j) and fluorescence contrast C UV (i, j); In this embodiment, these features and contrast information can be extracted using image processing algorithms in the prior art, such as Python (programming language) and OpenCV library (open source computer vision library that provides rich tools and algorithms for processing images and videos).

[0081] According to the extracted thermal contrast C IR (i,j) and fluorescence contrast C UV (i, j), assign weights to the two corresponding pixels in the infrared image and the ultraviolet fluorescence image. The weight of the pixel in the infrared image is ω IR (i, j), the weight of the pixel in the ultraviolet fluorescence image is ω UV (i,j), where for each pixel, ω IR (i,j) and ω UV The sum of (i, j) is 1, and the weight calculation formula is as follows:

[0082]

[0083] Where (i, j) is the coordinate label of the pixel point. Since the registered infrared image and ultraviolet fluorescence image have the same reference coordinate system and the pixels in the two images correspond one to one, the coordinate (i, j) can correspond to the coordinate of a pixel point in both images at the same time. The subscript IR represents the infrared image, and the subscript UV represents the ultraviolet fluorescence image.

[0084] Based on the weight distribution of different pixel positions, the pixels in the two images are fused in sequence to obtain a fused image, such as Figure 2 As shown in Figure 2, the fused image still uses the original reference coordinate system (i.e., the coordinate system used by the infrared image and the ultraviolet fluorescence image), and each pixel in the fused image contains a comprehensive feature vector. The comprehensive feature vector P fusion (i,j) is:

[0085] P fusion (i,j)=[ω IR ·T IR (i,j),ω UV I UV (i,j),ω IR ·F IR (i,j),ω UV ·F UV (i,j)]T

[0086] In the formula, [·] T It represents the transpose of the vector, and the subscript fusion represents the fusion feature. The fused image can comprehensively utilize the temperature information and texture information of the infrared image and the fluorescence intensity information and texture information of the ultraviolet fluorescence image to enhance the detection and recognition effect of the target.

[0087] Step 3: Based on the fused features, an oil film thickness calculation model is established. The oil film thickness corresponding to each pixel in the fused image is calculated using the oil film thickness calculation model, and the oil spill concentration corresponding to each pixel is calculated based on the oil film thickness.

[0088] Based on the comprehensive feature vector P fusion (i,j) Establish the oil film thickness calculation model:

[0089]

[0090] Where h is the oil film thickness, β0 is the intercept term, and β n is the regression coefficient, P fusion,n is the nth component of the comprehensive feature vector, ∈ is the error term, n is the label of the pixel in the fused image, and N is the total number of pixels in the fused image. Among them, n is the label of the pixel in the fused image, and the corresponding coordinates (i, j) are different only when n is different.

[0091] In this embodiment, in order to improve the accuracy of the oil film thickness calculation model, it is necessary to optimize the model parameters β0 and β n When the model was first established, the parameters β0 and β n The value of is the initial value assigned by the system. It is necessary to use optimization methods to continuously optimize these two parameters and finally get the optimal value. The model can be optimized using methods such as the least squares method (OLS) and the gradient descent method.

[0092] Optimize the oil film thickness calculation model and calculate the parameters β0 and β n The optimal value of , converts the above oil film thickness calculation model into:

[0093] H=Xβ+∈

[0094] Where H is the actual observation matrix of oil film thickness, X is the input feature matrix (including the column of intercept term with all 1s), β is the matrix including intercept term β0 and regression coefficient β n Assuming that a total of M samples are measured, H is the M×1 observation value vector, X is the M×(N+1) input feature matrix, β is the (N+1)×1 parameter vector, and the matrix H, matrix X and matrix β are as follows:

[0095]

[0096] In order to calculate β and optimize the model, we need to minimize the sum of squares of the error terms:

[0097]

[0098] In the formula, m is the label of the sample, ∈ T is the transpose of ∈, and the sum of squares of the above error terms min β Find the partial derivative of parameter β and set the derivative to 0, and we get:

[0099]

[0100] Solving the equation yields:

[0101] β=(X T X) -1 X T H

[0102] Finally, the intercept term β0 and regression coefficients β1, β2, …, β n , the value of the error term ∈ is equal to the difference between the actual observation value and the model prediction value. The error term can be used to evaluate the degree of fit of the model and further optimize the model. In this embodiment, the above process uses the least squares method to directly obtain a closed-form solution through matrix operations. This is a more direct optimization method. This formula minimizes the sum of squares of the error terms and uses matrix inversion to directly calculate the optimal parameter vector β.

[0103] Complete the intercept term β0 and regression coefficients β1, β2, …, β n After solving the problem, the comprehensive feature vector P of different pixels is fusion (i,j) and regression coefficient β n After the oil film thickness calculation model is brought in, the oil film thickness at different pixel points can be calculated.

[0104] According to the area S of the detection area corresponding to the fused image and the total number of pixels N in the fused image, the specific area s corresponding to each pixel in the detection area is calculated, s = S / N. The volume of water corresponding to each pixel is calculated based on the area s: In this embodiment, for more accurate inversion calculation, it is necessary to adjust the shooting height of the drone so that the area s is less than or equal to 1m 2 .

[0105] According to the volume of water corresponding to each pixel The oil film thickness value h corresponding to each pixel n , calculate the oil spill concentration C corresponding to each pixel in the detection area obs,n (Oil spill concentration refers to the amount of oil contained in a unit volume of water), as follows:

[0106]

[0107] In the formula, C obs,n is the oil spill concentration corresponding to the nth pixel of the fused image, that is, the observed concentration corresponding to the nth pixel. The coordinates of the pixel are (i, j). Each pixel in the fused image corresponds to a set of coordinates.

[0108] Step 4: Establish a convection-diffusion model for the target waters, introduce loss terms related to time and flow velocity into the convection-diffusion model, preset the initial oil spill location and initial oil spill time, divide the target waters based on the spatial coordinate system, discretize the convection-diffusion model and calculate the predicted concentration of each detection point in the spatial coordinate system through the discretized model.

[0109] Since the oil pollution will be dispersed during the flow of water, the oil pollution will gather into blocks and float on the water surface irregularly (most common oils will float on the water surface). Therefore, the oil spill concentration is needed to represent the oil pollution in the detection area. Most oils have very limited solubility in water and usually form an oil film floating on the water surface. In the target waters, some oil will evaporate, degrade, settle, dissolve and other reactions over time, which will further reduce the concentration of oil. Considering these reactions, it is necessary to introduce the loss term R(C) in the process of oil convection and diffusion, and the loss term R(C) should be related to the time and flow rate of the target waters. The following convection diffusion model is established:

[0110]

[0111] Where C is the oil spill concentration, is the velocity vector of the water body, D is the diffusion coefficient, R(C) is the loss term related to time and velocity used to describe reactions such as evaporation, degradation, sedimentation, and dissolution, t is the time, is the gradient operator, represents the rate of change of concentration C with time t, Represents the concentration C along the velocity vector The rate of change of direction, which indicates the change in concentration due to water flow, represents the diffusion effect of concentration C, that is, the concentration change caused by molecular diffusion, is the Laplace operator, which represents the sum of the second-order derivatives of the concentration in all directions.

[0112] The loss term R(C) is related to time t and flow velocity vector The specific formula is as follows:

[0113]

[0114] In the formula, α is the loss ratio constant, which represents the natural loss ratio of oil pollution in the absence of flow, and ε is the influence coefficient of flow velocity and time on the loss ratio. is the modulus of the water velocity, t0 is the initial time of the oil leakage, and (t-t0) is the time difference between the current time and the initial time. Among them, the loss proportional constant α and coefficient ε can be obtained through experiments. The oil pollution and clean water can be sampled in the target waters, and the simulated water velocity can be used in the laboratory instead of Measure multiple sets of data including concentration C, observed time difference (t-t0) and final loss value, and bring these data into the equation to calculate the proportional constant α and coefficient ε. This process can also be obtained by other conventional methods and will not be repeated here.

[0115] Establish a spatial coordinate system for the target water area, with any point in the target water area as the origin, the water flow direction as the x-axis, the direction from the center of the water body to the shore as the y-axis, and the vertical upward direction as the z-axis. Obtain the flow velocity vector in the target water area through actual measurement methods. Among them, u x 、u y and u z They are the actual measured flow velocities in the x-axis, y-axis and z-axis directions respectively; in this embodiment, instruments such as a flow meter and an ultrasonic Doppler flow meter can be used to measure the flow velocities in the x-axis, y-axis and z-axis directions in the spatial coordinate system of the target water area.

[0116] The above convection-diffusion model It can be expressed as:

[0117]

[0118] Where D x , D y , D z are the diffusion coefficients in the x-axis, y-axis, and z-axis directions respectively (which can be estimated from actual experiments).

[0119] According to the above-mentioned convection-diffusion model with the introduction of loss terms, when an initial oil spill position (x0, y0, z0) and an initial time t0 are preset, the oil spill concentration of each grid point in the spatial coordinate system of the target water area at different times can be calculated. In this embodiment, the spatial coordinate system of the target water area is established, and the target water area can be divided into a number of grid points according to the coordinate position. The spacing between any grid points is the same. Since oil pollution is detected in this embodiment and the observed concentration obtained is the oil spill concentration on the surface of the water surface, the flow velocity in the z-axis direction can be ignored during the calculation, and only the data in the x-axis and y-axis directions need to be calculated.

[0120] When using the convection-diffusion model with loss terms to calculate the concentration at different grid points, assume that the initial oil spill position is (x0, y0, z0), the initial time is t0, and the initial concentration of the current grid point is C0, that is, oil leakage occurs at the position with coordinates (x0, y0, z0) at time t0, and we get:

[0121] C(x,y,t0)=C0δ(x-x0)δ(y-y0)

[0122] Wherein, δ is the Dirac delta function, which means a function that has a value at point (x0, y0) and is zero at other points, that is, the oil pollution is released instantly at point (x0, y0) at time t0, and the initial concentration of pollutants at this point is C0. In this embodiment, the selection of the initial time and the initial oil spill location usually requires preliminary guesses based on prior knowledge or experimental data.

[0123] The convection-diffusion model is discretized, and the space where the water body is located is divided into several grid points according to the coordinates. The discrete concentration is defined at the grid points as Where (x, y) are the coordinates of the grid points, τ is the time step label, time is discretized, the time step is set to Δt, space is discretized, and the space step is set to Δx and Δy. The discretized convection-diffusion model is converted to:

[0124]

[0125] In the formula, is the concentration of the grid point (x, y) at time step τ, is the concentration of the grid point (x,y) at time step τ+1, is the discretization of the convection term, is the discretization of the diffusion term, is the discretization of the loss term. Through the above formula and the initial oil spill position (x0, y0, z0) and the initial time t0, the concentrations of other grid points at different time steps can be gradually calculated. After the coordinates of each detection point in the target waters spatial coordinate system are brought in, the predicted concentration C corresponding to each detection point can be calculated. mod .

[0126] Step 5: Establish an inversion model based on the observed concentration and the predicted concentration calculated by the convection-diffusion model, iteratively optimize the parameters of the inversion model, and obtain the optimal values ​​of the initial oil spill position and initial time.

[0127] The predicted concentration C of each detection point calculated according to the above diffusion model mod , and the observed concentration C obtained in step 3 obs,n Establish the inversion model:

[0128]

[0129] In the formula, C obs is the observed concentration, C mod is the model predicted concentration, θ is the parameter that needs to be inverted, parameter θ∈{t0,x0,y0,z0}, θ * is the optimized parameter, parameter That is, the preset initial time t0 and the guessed initial oil spill position (t0, x0, y0, z0), k is the label of the detection point in the spatial coordinate system of the target water area, K is the total number of detection points, where K = N, and all detection points are distributed in the detection area, and the detection points correspond to the pixels in the fused image one by one; the partial derivative of Φ with respect to the parameters (t0, x0, y0, z0) is calculated using the gradient descent rule, that is, the gradient, and the parameters (t0, x0, y0, z0) are updated as follows:

[0130]

[0131] in:

[0132]

[0133] In the formula, η is the learning rate. The iteration is stopped when the inversion model reaches the predetermined number of iterations. The gradient descent method can effectively optimize the inversion model to find the predicted concentration C mod As close as possible to the observed concentration C obs The best parameters After iterative optimization is completed, the initial oil spill position is output and initial time In this embodiment, after obtaining the predicted concentration at the current location, an inversion model is established by calculating the error between the predicted concentration and the observed concentration, and the initial oil spill location and time are continuously updated by optimization iteration, and the adjustment is repeated until the objective function inversion model Φ reaches a predetermined number of iterations, so that the error between the predicted concentration and the observed concentration is minimized, and the initial oil spill location is finally inverted. and initial time

[0134] In this embodiment, diesel is used as the most leaked oil, and a water tank simulating the target water area is used as an experimental water tank to provide power to keep the water flowing in the water tank. A certain amount of diesel is added at any time point and any position in the water tank. The time and position coordinates of the addition are used as the initial time and the initial oil spill position. A certain section of the water tank is selected as the detection area. At night or without obvious sunlight interference, infrared images and ultraviolet fluorescence images are taken using the methods of steps 1 to 3 in this application to obtain the observed concentration of each detection point. The predicted concentration of each detection point is calculated using the methods in steps 4-5. The initial time and initial oil spill position are inferred through the inversion model, and compared with the recorded initial time and initial oil spill position. The error is less than 5%.

[0135] In summary, the technical solution of the present application is described in detail above in combination with the accompanying drawings. The present application proposes a method for inversion of situation analysis of water pollution incidents, including: step 1, obtaining an infrared image and an ultraviolet fluorescence image of the detection area, and performing image preprocessing respectively, and aligning the infrared image and the ultraviolet fluorescence image that have completed the preprocessing, so that the pixel positions of the two images after the alignment are aligned; step 2, extracting the features and contrast information of the corresponding pixel points in the infrared image and the ultraviolet fluorescence image, assigning weights according to the contrast information of the corresponding pixel points in the two images, and performing feature fusion of the corresponding pixel points in the two images based on the weights; step 3, establishing an oil film thickness calculation model based on the fused features, and using the oil film thickness The computational model calculates the oil film thickness corresponding to each pixel in the fused image, and calculates the oil spill concentration corresponding to each pixel according to the oil film thickness as the observed concentration; Step 4, establishes a convection-diffusion model for the target waters, introduces loss terms related to time and flow velocity into the convection-diffusion model, presets the initial oil spill position and initial oil spill time, divides the target waters based on the spatial coordinate system, discretizes the convection-diffusion model, and calculates the predicted concentration of each detection point in the spatial coordinate system through the discretized model; Step 5, establishes an inversion model based on the observed concentration and the predicted concentration calculated by the convection-diffusion model, iteratively optimizes the parameters of the inversion model, and obtains the optimal values ​​of the initial oil spill position and initial time.

[0136] The technical solution in this application fuses the features of infrared images and ultraviolet fluorescence images, and introduces loss terms related to time and flow rate in the convection-diffusion model. This can solve the problem in the prior art that the subsequent inversion calculation cannot accurately infer the oil leakage point and leakage time due to inaccurate oil film thickness detection or inaccurate description of the amount of oil lost during water flow, thereby improving the inversion accuracy of leakage location and leakage time.

[0137] The steps in this application can be adjusted in order, combined, and deleted according to actual needs.

[0138] The units in the device of the present application can be combined, divided and deleted according to actual needs.

[0139] Although the present application is disclosed in detail with reference to the accompanying drawings, it should be understood that these descriptions are merely exemplary and are not intended to limit the application of the present application. The scope of protection of the present application is defined by the appended claims and may include various modifications, alterations and equivalents made to the invention without departing from the scope and spirit of the present application.

Claims

1. A method for inversion of water pollution incident situation analysis, characterized in that: The water pollution incident situation analysis inversion method comprises: Step 1, obtaining an infrared image and an ultraviolet fluorescence image of the detection area, and performing image preprocessing respectively, and registering the preprocessed infrared image and ultraviolet fluorescence image so that the pixel positions of the two registered images are aligned; Step 2, extracting features and contrast information of corresponding pixels in the infrared image and the ultraviolet fluorescence image, assigning weights according to the contrast information of corresponding pixels in the two images, and performing feature fusion on corresponding pixels in the two images based on the weights; Step 3, establishing an oil film thickness calculation model based on the fused features, using the oil film thickness calculation model to calculate the oil film thickness corresponding to each pixel in the fused image, and calculating the oil spill concentration corresponding to each pixel according to the oil film thickness as the observed concentration; Step 4, establishing a convection diffusion model for the target waters, introducing loss terms related to time and flow velocity into the convection diffusion model, presetting the initial oil spill location and initial oil spill time, dividing the target waters based on the spatial coordinate system, discretizing the convection diffusion model, and calculating the predicted concentration of each detection point in the spatial coordinate system through the discretized model; Step 5: establishing an inversion model based on the observed concentration and the predicted concentration calculated by the convection-diffusion model, iteratively optimizing the parameters of the inversion model, and obtaining the optimal values ​​of the initial oil spill position and initial time.

2. The method for inversion of water pollution incident situation analysis according to claim 1, characterized in that: The step 2 specifically includes: Extract the radiation temperature T of each pixel in the infrared image IR (i, j), texture feature F IR (i, j) and thermal contrast C IR (i, j), extract the fluorescence intensity I of each pixel in the ultraviolet fluorescence image UV (i, j), texture feature F UV (i,j) and fluorescence contrast C UV (i,j), according to the extracted thermal contrast C IR (i,j) and fluorescence contrast C UV (i, j) assigns weights to the two corresponding pixels in the infrared image and the ultraviolet fluorescence image: In the formula, ω IR (i,j) is the weight of the pixel in the infrared image, ω UV (i, j) is the weight of the pixel in the UV fluorescence image, (i, j) is the coordinate label of the pixel, ω IR (i,j) and ω UV The sum of (i, j) is 1, and the registered infrared image and ultraviolet fluorescence image have the same reference coordinate system; Based on the weight distribution of different pixel positions, the pixels in the two images are fused in turn to obtain a fused image. Each pixel in the fused image corresponds to a comprehensive feature vector P fusion (i,j): P fusion (i,j)=[ω IR ·T IR (i,j),ω UV ·I UV (i,j),ω IR ·F IR (i,j),ω UV ·F UV (i,j)] T In the formula, [·] T Represents the transpose of a vector.

3. The method for inversion of water pollution incident situation analysis according to claim 2, characterized in that: The step 3 of establishing the oil film thickness calculation model based on the fused features specifically includes: Based on the comprehensive feature vector P fusion (i,j) Establish the oil film thickness calculation model: Where h is the oil film thickness, β0 is the intercept term, and β n is the regression coefficient, P fusion,n is the nth component of the comprehensive feature vector, ∈ is the error term, n is the label of the pixel in the fused image, and N is the total number of pixels in the fused image; Optimize the oil film thickness calculation model and calculate the parameters β0 and β n The optimal value of , converts the oil film thickness calculation model into: H=Xβ+∈ Where H is the actual observation value matrix of oil film thickness, X is the input feature matrix, and β is the intercept term β0 and the regression coefficient β n Assuming that a total of M samples are measured, H is the M×1 observation value vector, X is the M×(N+1) input feature matrix, β is the (N+1)×1 parameter vector, and the matrix H, matrix X and matrix β are as follows: Minimize the sum of squared error terms: In the formula, m is the label of the sample, ∈ T is the transpose of ∈, and the sum of squares of the error terms min β Find the partial derivative of parameter β and set the derivative to 0, and we get: Solving the equation yields: β=(X T X) -1 X T H Finally, the intercept term β0 and regression coefficients β1, β2, …, β n ; The comprehensive feature vector P of different pixels fusion (i,j) and regression coefficient β n Bring in the oil film thickness calculation model to calculate the oil film thickness h at different pixel points n .

4. The method for inversion of water pollution incident situation analysis according to claim 3, characterized in that: Calculating the oil spill concentration corresponding to each pixel point according to the oil film thickness in step 3 specifically includes: According to the area S of the detection area corresponding to the fused image and the total number of pixels N in the fused image, the specific area s corresponding to each pixel in the detection area is calculated, s = S / N. The volume of water corresponding to each pixel is calculated according to the area s: According to the volume of water corresponding to the pixel and oil film thickness h n , calculate the oil spill concentration C corresponding to the pixel point obs,n : In the formula, C obs,n is the observed concentration corresponding to the nth pixel of the fused image.

5. The method for inversion of water pollution incident situation analysis according to claim 4, characterized in that: In step 4, a convection-diffusion model of the target waters is established, and loss terms related to time and flow velocity are introduced into the convection-diffusion model, specifically including: A convection-diffusion model is established, and the loss term R(C) is introduced: Where C is the oil spill concentration, is the velocity vector of the water body, D is the diffusion coefficient, t is the time, is the gradient operator, represents the rate of change of concentration C with time t, Represents the concentration C along the velocity vector The rate of change of direction, represents the diffusion effect of concentration C, represents the sum of the second-order derivatives of concentration in all directions; The loss term R(C) is related to time t and flow velocity vector The formula is as follows: In the formula, α is the loss proportional constant, ε is the influence coefficient of flow velocity and time on the loss proportion, is the model of water flow velocity, t0 is the initial time of oil leakage, and (t-t0) is the time difference between the current time and the initial time; oil pollution and clean water are sampled in the target water area, and the simulated water flow velocity is used in the laboratory instead of Measure multiple sets of data including concentration C, observed time difference (t-t0) and final loss value, and put these data into the equation to calculate the proportional constant α and coefficient ε; Establish a spatial coordinate system for the target water area, take any point in the target water area as the origin, take the water flow direction as the x-axis, take the direction from the center of the water body toward the shore as the y-axis, and take the vertical upward direction as the z-axis, and measure the flow velocity vector in the target water area. Among them, u x 、u y and u z are the actual measured flow velocities in the x-axis, y-axis and z-axis directions respectively. The convection-diffusion model is expressed as: Where D x , D y , D z are the diffusion coefficients in the x-, y-, and z-axis directions, respectively.

6. The method for inversion of water pollution incident situation analysis according to claim 5, characterized in that: In step 4, the target water area is divided based on the spatial coordinate system, the convection diffusion model is discretized, and the predicted concentration of each detection point in the spatial coordinate system is calculated through the discretized model, which specifically includes: Assume that the initial oil spill position is (x0, y0, z0), the initial time is t0, the initial concentration of the current grid point is C0, the oil floats on the surface of the water, and the flow velocity in the z-axis direction is ignored, and we get: C(x,y,t0)=C0δ(x-x0)δ(y-y0) In the formula, δ represents the instantaneous release of oil pollution at the point (x0, y0) at time t0, and the initial concentration of oil pollution at this point is C0; Discrete convection-diffusion model divides the space where the water body is located into several grid points according to the coordinates, and defines the discrete concentration at the grid points as Where (x, y) are the coordinates of the grid points, τ is the time step label, time is discretized, the time step is set to Δt, space is discretized, and the space step is set to Δx and Δy. The discretized convection-diffusion model is converted to: In the formula, is the concentration of the grid point (x, y) at time step τ, is the concentration of the grid point (x,y) at time step τ+1, is the discretization of the convection term, is the discretization of the diffusion term, is the discretization of the loss term. Based on the initial oil spill position (x0, y0, z0) and the initial time t0, the coordinates of each detection point in the target waters spatial coordinate system are brought in, and the predicted concentration C of each detection point at different time steps is gradually calculated. mod .

7. The method for inversion of water pollution incident situation analysis according to claim 6, characterized in that: The step 5 specifically includes: The predicted concentration C of each detection point calculated according to the convection-diffusion model mod , and the observed concentration C obtained in step 3 obs,n Establish the inversion model: In the formula, C obs is the observed concentration, C mod is the model predicted concentration, θ is the parameter that needs to be inverted, parameter θ∈{t0,x0,y0,z0}, θ * is the optimized parameter, parameter k is the label of the detection point in the target water area spatial coordinate system, K is the total number of detection points, where K = N, and all detection points are distributed in the detection area. The detection points correspond to the pixels in the fused image one by one. The gradient descent rule is used to calculate the gradient of Φ to the parameters (t0, x0, y0, z0) to update the parameters (t0, x0, y0, z0), as follows: in: In the formula, η is the learning rate. The iteration is stopped when the inversion model reaches the predetermined number of iterations. After the iteration is completed, the initial oil spill position is output. and initial time 8. The method for inversion of water pollution incident situation analysis according to claim 1, characterized in that: The step 1 specifically includes: Use a drone equipped with an infrared sensor, an ultraviolet light source, and a fluorescent imaging device to collect infrared images and ultraviolet fluorescent images of the detection area at night or without obvious sunlight interference, wherein the infrared sensor is used to collect infrared images, which captures and records the infrared radiation intensity of each pixel point to generate an infrared image; the ultraviolet light source and the fluorescent imaging device are used to collect ultraviolet fluorescent images, aim the ultraviolet light source at the detection area, and the oil film in the detection area generates a fluorescent signal under the excitation of ultraviolet light, and the fluorescent imaging device captures and records the fluorescent signal to generate an ultraviolet fluorescent image; The collected infrared image and ultraviolet fluorescence image are subjected to image denoising respectively, and the denoised infrared image and ultraviolet fluorescence image are subjected to image registration so that the pixels in the infrared image and ultraviolet fluorescence image are aligned in spatial position. After registration, the two images have the same spatial reference coordinate system.

Images