CN101900687B - Method for monitoring and early warning water bloom in small water area based on image processing - Google Patents

Abstract

An image processing-based monitoring and early warning method for algal blooms in small water areas relates to the monitoring and early warning methods for algal blooms in small water areas. The invention utilizes a camera and a computer to perform calculation through a program to realize monitoring and early warning of algae blooms in small water areas. The invention has the advantages of simple monitoring equipment and strong adaptability to the environment; the calculation amount of the algorithm is small, the speed is fast, and the hardware requirements are low, and it can be directly implanted into the existing water surface video monitoring system or developed into an embedded device or run directly on the computer. All can carry out early warning of water bloom quickly and effectively in real time; flexible and convenient, low cost, easy to maintain and so on. The invention can be widely used in monitoring and early warning of water blooms in rivers, upstream reservoirs of rivers, reservoirs of urban drinking water sources, river sections of urban water intake areas, or small water areas such as rivers, rivers, lake areas and landscape waters where algal blooms often occur.

Description

A monitoring and early warning method for algal blooms in small waters based on image processing

technical field

The invention belongs to the technical field of water environment monitoring, and in particular relates to small water areas (namely, rivers, upstream reservoirs of rivers, reservoirs of urban drinking water sources, river sections of urban water intake areas, rivers, rivers, lake areas and landscape water areas where algal blooms often occur) Water bloom monitoring and early warning method.

Background technique

In recent years, with the rapid development of industry and science and technology, a large number of industrial, agricultural and domestic wastes containing a large amount of nitrogen and phosphorus have been discharged into the water, resulting in eutrophication of fresh water in rivers, rivers, lakes and other waters and the outbreak of algal blooms ( Such as green algae bloom, cyanobacteria bloom, etc.) Water blooms have occurred in Taihu Lake, Dianchi Lake, Chaohu Lake, and Hongze Lake in China. The biggest harm caused by water blooms is: ① algae toxins indirectly affect human health or algae through the food chain The carcinogens produced directly threaten human health and survival, and pollute water sources; ②The filter devices of waterworks are filled with algae "blooms", affecting water intake; ③"Water blooms" floating on the water surface affect the landscape, and there are unpleasant Odor, pollute the environment, etc. It can be seen that algae blooms occur frequently, which not only causes damage to large areas of water environment, but also brings huge economic losses. Therefore, monitor fresh water sources such as source waters in the upper reaches of rivers, local lakes, and specific functional waters (urban water intake areas, drinking water sources, landscape waters, etc.), and monitor the occurrence, changes, and disaster trends of algal blooms in these water areas Early warning is very meaningful.

Existing water bloom monitoring and early warning methods, as published on August 6, 2008, the publication number is CN101236519A "buoy for cyanobacteria monitoring and cyanobacteria bloom early warning" patent, the disclosed buoy is composed of a buoy carrier, an instrument cabin, a monitoring sensor The assembly (composed of five different sensors), communication antennas, and solar panels. In application, a certain number of buoys are arranged at specific positions in the measured waters to monitor different depths of the water body, and the monitored data are transmitted to the monitoring center through the antenna for real-time early warning after data analysis. The main disadvantages of this patent are:

1. The limitations of monitoring are large. This patent only arranges a certain number of buoys at specific positions in the measured waters to monitor different depths in the water (that is, underwater monitoring). It cannot monitor large areas of water, nor can it monitor the spread of pollutants on the water surface Regularity, distribution range, pollution degree, etc.; buoys need to be in direct contact with water bodies during monitoring, and occupy a certain amount of water space, so they are not suitable for direct application in water areas with specific functions (such as drinking water sources, landscape water areas, etc.). Therefore, the limitations of buoy monitoring are relatively large.

2. High maintenance costs and high monitoring costs. In the practical application of this patent, specialized technicians are required to regularly go to the monitoring site to maintain the buoy, which is time-consuming and manpower-consuming, thereby greatly increasing the maintenance cost. In order to obtain accurate water quality information, expensive high-precision monitoring sensors must be used, and the sensors in this patent have a large quantity and various types, thereby increasing the monitoring cost.

3. Poor popularity. Because the cost of this patent is high, the structure is complex, and the use, operation and maintenance are highly specialized, it is not conducive to popularization and use in management departments, enterprises and scientific research units.

Contents of the invention

The purpose of the present invention is to address the shortcomings of the existing water bloom monitoring and early warning methods, and provide a small-scale water bloom monitoring and early warning method based on image processing, which has low cost, strong environmental adaptability, real-time acquisition of water surface pollution information and automatic It can identify algal blooms and give an alarm, which is convenient for popularization and application.

Mechanism of the present invention: the present invention is at first set up color prior model with the historical algae bloom region image of measured waters as sample figure, is about to convert each pixel of sample figure from three primary colors (RGB) space to color, brightness separation (HSV) Space, in order to lock the bloom area under various lighting conditions, the saturation (S) and brightness (V) components in the HSV color space are discarded, and only the hue (H) component is used. The HSV model is the closest to the human perception of color, and the color information is sensitive to the H (hue) component of this space, and the H (hue) component largely removes the influence of light, so the established color prior model Unaffected by changes in ambient light.

First establish the color prior model of the water bloom area, set {x′ _j } _{j=1, 2...n} as the n pixels of the water bloom area of the historical sample map, in the color model, when the H (hue) component The discretization value of is i (i=1, 2..., 360) and the hue probability is:

p={p _i } _{i=1, 2...m} ; $p_i=\frac{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\Σ}}}{\left|\right|x_j^{\′}\left|\right|}^2\mathrm{\δ}[b\left(x_j^{\′}\right)-i]}{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\Σ}}}\left({\left|\right|x_j^{\′}\left|\right|}^2\right)};$ $\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{p}_i=1---\left(1\right)$

In formula (1): p is the color model, p _i is the probability of the hue value i, δ is the δ function, and the function b(x′ _j ) is the index of the space R ² →{1, 2…m}, that is, in The index of the pixel at position x′ _j to the histogram quantization feature space.

Then use the continuous adaptive mean shift (Camshift) algorithm to detect the bloom color block, the continuous adaptive mean shift algorithm is evolved from the mean shift (MeanShift) algorithm, the process of the mean shift (MeanShift) algorithm is through the kernel A recursive process of updating the center of the kernel G by sampling the mean translation vector M _{h on G(x)} (ie, the estimate of the density gradient on the kernel K). When the mean shift algorithm is applied to a continuous sequence, the mean shift vector is continuously calculated from the initial center of the position of the kernel G, and the target position is updated iteratively until it converges to the optimal matching point.

M _{h, G(x)} (that is, the estimation of the density gradient on the kernel K) is calculated as follows:

${\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&dtri;</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; C</span>}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; G</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

Among them: h is the band width, C is the normalization constant, it can be seen that the sampling mean translation vector on the kernel G is the estimation of the density gradient on the kernel K. is the multivariate density estimate of band width h over kernel K:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; no</span>}{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; d</span>}}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\frac{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The set {x′ _j } _{j=1, 2...n} is n points in the d-dimensional Euclidean space R ^d , and k(x) represents the kernel function of the pixel point:

Among them: c _d is the volume of the d-dimensional unit sphere.

The continuous adaptive mean shift algorithm is to perform directional projection and mean shift calculation in the processing area of the video image, take the centroid and area of the search window of the current frame as the initial value of the search window of the next frame, and then iterate to achieve the goal. Detection tracking.

Then use the improved continuous adaptive mean shift algorithm to filter the interference color blocks, that is, before applying the Camshift algorithm, firstly corrode and expand the probability distribution map to remove the noise blocks, and then detect the inner and outer contours of the target color block to find the target color block. Holes in the blocks are filled. However, after the above preprocessing, there may still be interference color blocks, so it is necessary to improve the adaptability of the Camshift algorithm when it is applied to the detection of water blooms.

It usually takes about 15 days from the appearance of algae on the water surface to the full-scale outbreak of the water bloom. Therefore, in a short period of time, the change of the green area on the water surface is slow. In comparison, the time interval between two adjacent frames of images is very Short (tens of milliseconds), then the change in the area of the target color block between frames is very small. Based on this change rule, we can filter out the interference color blocks by predicting the area of the feature.

Suppose the area of the color block in the current frame image is S _n , and the predicted area of the next frame is S′ _n+1 , then

S' _n+1 =S _n +A (5)

Among them, the constant A is an empirical value (set according to the specific situation), which represents the increment of the area of the target color block between frames. According to formula (5), the area of the target color block in the next frame image can be predicted, if it satisfies

|S _n+1 -S′ _n+1 |≤ξ (6)

S _n+1 represents the actual measured area of the color block in the next frame, and ξ is the allowable error range, indicating that the currently detected color block conforms to the change rule of the color block area between frames; if |S _n+1 -S′ _{n+ 1} | When it exceeds the allowable error range, it is considered that the area of the color block actually detected in the next frame does not conform to the change law of the area of the bloom area between frames, and the excess part is defined as the interference color block and filtered out.

The technical scheme that realizes the object of the present invention is: a kind of small-scale water area algae bloom monitoring and early warning method based on image processing, utilizes camera and computer, by program, establishes priori color model according to historical data, uses improved continuous self-adaptive mean drift ( Camshift) algorithm filters out the interference color blocks and detects the algae bloom area, analyzes the disaster trend through the change of the area area, and establishes an early warning model similar to meteorological disasters to monitor and warn of algal blooms. Specific steps are as follows:

(1) Establish a color prior model

Firstly, based on the sample image or video data of the historical algal bloom area of the measured water area, a color prior model is established, namely:

Let {x′ _j } _{j=1, 2...n} be n pixels of the sample image or video data, convert each pixel from RGB (three primary colors) space to HSV (color, bright separation) space, and convert H (hue ) component's discretization value is i (i=1, 2..., 360), and the range [0, 360] of the H (hue) component is scaled to [0, 255], so that the value of the value range can be Represented by a byte (byte). The hue probability of the H component is calculated by the computer according to the following formula:

p={p _i } _i=12...m ; $p_i=\frac{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}{\left|\right|x_j^{\&prime;}\left|\right|}^2\mathrm{\&delta;}[b\left(x_j^{\&prime;}\right)-i]}{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}\left({\left|\right|x_j^{\&prime;}\left|\right|}^2\right)};$ $\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_i=1---\left(1\right)$

Where: p is the color model, p _i is the probability of the hue value i, δ is the δ function, and the function b(x′ _j ) is the index of the space R ² →{1, 2…m}, that is, at the position x′ _j The pixels of are indexed into the histogram quantization feature space.

(2) Detect the bloom area

①Monitoring of small waters

After step (1) is completed, in small waters (that is, upstream reservoirs of rivers or rivers or reservoirs of urban drinking water sources or river sections of urban water intake areas or rivers, rivers, lakes or landscape waters where algal blooms often occur) A camera is set up on the shore, and is connected with a computer through a video transmission line. The specific position, quantity and height of the camera installation shall be determined according to the specific conditions of the water area to be measured. It is used to monitor the water surface conditions of the measured water area (the most significant change due to the outbreak of algal blooms in the water area is the change of the visual characteristics of the water surface) and capture video images of the water surface to achieve the purpose of real-time monitoring of water surface images.

② Reverse projection processing

After the (2)-① step is completed, carry out back projection processing (BackProjection) to each frame of video image input in the (2)-① step, that is, for each pixel in the video image processing area, query the image by computer The degree of matching between the pixel and the color prior model established in step (1) (i.e., the H-component tone probability model) can obtain the probability that the pixel is the target pixel (the probability of other areas outside this area is 0, and each The value of a pixel becomes a discretized measure of the possibility of the target color information appearing here, the greater the possibility of appearing here, the larger the value of the pixel, and vice versa). After the above processing, the target color back-projection map of each frame image is obtained.

③Detect and segment the water bloom color blocks and filter out the interfering color blocks

After step (2)-② is completed, detect the bloom color block and filter out the interference color block processing, namely:

First, corrode and expand the back projection image obtained in step (2)-②, detect the inner and outer contours of the target color block and find the holes in the target color block for filling and other preprocessing, so that part of the interfering color blocks can be removed. Avoid noise interference to a certain extent.

Then detect and segment the water bloom color block area, that is, calculate the mean translation vector M _{h, G} (x) continuously to the back projection image after the above-mentioned preprocessing through the computer according to the formula:

${\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="2"\; label="h"\; filenames="HSA00000171985400011.png"\; state="\{\{state\}\}">h</figure-callout></span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&dtri;</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; C</span>}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; G</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

为核K上带宽度为h的多变量密度估计：Where: h is the band width, C is the normalization constant,

is the multivariate density estimate of band width h over kernel K:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; no</span>}{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; d</span>}}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\frac{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The set {x′ _j } _{j=1, 2...n} is n points in the d-dimensional Euclidean space R ^d , and k(x) represents the kernel function of the pixel point:

Among them: c _d is the volume of the d-dimensional unit sphere. It is used to iteratively update the position of the target color block until it converges to the optimal matching point, and the area of the matching point is the target color block (that is, the water bloom color block area) in the current frame image.

Judgment: When the target color block is not detected, it is considered that there is no bloom on the water surface, and return to step (2)-① to obtain the next frame to continue detection; when the target color block is detected, the current frame is the most The centroid and area of the optimal matching point are used as the initial value of the search window in the next frame, and the iterative calculation is continued to find the optimal matching point. This loop iterates until the end of the program, and the target color block of each frame image is detected.

Then filter the interference color block, that is, calculate the area of the target color block in the current frame and the next frame obtained in the previous step, and use the area of the target color block in the current frame to predict the image of the next frame according to the inter-frame variation law of the area of the bloom area The area of the target color block in , that is, calculate the difference between the actual detection area of the color block in the next frame and the predicted area of the next frame (i.e. the redundant part), and judge according to the following formula:

|S _n+1 -S′ _n+1 |≤ξ (5)

Among them: S _n+1 is the actual measured area of the color block in the next frame, S′ _n+1 is the predicted area of the next frame, and ξ is the allowable error range.

When |S _n+1 -S′ _n+1 |≤ξ, it is considered that the area of the color patch actually detected in the next frame conforms to the inter-frame variation law of the area of the bloom area and takes it as the detection result; when |S _{n +1} -S′ _n+1 |>ξ, that is, when the difference exceeds the allowable error range, it is considered that the area of the color block actually detected in the next frame does not conform to the change law of the area of the bloom area between frames, and the redundant part Define it as an interference color block and filter it out.

After the above-mentioned processing, the water bloom color block after filtering out the interference color block is detected and segmented, and the pixel area of the water bloom color block in the image is obtained. Then the water bloom color block area in each frame of image is calibrated and the detection result is output.

(3) Water bloom warning

① Calculate the actual area of the algae bloom area:

After step (2) is completed, first calculate the percentage of the pixel area of the water bloom color block obtained in step (2) in the image to the total area of the field of view pixels, and multiply it by the actual total area of the field of view to obtain the actual area of the water bloom area . Then through the camera calibration, the actual total area of the field of view is calculated, that is, the reference object with a known area is placed in the water area under the camera, and the ratio of its pixel area in the image to the total area of the field of view pixels is calculated, and the actual area of the known reference object is used. Divide the area by the ratio to obtain the actual total area of the field of view, and output the obtained actual area of the algae bloom area to the computer software interface, so that the user can intuitively understand the actual area of the algae bloom area.

②Establish an early warning model:

After the (3)-1. step is completed, according to the actual area of the algae bloom area obtained by the (3)-1. step, the mean value of the actual area of the algae bloom area in each time interval is calculated at an interval of 6 to 10 hours. Then, analyze the change law of the mean value of the actual area by computer to establish a linear regression forecasting model to predict the mean value of the actual area of the next time interval, and calculate the percentage of the mean value of the actual area in the total area of the water area, and then represent it with this percentage Based on the degree of algae bloom disaster, the early warning levels of algae bloom outbreaks are divided into four levels: blue, yellow, orange and red, and the corresponding early warning intervals are general, heavy, serious and particularly severe, and judged : when the percentage of the representative algal bloom disaster degree is less than the early warning value (the size of the early warning value is determined according to the specific situation and application requirements of the measured water area, and the present invention gets 5%), then return to return (2)-1. step acquisition Continue to detect in the next frame; when the percentage of the representative algal bloom disaster degree is greater than the warning value, the warning will be given according to the divided warning level. The specific classification of algal bloom outbreak warning levels is as follows:

General blue warning (i.e. sporadic algae bloom): sporadic accumulation of algal blooms, algae biodensity in the main water area is less than 30 million/L, and the area of algae bloom is greater than or equal to 5% of the total area of the water body.

Heavy yellow warning (localized algae bloom): Algae gather in local waters, the density of algae in the main water area is between 30 million and 50 million/L, and the area of algae bloom is greater than or equal to 10% of the total area of the water body.

Severe orange warning (i.e. regional algal bloom): When regional algal bloom occurs, the density of algae in the main water area is between 50 and 80 million/L, and the area of algae bloom is greater than or equal to 40% of the total area of the water body.

Especially serious red warning (i.e. comprehensive algae bloom): Algae blooms are all-round outbreaks, the algae density in the main water area is greater than 80 million/L, and the area of algal blooms is greater than or equal to 60% of the total area of the water body.

After the present invention adopts above-mentioned technical scheme, mainly have following effect:

1. It has good popularity. The invention provides a fast and effective water bloom monitoring and early warning method, which can automatically identify the water bloom area and calculate the actual area of the area in real time, and can perform early warning according to the early warning level so as to notify relevant departments to understand the degree of water bloom disaster in real time and take corresponding measures. governance measures. In practical applications, it can be directly implanted into the traditional water surface video monitoring system or developed into an embedded device to realize the monitoring and early warning of algal blooms. It has the characteristics of intelligence and easy operation, and is easy to popularize in management departments, enterprises and scientific research units. use.

2. Low cost. The algorithm adopted in the present invention has small calculation amount and high speed, and has low requirements on hardware. In practical applications, whether it is run directly on a computer or developed into an embedded device, it can quickly and effectively realize algal bloom monitoring and early warning in real time, which greatly saves monitoring costs and is flexible and convenient.

3. Strong ability to adapt to the environment. The present invention uses a camera to monitor the water surface and removes the influence of light to a large extent, which can effectively avoid the noise interference caused by the complexity of the water surface environment and the change of ambient light, so that it can well adapt to environmental changes, so the monitoring equipment is simple, Low cost, easy maintenance, and effective real-time monitoring.

The present invention can be directly implanted into a traditional water surface video monitoring system or developed into an embedded device, and is widely used in rivers, upstream reservoirs of rivers, reservoirs of urban drinking water sources, river sections in urban water intake areas or rivers where algal blooms often occur, Water bloom monitoring and early warning of small waters such as rivers, lakes and landscape waters.

Description of drawings

Fig. 1 is a program flow diagram of the method of the present invention.

Detailed ways

The present invention will be further described below in combination with specific embodiments.

Example

As shown in Figure 1, a method for monitoring and early warning of water blooms in small waters based on image processing is used to monitor and early warning of water blooms in the Pengxi River (the first-level tributary of the Yangtze River located in the heart of the Three Gorges Reservoir Area) in Kaixian County, Chongqing City. The total area of the water area is about 63,000 square meters. The specific steps are as follows:

(1) Establish a color prior model

Firstly, based on the sample image or video data of the historical algal bloom area of the measured water area, a color prior model is established, namely:

Let {x′ _j } _{j=1, 2...n} be n pixels of the sample image or video data, convert each pixel from RGB (three primary colors) space to HSV (color, bright separation) space, and convert H (hue ) component's discretization value is i (i=1, 2..., 360), and the range [0, 360] of the H (hue) component is scaled to [0, 255], so that the value of the value range can be Represented by a byte (byte). The hue probability of the H component is calculated by the computer according to the following formula:

p={p _i } _{i=1, 2...m} ; $p_i=\frac{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}{\left|\right|x_j^{\&prime;}\left|\right|}^2\mathrm{\&delta;}[b\left(x_j^{\&prime;}\right)-i]}{\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}\left({\left|\right|x_j^{\&prime;}\left|\right|}^2\right)};$ $\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_i=1---\left(1\right)$

Where: p is the color model, p _i is the probability of the hue value i, δ is the δ function, and the function b(x′ _j ) is the index of the space R ² →{1, 2…m}, that is, at the position x′ _j The pixels of are indexed into the histogram quantization feature space.

(2) Detect the bloom area

①Monitoring of small waters

After step (1) is completed, set up a camera on the water tower on the bank of the small water area of the Pengxi River Reservoir in Kaixian County, Chongqing City. The camera field of view covers the entire water area and is connected to the computer through a video transmission line. The specific position, quantity and height of the camera installation shall be determined according to the specific conditions of the water area to be measured. It is used to monitor the water surface conditions of the measured water area (the most significant change due to the outbreak of algal blooms in the water area is the change of the visual characteristics of the water surface) and capture video images of the water surface to achieve the purpose of real-time monitoring of water surface images.

② Reverse projection processing

After the (2)-① step is completed, carry out back projection processing (BackProjection) to each frame of video image input in the (2)-① step, that is, for each pixel in the video image processing area, query the image by computer The degree of matching between the pixel and the color prior model established in step (1) (i.e. the H-component tone probability model) can obtain the probability that the pixel is the target pixel (the probability of other areas outside this area is 0, and each The value of a pixel becomes a discretized measure of the possibility of the target color information appearing here, the greater the possibility of appearing here, the larger the value of the pixel, and vice versa). After the above processing, the target color back-projection map of each frame image is obtained.

③Detect and segment the water bloom color blocks and filter out the interfering color blocks

After step (2)-② is completed, detect the bloom color block and filter out the interference color block processing, namely:

First, corrode and expand the back projection image obtained in step (2)-②, detect the inner and outer contours of the target color block and find the holes in the target color block for filling and other preprocessing, so that part of the interfering color blocks can be removed. Avoid noise interference to a certain extent.

Then detect and segment the water bloom color block area, that is, calculate the mean translation vector M _{h, G(x)} continuously for the back projection image after the above preprocessing through the computer according to the following formula:

${\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="2"\; label="h"\; filenames="HSA00000171985400011.png"\; state="\{\{state\}\}">h</figure-callout></span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&dtri;</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; C</span>}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; G</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

为核K上带宽度为h的多变量密度估计：Where: h is the band width, C is the normalization constant,

is the multivariate density estimate of band width h over kernel K:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; no</span>}{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; d</span>}}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\frac{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The set {x′ _j } _{j=1, 2...n} is n points in the d-dimensional Euclidean space R ^d , and k(x) represents the kernel function of the pixel point:

Among them: c _d is the volume of the d-dimensional unit sphere. It is used to iteratively update the position of the target color block until it converges to the optimal matching point, and the area of the matching point is the target color block (that is, the water bloom color block area) in the current frame image.

Judgment: When the target color block is not detected, it is considered that there is no bloom on the water surface, and return to step (2)-① to obtain the next frame to continue detection; when the target color block is detected, the current frame is the most The centroid and area of the optimal matching point are used as the initial value of the search window in the next frame, and the iterative calculation is continued to find the optimal matching point. This loop iterates until the end of the program, and the target color block of each frame image is detected.

Then filter the interference color block, that is, calculate the area of the target color block in the current frame and the next frame obtained in the previous step, and use the area of the target color block in the current frame to predict the image of the next frame according to the inter-frame variation law of the area of the bloom area The area of the target color block in , that is, calculate the difference between the actual detection area of the color block in the next frame and the predicted area of the next frame (i.e. the redundant part), and judge according to the following formula:

|S _n+1 -S′ _n+1 |≤ξ (5)

Among them: S _n+1 is the actual measured area of the color block in the next frame, S′ _n+1 is the predicted area of the next frame, and ξ is the allowable error range.

When |S _n+1 -S′ _n+1 |≤ξ, it is considered that the area of the color patch actually detected in the next frame conforms to the inter-frame variation law of the area of the bloom area and takes it as the detection result; when |S _{n +1} -S′ _n+1 |>ξ, that is, when the difference exceeds the allowable error range, it is considered that the area of the color block actually detected in the next frame does not conform to the change law of the area of the bloom area between frames, and the redundant part Define it as an interference color block and filter it out.

After the above-mentioned processing, the water bloom color block after filtering out the interference color block is detected and segmented, and the pixel area of the water bloom color block in the image is obtained. Then the water bloom color block area in each frame of image is calibrated and the detection result is output.

(3) Water bloom warning

① Calculate the actual area of the algae bloom area:

After step (2) is completed, first calculate the percentage of the pixel area of the water bloom color block obtained in step (2) in the image to the total area of the field of view pixels, and multiply it by the actual total area of the field of view to obtain the actual area of the water bloom area . Then through the camera calibration, the actual total area of the field of view is calculated, that is, the reference object with a known area is placed in the water area under the camera, and the ratio of its pixel area in the image to the total area of the field of view pixels is calculated, and the actual area of the known reference object is used. Divide the area by the ratio to obtain the actual total area of the field of view, and output the obtained actual area of the algae bloom area to the computer software interface, so that the user can intuitively understand the actual area of the algae bloom area.

②Establish an early warning model:

After the (3)-① step is completed, according to the actual area of the algae bloom area obtained in the (3)-① step, the average value of the actual area of the algae bloom area in each time interval is calculated at an interval of 8 hours, and then passed The computer analyzes the change law of the mean value of the actual area to establish a linear regression prediction model to predict the mean value of the actual area in the next time interval, and calculate the percentage of the mean value of the actual area in the total area of the water area, and then represent the bloom with this percentage Based on the degree of disaster, the early warning levels of algae bloom outbreaks are divided into four levels: blue, yellow, orange and red. When the percentage of the representative algal bloom disaster degree is less than the early warning value (the size of the early warning value is determined according to the specific conditions and application requirements of the measured water area, the present invention gets 5%), then return to the return (2)-1. step to obtain the next step The frame continues to be detected; when the percentage of the represented algal bloom disaster degree is greater than the warning value, the warning is given according to the divided warning level. The specific classification of algal bloom outbreak warning levels is as follows:

General blue warning (i.e. sporadic algae bloom): sporadic accumulation of algal blooms, algae biodensity in the main water area is less than 30 million/L, and the area of algae bloom is greater than or equal to 5% of the total area of the water body.

Heavy yellow warning (localized algae bloom): Algae gather in local waters, the density of algae in the main water area is between 30 million and 50 million/L, and the area of algae bloom is greater than or equal to 10% of the total area of the water body.

Severe orange warning (i.e. regional algal bloom): When regional algal bloom occurs, the density of algae in the main water area is between 50 and 80 million/L, and the area of algae bloom is greater than or equal to 40% of the total area of the water body.

Especially serious red warning (i.e. comprehensive algae bloom): Algae blooms are all-round outbreaks, the algae density in the main water area is greater than 80 million/L, and the area of algal blooms is greater than or equal to 60% of the total area of the water body.

After testing the image processing-based algal bloom monitoring and early warning method in small water areas, the following conclusions can be drawn:

①During the monitoring process, algae bloom occurred on the water surface of the measured water area (Pengxi River Reservoir, Kaixian County, Chongqing City). The method of the present invention can detect the area of the algae bloom color block in real time and calculate the actual area of the area, and can accurately carry out early warning according to the early warning level , in the early stage of algae bloom outbreak, Chongqing Municipal Environmental Protection Bureau was able to grasp the algae bloom disaster situation in the water area in time and notify relevant departments to take corresponding measures to control it (salvage in sporadic algae bloom), effectively cooperating with disaster prevention and mitigation work.

②The method of the present invention can quickly and effectively realize algae bloom monitoring and early warning by only using a computer and a camera, which greatly saves monitoring costs, is convenient for users to use, and has low maintenance costs. The monitoring effect can meet the actual needs of the Chongqing Municipal Environmental Protection Bureau , so as to facilitate the promotion and application.

③ The present invention adopts the improved continuous self-adaptive mean shift algorithm to detect and divide the water bloom color blocks and filter the interference color blocks, which can effectively suppress the noise interference caused by the complexity of the water surface background and the change of ambient light. The algorithm has a small amount of calculation and high speed, which is the basis for ensuring the real-time monitoring and early warning of algal blooms.

The above conclusions illustrate that the monitoring and early warning of water blooms in small water areas based on image processing realized by the method of the present invention can realize real-time monitoring and early warning of water blooms. Low, strong adaptability to the environment, easy to popularize and apply, so the present invention can be applied in actual projects.

Claims (1)

1. the method for monitoring and early warning water bloom in small water area based on Flame Image Process is characterized in that utilizing camera and computing machine, calculates through program, and its concrete steps are following:

(1) sets up the color prior model

At first, set up the color prior model, that is: according to the historical wawter bloom in tested waters regional sample figure or video data

If x ' _j} _{J=1,2 ... N}Be n pixel of sample figure or video data, with each pixel from the three primary colors space conversion to look, show the score from the space, be i with the discretize value of tone component, i.e. i=1,2..., 360, zoom to [0,255] to the scope of tone component [0,360]; Be calculated as follows the tone probability of tone component through computing machine:

p＝{p _i} _{i＝1，2…m}； $p_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left|\right|x_j^{\&prime;}\left|\right|}^2\mathrm{\&delta;}[b\left(x_j^{\&prime;}\right)-i]}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left({\left|\right|x_j^{\&prime;}\left|\right|}^2\right)};$ $\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_i=1---\left(1\right)$

Wherein: p is a color model, p _iBe that tone value is the probability of i, δ is the δ function, function b (x ' _j) be space R ²→ { 1,2 ... The index of m} promptly is positioned at position x ' _jPixel to the index in histogram quantization characteristic space;

(2) the wawter bloom zone is detected

1. small-sized waters is monitored

(1) step was set up camera on the bank in small-sized waters after accomplishing, and was connected with computing machine through video transmission line; The particular location that camera sets up, quantity and height are confirmed according to the concrete condition in tested waters;

2. back projection is handled

After the completion of (2)-1. step; Each frame video image to the input of (2)-1. step carries out the back projection processing; Promptly to each pixel in the video image processing region; The matching degree of the color prior model of setting up through this pixel of computer inquery and (1) step just obtains the color of object back projection figure of every two field picture;

3. detect and cut apart the wawter bloom color lump and filter the interference color lump

(2)-2. step was detected the wawter bloom color lump and disturbs color lump to handle with filtering, that is: after accomplishing

The back projection figure that earlier (2)-2. step was obtained carries out burn into and expands, and detects the inside and outside contour of target color lump and finds the cavity in the target color lump to fill pre-service;

Detect again and cut apart the wawter bloom color block areas, promptly press following formula constantly to the back projection figure computation of mean values translation vector M after the above-mentioned pre-service of process through computing machine _{H, G (x)}:

$M_{h,G\left(x\right)}=\frac{h^2\&dtri;{\hat{f}}_K\left(x\right)}{2/C{\hat{f}}_G\left(x\right)}---\left(2\right)$

Wherein: h is a bandwidth; C is a normaliztion constant,

be that the last bandwidth of nuclear K is the multivariate density Estimation of h:

${\hat{f}}_K\left(x\right)=\frac{1}{n{h}^d}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}k\left({\left|\right|\frac{x-x_i}{h}\left|\right|}^2\right)---\left(3\right)$

Set x ' _j} _{J=1,2 ... N}Be d dimension Euclidean space R ^dN point, k (x) representes the kernel function of this pixel:

Wherein: c _dBe d dimension unit spheroid volume; In order to the position of iteration renewal target color lump, up to the Optimum Matching point that converges on, the zone of this match point is exactly the target color lump in the current frame image;

Judge: when not detecting the target color lump, then returned for (2)-1. step and obtain next frame and proceed to detect; When detecting the target color lump; Then with the barycenter of present frame Optimum Matching point and area initial value as the next frame search window; Proceed iterative computation and find the Optimum Matching point, so till loop iteration to the EOP (end of program), just detect the target color lump that is partitioned into every two field picture;

Filter then and disturb color lump; The area of the present frame that the last step of i.e. calculating obtains and the target color lump of next frame; Interframe Changing Pattern according to the wawter bloom region area; Utilize the target color lump area of present frame to predict the target color lump area in the next frame image, promptly calculate the difference between next frame color lump actual detected area and the next frame prediction area, judge according to following formula:

|S _n+1-S′ _n+1|≤ξ (5)

Wherein: S _N+1Be the actual measurement area of next frame color lump, S ' _N+1Be the prediction area of next frame, the error range of ξ for allowing;

When | S _N+1-S ' _N+1| during≤ξ, then think the next frame actual detected to the color lump area meet the interframe Changing Pattern of wawter bloom region area and with it as testing result; When | S _N+1-S ' _N+1|＞ξ, when promptly difference exceeds in the error range of permission, then think the next frame actual detected to the color lump area do not meet the Changing Pattern of interframe wawter bloom region area, redundance is defined as disturbs color lump and it is filtered; Detection is partitioned into and filters the wawter bloom color lump that disturbs behind the color lump fixed and output testing result of rower of going forward side by side;

(3) bloom prealarming

1. calculate wawter bloom zone real area:

(2) step calculated (2) earlier and goes on foot the number percent that the elemental area of wawter bloom color lump in image that obtains accounts for the visual field pixel total area after accomplishing, and the actual total area that multiply by the visual field just can obtain wawter bloom zone real area; Pass through camera calibration again; Calculate the actual total area in the visual field; The object of reference that is about to known area places the position, waters under the camera; Calculate the ratio that its elemental area in image accounts for the visual field pixel total area, just obtain the actual total area in the visual field divided by this ratio, and the wawter bloom zone real area that obtains is outputed on the computer software interface with the real area of this known object of reference;

2. set up Early-warning Model:

After the completion of (3)-1. step; Go on foot the wawter bloom zone real area that obtains according to (3)-1.; Be at interval by 6～10 hours earlier; Calculate the average of the wawter bloom zone real area in each time interval, set up the one-variable linear regression forecast model through the Change in Mean rule of this real area of Computer Analysis again, predict the average of the real area in the next time interval; And the average of calculating this real area accounts for the number percent of the waters total area; Wawter bloom disaster degree with the representative of this number percent is a foundation then, and the advanced warning grade of breakout of water bloom is divided into blueness, yellow, orange and red four grades, corresponding early warning interval be respectively generally, heavier, serious and four grades of especially severe; And judge:, then turn back to and returned for (2)-1. step and obtain next frame and proceed to detect when the number percent of wawter bloom disaster degree of representative during less than early warning value; When the number percent of wawter bloom disaster degree of representative during, then carry out early warning according to the advanced warning grade of dividing greater than early warning value;

The concrete division of breakout of water bloom advanced warning grade is following:

General blue early warning: wawter bloom is sporadicly gathered, and main waters district algae bio density is less than 3,000 ten thousand/L, and the wawter bloom area is more than or equal to 5% of the water body total area;

Heavier yellow early warning: algae gathers in local waters, and main waters district algae bio density is between 3000～5,000 ten thousand/L, and the wawter bloom area is more than or equal to 10% of the water body total area;

Serious orange early warning: when generation area property wawter bloom, main waters district algae bio density is between 5000～8,000 ten thousand/L, and the wawter bloom area is more than or equal to 40% of the water body total area;

The especially severe red early warning: wawter bloom is comprehensive to be broken out, and the algae bio density in main waters district is greater than 8,000 ten thousand/L, and the wawter bloom area is more than or equal to 60% of the water body total area.

Images