CN107590816B - Water body information fitting method based on remote sensing image - Google Patents
Water body information fitting method based on remote sensing image Download PDFInfo
- Publication number
- CN107590816B CN107590816B CN201710806824.1A CN201710806824A CN107590816B CN 107590816 B CN107590816 B CN 107590816B CN 201710806824 A CN201710806824 A CN 201710806824A CN 107590816 B CN107590816 B CN 107590816B
- Authority
- CN
- China
- Prior art keywords
- water body
- fitting
- image
- function
- reflectivity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
The invention relates to a water body information extraction and fitting method, in particular to a water body information fitting method based on a remote sensing image, and provides a water body information fitting method based on a remote sensing image, aiming at solving the defects that the existing conventional measurement is difficult to master the water body change and water quality change conditions in time, and the conventional measurement possibly cannot find some pollution sources and the characteristics of the pollution sources, and the method comprises the following steps: processing the remote sensing image by using a water body index method to obtain a processed image; performing two-dimensional Otsu threshold segmentation on the processed image to obtain a threshold segmentation result; taking intersection of the threshold segmentation result and the processed image, and determining the light reflectivity of the image after the intersection is taken; selecting a function model, and calculating the dissolved oxygen and the permanganate by using the light reflectivity and the function model respectively to obtain the fitting degree of the dissolved oxygen and the fitting degree of the permanganate; and respectively selecting the fitting with the optimal fitting degree according to a preset standard for fitting. The method is suitable for fitting the water body information.
Description
Technical Field
The invention relates to a water body information extraction and fitting method, in particular to a water body information fitting method based on remote sensing images.
Background
The conventional water quality analysis method is field sampling, namely continuously monitoring and recording the types and the concentrations of pollutants in a water body for a long time, analyzing the water pollution change trend according to the recorded data, and evaluating the water quality condition in a grading manner. In order to observe different areas at different time, multiple sampling is needed, in order to deeply understand the water pollution change trend of the observed area, corresponding processing needs to be carried out on the sampled data points, such as abnormal points removal, and the like, so that huge manpower, material resources and financial resources are consumed; along with population growth, urban area extension, industrial and agricultural pollution and the like, the water pollution degree is continuously worsened, water pollution areas are continuously diffused, and partial water samples after sampling cannot represent the distribution characteristics of water pollution and the water quality change conditions of all the areas; surface water is easily affected by human factors such as climate change, land coverage and the like, and the water body change and the water quality change condition are difficult to master in time through conventional measurement due to the characteristics of high dynamic state, uneven distribution and the like of the water body; with the increasing variety of pollution sources, some pollution sources and the characteristics of the pollution sources may not be found through conventional measurement.
Disclosure of Invention
The invention aims to provide a water body information fitting method based on remote sensing images, which aims to solve the problems that the existing conventional measurement is difficult to master the water body change and the water quality change condition in time, and some pollution sources and pollution source characteristics can not be found in the conventional measurement.
A water body information fitting method based on remote sensing images comprises the following steps:
step one, processing a remote sensing image by using a water body index method to obtain a processed image;
secondly, performing two-dimensional Otsu threshold segmentation on the processed image to obtain a threshold segmentation result;
thirdly, taking intersection of the threshold segmentation result and the processed image, and determining the light reflectivity of the image after the intersection is taken;
sequentially selecting function models from a preset function set, and calculating the dissolved oxygen and the permanganate by using the light reflectivity and the function models respectively to obtain the fitting degree of the dissolved oxygen and the fitting degree of the permanganate corresponding to each function model;
and fifthly, selecting the fitting degree with the optimal fitting degree from the fitting degree of the dissolved oxygen and the fitting degree of the permanganate according to preset standards to perform fitting.
The invention has the beneficial effects that: 1. the water body change and the water quality change can be mastered in time; 2. the pollution source and the characteristics of the pollution source are easier to extract; 3. the invention has low cost and is easy to realize.
Drawings
FIG. 1 is a flow chart of a water body information fitting method based on remote sensing images according to the invention;
FIG. 2 is a partition of a two-dimensional histogram from a threshold segmentation vector of the present invention;
fig. 3 is a comparison graph of extraction results of an improved normalized difference water body index method and a two-dimensional Otsu threshold segmentation method, wherein fig. 3a is an image after RGB synthesis, fig. 3b is an image after processing by the improved normalized difference water body index method, and fig. 3c is an image after processing by the two-dimensional Otsu threshold segmentation method;
FIG. 4 is a graph of a fifth embodiment permanganate fit;
FIG. 5 is a graph of an inverse model of dissolved oxygen in accordance with a fifth embodiment;
fig. 6 is a remote sensing image after registration in the seventh embodiment, wherein fig. 6a is an image of registration of 5 months 2014 and 5 months 2015 of pinkish; FIG. 6b is an image of the water region registered at 5 months 2015 and 6 months 2015; FIG. 6c is an image of the water region registered at 5 months 2015 and 7 months 2015; FIG. 6d is an image of the water region registered at 5 months 2015 and 10 months 2015;
fig. 7 is a graph of the change region of the water area in fig. 6, wherein fig. 7a is a graph of change of 5 months in 015 versus 5 months in 2014, and fig. 7b is a graph of change of 6 months in 2015 versus 5 months in 2015; FIG. 7c is a graph of the change of month 7 in 2015 to month 5 in 2015, and FIG. 7d is a graph of the change of month 10 in 2015 to month 5 in 2015;
FIG. 8 is a graph of the results of the 2016 inversion, wherein FIG. 8a is a graph of the results of the 2016 inversion at 9 months, and FIG. 8b is a graph of the results of the 2016 inversion at 10 months;
fig. 9 is a graph of inversion results of 2015, wherein fig. 9a is a graph of inversion results of 2015 at 4 months, fig. 9b is a graph of inversion results of 2015 at 7 months, and fig. 9c is a graph of inversion results of 2015 at 10 months;
FIG. 10 is a diagram showing the distribution of permanganate content in the Songhua river water area;
FIG. 11 is a diagram showing the distribution of the dissolved oxygen content in Songhua river water.
Detailed Description
The first embodiment is as follows: the water body information fitting method based on the remote sensing image comprises the following steps:
step one, processing the remote sensing image by using a water body index method to obtain a processed image.
And step two, performing two-dimensional Otsu threshold segmentation on the processed image to obtain a threshold segmentation result.
In particular, obtaining water body information from different remote sensing images is often the first step in water body related research. The traditional water body extraction method belongs to the water body extraction based on spectral characteristics in principle, namely, the water body and other spectral characteristics with obvious ground object categories (the water body has high absorption rate in a near infrared band and low reflectivity in a visible light band) are utilized, and the water body index method is adopted to process multispectral remote sensing images, so that the water body is enhanced, non-water bodies are restrained, and a proper threshold value is manually selected to extract the water body on the basis.
According to the method, the two-dimensional Otsu threshold segmentation is applied to water body extraction on the basis of improving the normalized difference water body index, the method automatically selects the appropriate threshold according to the two-dimensional histogram of the image, and compared with the manual threshold selection, the method is higher in applicability and higher in precision.
And thirdly, taking intersection of the threshold segmentation result and the processed image, and determining the light reflectivity of the image after the intersection is taken.
And step four, sequentially selecting function models from a preset function set, and calculating the dissolved oxygen and the permanganate by using the light reflectivity and the function models respectively to obtain the fitting degree of the dissolved oxygen and the fitting degree of the permanganate corresponding to each function model.
The water quality monitoring usually selects one or more wave bands as dependent variables according to the correlation between the spectral reflectivity and the actually measured water quality parameter concentration, and performs linear or nonlinear regression to obtain a regression model. The data required for water quality monitoring typically includes the spectral reflectance of the remote sensing image and the data of the on-site sampling point. The time of the remote sensing data is preferably kept consistent with that of the field sampling data so as to ensure timeliness; secondly, the sampling point data is not too small to ensure the accuracy and applicability of the model. The ammonia nitrogen content, the oxygen content and the permanganate content in water have certain correlation with the single-band reflectivity or the band combination reflectivity. Whereas in their sensitive bands the correlation coefficient is larger.
Based on the thought, the method adopts nearly 4 years Landsat 8 remote sensing images and measurement data of the origin-causing sites of Heilongjiang province (east longitude 124 degrees 59 '20', north latitude 45 degrees 28 '15'), analyzes the correlation between the spectral reflectivity and the measured data, finds out the sensitive wave band, analyzes the fitting condition of an inversion model of permanganate and the dissolved oxygen in water in the forms of exponential function, power function and polynomial function, and selects a proper inversion model according to the square Sum of the fitted Errors (Sum of the Squared Errors), the deterministic coefficient (R-square), the Root-Mean-square Error (Root Mean Squared Error) and the corrected deterministic coefficient (Adjusted R-square). If the Adjusted R-square is negative, the fitting is unsuccessful, and the closer to 1, the better; the R-square interval is [ -1,1], the closer the absolute value is to 1, the better; the closer the SSE and RMSE are to 0, the better.
And fifthly, selecting the fitting degree with the optimal fitting degree from the fitting degree of the dissolved oxygen and the fitting degree of the permanganate according to preset standards to perform fitting.
The invention utilizes the remote sensing technology to monitor water quality, has the advantages of large monitoring area, high speed, convenience for long-term monitoring and the like, and in order to more accurately master the water body and the water quality change trend, the invention applies two-dimensional Otsu threshold segmentation to extract the water body on the improved normalized difference water body index, establishes an inversion model of water quality parameters according to the correlation between the spectral reflectivity of the remote sensing image and the field reference data on the basis of water body extraction, and analyzes the water quality of the water body in the Harbin region through the established model.
The second embodiment is as follows: the first difference between the present embodiment and the specific embodiment is: in the first step, the calculation method of the normalized difference water body index in the water body index method comprises the following steps:
MNDWI＝(Green-SWIR1)/(Green+SWIR1)
wherein MNDWI is normalized difference water body index, and SWIR1 is short wave infrared band.
Particularly, among various water body extraction methods, the most widely used is the water body index method. It highlights the body of water of interest by enhancing the difference between the body of water and adjacent pixels, while suppressing most of the noise.
On the basis of an improved normalized difference water body index (NDWI), the Xuequauu finds that the short-wave infrared band is used for replacing a near-infrared band, the contrast between a water body and a background can be enhanced more than the NDWI, and the improved normalized difference water body index (MNDWI) is provided:
MNDWI＝(Green-SWIR1)/(Green+SWIR1) (0-1)
other steps and parameters are the same as those in the first embodiment.
The third concrete implementation mode: the present embodiment differs from the first or second embodiment in that: in the second step, the specific process of Otsu threshold segmentation is as follows:
according to the formula
F (x, y) is less than or equal to s and g (x, y) is more than t; and pixels where f (x, y) > s and g (x, y) ≦ t are ignored and set to 0 or 1.
The two-dimensional Otsu thresholding method is to divide a two-dimensional histogram region with gradient gray levels and apply the conventional Otsu thresholding twice on two projection histograms in order to separate the region of interest from the background. The method has strong robustness to noise and high operation speed.
Assuming that the size of a gray image f (x, y) (x is more than or equal to 1 and less than or equal to M, y is more than or equal to 1 and less than or equal to N) is M × N, calculating the average gray value of N × M neighborhood of each pixel point f (x, y) in the image to obtain a smooth g (x, y), and the gray levels of the two are 0, 1. Where an n × m neighborhood may be characterized by the following template:
l is typically 255, n and m are numbers greater than 1, which may not be equal. There are many templates for selecting the neighborhood in the current relevant document, and an 8-neighborhood template or a 4-neighborhood template is usually adopted. In the invention, 8 neighborhood templates are selected.
For each pixel in the image, we can obtain a pair of (i, j), where i is the original gray level appearing in f (x, y) and j is the neighborhood average gray level appearing in g (x, y). Let cijThe probability of (i, j) is expressed as the joint probability
Wherein
j＝0，1，...，L-1。
Given an arbitrary threshold vector (s, t). The areas a and D represent the object and the background, respectively. Regions B and C represent edges and noise. Let two classes C0And C1Respectively representing an object and a background; omegaiAnd muiRespectively represent class CiProbability and mean vector of (2).
C0And C1The probability of (d) can be expressed as:
C0and C1The vector of mean values of (a) may be expressed as:
the total average vector of the 2D histogram is:
in most cases, the probability of being far from the diagonal is negligible. It is easy to find: omega0+ω 11, and μT≈ω0μ0+ω1μ1。
The inter-class dispersion matrix is defined as:
the traces of the discrete matrix are represented as:
thus, the optimal threshold value may be expressed as:
the conditions are satisfied: s is less than or equal to f (x, y)*And g (x, y) > t*(ii) a And f (x, y) > s*And g (x, y) is less than or equal to t*The pixel of (1) is ignored and set to 0 or 1. Fig. 2 shows the partitioning of a threshold segmentation vector into two-dimensional histograms.
The water body information extraction result of the present embodiment is analyzed by the following experiment:
in order to compare whether the extraction result is accurate, color synthesis is a better method, and comparison is convenient. The invention adopts Landsat 81T grade product obtained from USGS for experiment. In order to enhance the water body, Green wave bands are regenerated by using weighting processing of Green wave bands, Red wave bands and NIR wave bands. Namely R: red, G: (Green + Red + NIR)/3, B: green, obtaining a new waveband, synthesizing the three wavebands, displaying RGB, and obtaining a synthetic result as shown in figure 3 a. In the Landsat 8 remote sensing image, Green, Red and NIR respectively correspond to a third wave band, a fourth wave band and a fifth wave band.
The extraction results of the improved normalized difference water body index method and the two-dimensional Otsu threshold segmentation method are respectively shown in fig. 3b and fig. 3 c.
In order to compare the two extraction precisions more intuitively, the extraction result of each method is compared with the water body information in the truth-value chart, and according to the definition of the correct detection rate P (TA), the false alarm rate P (FA) and the undetected rate P (MA):
and calculating the correct detection rate, the false alarm rate and the omission factor, and analyzing the precision, wherein the analysis precision is shown in a table 1-1.
TABLE 1-1 extraction accuracy
Other steps and parameters are the same as those in the first or second embodiment.
The fourth concrete implementation mode: the difference between this embodiment mode and one of the first to third embodiment modes is: in the fourth step, the preset function set includes: a first order exponential function, a second order exponential function, a first order polynomial, a second order polynomial, a third order polynomial, a first order power function, and a second order power function.
In general, there are two ways to build a mathematical model: the first is a mechanistic analysis method, and the other is a test analysis method. Models established by mechanistic analysis methods are generally of clear significance. In the test analysis, under the condition that the internal mechanism of a research object cannot be directly analyzed, the research object is regarded as a black box, input and output data passing through the black box are measured, and a model with the best fitting effect with the measured data is selected from a certain type of determined models by using a statistical analysis method on the basis.
The establishment of the inversion model in the invention belongs to a test analysis method. Firstly, a commonly used inversion model function form is given, secondly, the correlation between the spectral reflectivity and the reference data concentration is analyzed, one or two wave bands with the strongest correlation are selected for regression analysis, and a model with the best fitting effect is selected according to the fitted mean square error and the like.
In statistics, regression analysis refers to a statistical analysis method that analyzes interdependent relationships, i.e., correlations, between two or more variables. According to the quantity of variables, the method can be divided into univariate regression analysis and multivariate regression analysis; according to the relationship between independent variables and dependent variables, there are classified into linear regression analysis and nonlinear regression analysis.
The water quality parameter inversion model function form selected in the invention is shown in the table 2-1.
And analyzing the correlation according to the spectral reflectivity and the water quality parameter concentration, selecting a wave band or a plurality of wave bands with strong correlation for regression analysis, and establishing an inversion model. The spectral reflectance and the reference data concentration are shown in tables 2-2 and 2-3.
TABLE 2-2 spectral reflectance
Table 2-3 reference data concentrations
Other steps and parameters are the same as those in one of the first to third embodiments.
The fifth concrete implementation mode: the difference between this embodiment and one of the first to fourth embodiments is: in the fourth step, the fitting degree specifically includes: sum of squares of errors, deterministic coefficients, root mean square error, modified deterministic coefficients.
In the fourth embodiment, the fitting degree of the curve is further determined by the sum of squares of errors, the deterministic coefficient, the root mean square error, and the corrected deterministic coefficient, and the function model with the best fitting degree is selected as the final choice.
The specific process of one embodiment is as follows:
and (4) carrying out correlation statistics on the reflectivity and the reference data to find the degree of correlation between the inverted high-manganese acid salt index and the spectral reflectivity. The correlation coefficients are shown in tables 2-4.
TABLE 2-4 correlation coefficient between reflectance and reference data
It can be seen that the first few bands have strong correlation with permanganate, and are fitted by using the reflectivity of the second band, the reflectivity of the third band, the product of the reflectivity of the second band and the third band, and the sum of the reflectivity of the second band and the third band as independent variables, and the corresponding permanganate index content as dependent variables, and using function forms such as power function, exponential function, first-order polynomial, second-order polynomial, and the like. Due to the fact that the function form order is too high, overfitting can be caused, the function is changed violently, and the change trend of the actual water quality parameters is not met; too low an order of the functional form may result in some points being fitted unsuccessfully and with reduced accuracy. Therefore, the argument is (ρ)2+ρ3) Namely the sum of the reflectivities of the second and third wave bands, the dependent variable is the content of permanganate, and the function form is a second-order exponential form, so that the fitting effect is better. The fitted curve of permanganate is shown in fig. 4.
Functional form: y is a.ebx+c·edx
Parameter values: a is 13.17; -8.742; c is 1.321; d-4.726
The degree of fit: SSE 0.1098; r-square 0.9596; adjusted R-square 0.9191;
RMSE＝0.1914
and (4) carrying out correlation statistics on the reflectivity and the reference data to find the correlation degree between the dissolved oxygen and the spectral reflectivity. The correlation coefficients are shown in tables 2-5.
TABLE 2-5 correlation coefficient between reflectance and reference data
It can be seen that the correlation between the dissolved oxygen and each band is not strong, and 0.5942 is the strongest correlation with the first band. The spectral reflectivities of a plurality of wave bands can be combined according to a certain mathematical relationship and then are subjected to inversion operation with dissolved oxygen, so that the reliability is greatly increased. Fitting by using the first waveband reflectivity, the second waveband reflectivity, the sum of the first waveband reflectivity and the second waveband reflectivity, the product of the first waveband reflectivity and the second waveband reflectivity as independent variables, using the corresponding dissolved oxygen content as a dependent variable, and using function forms such as power function, exponential function, first-order polynomial, second-order polynomial and the like. The experimental result shows that the fitting degree of the function type of the third-order polynomial is best by taking the product of the reflectivities of the first wave band and the second wave band as an independent variable, the dissolved oxygen content as a dependent variable and the function type as the third-order polynomial. The fitted curve of the dissolved oxygen amount is shown in fig. 5.
Functional form: y ═ p1·x3+p2·x2+p3·x+p4
Parameter values: p is a radical of1＝1195000；p2＝-91810；p3＝2183；p4＝-8.265
The degree of fit: SSE 5.294; r-square 0.7923; adjusted R-square 0.5846; RMSE ═ 1.328
Other steps and parameters are the same as in one of the first to fourth embodiments.
The sixth specific implementation mode: the difference between this embodiment and one of the first to fifth embodiments is:
fifthly, establishing four evaluation indexes, specifically:
the closer the corrected decision coefficient is to 1, the better the fitting degree is;
the closer the absolute value of the determinant coefficient is to 1, the better the fitting degree is;
the closer the sum of squared errors is to 0, the better the degree of fit;
the closer the root mean square error is to 0, the better the fitting degree is;
step two, calculating the absolute value, the error sum of squares and the root mean square error of the corrected decision coefficient, the absolute value, the error sum of squares and the root mean square error which correspond to each function model for the dissolved oxygen and the permanganate content;
step three, counting how many evaluation indexes of each function model are higher than those of other function models; and selecting the function model which meets the conditions and has the most evaluation indexes.
That is, in the embodiment, 4 indexes are used to determine whether the fitting degree is good, and a function that meets the indexes is selected to perform fitting. And if the number of the qualified indexes of the plurality of functions is the largest at the same time, selecting all the functions, fitting all the functions, and selecting the image which is most consistent with the actual terrain condition from the images obtained by fitting as a final selection.
The embodiment provides a method for selecting a function model with optimal fitting degree, and in the actual process, in order to ensure that the fitting result is more consistent with the actual situation, final verification and approval can be manually carried out, and the model which is more consistent with the objective conditions of the water body in an intuitive way is screened out.
Other steps and parameters are the same as those in one of the first to fifth embodiments.
The seventh embodiment: the difference between this embodiment and one of the first to sixth embodiments is: before the first step is executed, the method further comprises a multi-temporal remote sensing image registration step, and specifically comprises the following steps:
a, performing Harris angular point feature extraction on a multi-temporal remote sensing image;
b, eliminating mismatching and mismatching vectors by using a normalized cross-correlation matching algorithm to obtain a registered image;
and C, processing the registered image as a remote sensing image in the first step.
The embodiment aims at water quality analysis of multi-temporal remote sensing images, and aims to perform rocker image registration, and registered images can be input as one step of images.
Specifically, in order to visually recognize the change of the water body region, the change detection of the multi-temporal remote sensing images is required, and the registration between the remote sensing images is required before the change detection. Firstly, Harris angular point feature extraction is carried out, then rough matching is carried out by utilizing a normalized cross-correlation matching algorithm (NCC), mismatching and mismatching vectors are eliminated, and registration errors are calculated based on gray scale correlation coefficients, so that registered images are obtained.
(1) Harris angular point detection algorithm
Harris proposed a Harris corner detection algorithm in 1988, and if considering the difference in contribution weights of different points within a window, the corner response function can be written as:
wherein the content of the first and second substances,
is a two-dimensional gaussian window function. If the diversity of the offset directions is considered, the corner response function of the Harris algorithm can be obtained from the formula (4-1) by using a first-order Taylor approximation:
wherein, IxAnd IyIs the first derivative of the pixel (u, v) in the horizontal and vertical directions:
three values are calculated for the pixel points:
then a gradient matrix of pixel points can be derived
The Corner Response Function (CRF) expression of Harris is derived from this:
CRF(u,v)＝det(M)-k(trace(M))2＝(AB-C2)-k(A+B)2 (0-14)
wherein det (M) AB-C2Trace (M) is usually 0.04 to 0.06, where k is a constant. When the CRF value of the target is larger than or equal to the threshold value, the pixel point is judged as a corner point.
(2) NCC matching algorithm
The Normalized Cross Correlation method (NCC) matching algorithm is a classical statistical matching algorithm, and determines the degree of matching by calculating the Cross Correlation values of the template and the match.
According to the definition of vector dot product:
a·b＝|a|·|b|·cosθ (0-15)
if the two vectors are similar, the directions of the two vectors are the same, and the included angle is theta, so that the similarity of the two vectors can be judged according to the value of cos theta. Generalizing it into two dimensions, then
Wherein R (u, v) is a normalized correlation coefficient of the position point (u, v); n is a radical of1×N2Is the matching template size; x is the number ofi+u,j+v，yi,jThe gray values at (i + u, j + v), (i, j) in the two frames to be matched are respectively. The larger the value of R (u, v), the higher the similarity between the two is demonstrated.
The result after registration is shown in fig. 6.
And intercepting the common area in each registered group, extracting the water body, and calculating the change area. Since performing water extraction is equivalent to binarizing the gray level, the change region is calculated by the following formula:
change＝abs(imt1-imt2) (0-17)
the imt1 and the imt2 are respectively a first remote sensing image and a second remote sensing image, and the first remote sensing image and the second remote sensing image are subjected to subtraction and absolute value taking to obtain a change area. The results of the detection are shown in FIG. 7.
It can be found by referring to the data on the spot that pinkish river is classified according to the dissolved oxygen, and generally only water bodies of class i and class ii are available, so that the dissolved oxygen is classified more finely and displayed for more detailed display analysis.
3-5mg/L for IV water bodies and 5-6mg/L for III water bodies, and uniformly displaying the III water bodies and the IV water bodies in green; for the II-class water body, the water body is more finely divided into 6-6.5mg/L, 6.5-7mg/L and 7-7.5mg/L which are respectively displayed by blue, yellow and pink colors; for class I water, the water content is more finely divided into 7.5-10mg/L and more than 10mg/L, and the colors are respectively displayed by orange and cyan. Wherein, the lower the ammonia nitrogen content, the worse the water quality, from class I to class V to poor class V (the dissolved oxygen concentration is less than 2mg/L), the worse the water quality. For 2016, the inversion results are shown in FIG. 8.
And (3) analyzing an inversion result: the concentration of the permanganate in the water body of the Binhahar zone of the Songhua river is concentrated in the range of the IV class water body, namely 6-10 mg/L. After inversion, the permanganate content was classified and the regions of different permanganate content were displayed in different colors.
Uniformly displaying the class I water body, the class II water body and the class III water body in a green way; for IV water bodies, the IV water bodies are more finely divided into 6-6.5mg/L, 6.5-7mg/L, 7-7.5mg/L and 7.5-10mg/L which are respectively displayed by blue, yellow, orange and pink colors; the class V water body and the inferior class V water body (the concentration of permanganate is more than 15mg/L) are uniformly expressed by cyan. Wherein, the lower the permanganate content is, the better the water quality is, from type I to type V to poor type V, the water quality is increasingly poor. The inversion results are shown in fig. 9, for 2015 as an example.
The experimental results were analyzed:
the water body ratios of different permanganate contents and oxygen contents are calculated by counting the water bodies and are shown in figures 10 and 11,
other steps and parameters are the same as those in one of the first to sixth embodiments.
From the analysis in fig. 10, the permanganate content in the watershed of the piny Jianghahar segment generally belongs to the IV-class water body, a few parts of the watershed belong to the III-class water body, and the I-class water body, the II-class water body, the V-class water body and the inferior V-class water body do not exist. In IV-class water bodies, the content of permanganate mostly falls within the range of 6.5-7 mg/L. Most of permanganate is concentrated in the range of 6-6.5mg/L in about 9 and 10 months every year, and the water quality is improved; in summer, the content of high-manganese acid salt is increased to about 7mg/L, and the water quality is deteriorated. And the water quality at the edge of the river is worse than the water quality inside the water body. From the analysis of permanganate, in the water pollution of Songhua Jianghha Bin river basin, the permanganate pollution is serious, and the water quality is poor.
According to the analysis of FIG. 11, the dissolved oxygen of the waterfront of Haerbin Songhuang is concentrated in the range of 7.5-10mg/L, and belongs to I-class water body. The few parts of the water bodies belong to the class V water body, namely the concentration is in the range of 3-5 mg/L. From the oxygen content angle analysis, the water body of the Songhua Jianghua Harbin basin has no pollution and has better water quality.
As can be predicted from FIGS. 10 and 11, the permanganate content in the water in the waterfront of the Harbin basin of Songhua river will float around 7mg/L in the 6-10 months of the year, and will drop to around 6mg/L in the rest months, but in general, the water still belongs to IV-class water, and measures should be taken to improve the water quality; in 10 or 11 months each year, the maximum dissolved oxygen in the water body of the Harbin basin of Songhua river is 7.5-10mg/L, the water quality is better, and the rest months can float in the range of about 7mg/L, but in general, the water body belongs to I or II water bodies, and the water quality is better.
The present invention is capable of other embodiments and its several details are capable of modifications in various obvious respects, all without departing from the spirit and scope of the present invention.
Claims (2)
1. A water body information fitting method based on remote sensing images is characterized by comprising the following steps:
step one, processing a remote sensing image by using a water body index method to obtain a processed image;
the water body index method refers to a normalized difference water body index, and the calculation method of the normalized difference water body index comprises the following steps:
MNDWI＝(Green-SWIR1)/(Green+SWIR1)
wherein MNDWI is normalized difference water body index, SWIR1 is short wave infrared band;
secondly, performing two-dimensional Otsu threshold segmentation on the processed image to obtain a threshold segmentation result;
the specific process of Otsu threshold segmentation is as follows:
according to the formula
F (x, y) is less than or equal to s*And g (x, y)>t*(ii) a And f (x, y)>s*And g (x, y) is less than or equal to t*Is ignored and set to 0 or 1;
assuming that the size of a gray image f (x, y) (1 ≦ x ≦ M,1 ≦ y ≦ N) is mxn, the average gray value of N × M neighborhood is calculated for each pixel point f (x, y) in the image to obtain a smooth g (x, y), both gray levels are 0,1,. L, where N × M neighborhood is characterized by the following template:
taking L as 255, taking n and m as numbers larger than 1, and selecting a neighborhood template as 8;
for each pixel in the image, obtaining a pair of (i, j), where i is the original gray level appearing in f (x, y) and j is the neighborhood average gray level appearing in g (x, y); let cijThe probability of (i, j) is expressed as the joint probability
Wherein
Given an arbitrary threshold vector (s, t), regions A and D represent object and background, respectively, regions B and C represent edges and noise, leaving two classes C0And C1Respectively representing an object and a background; omegaiAnd muiRespectively represent class CiProbability and mean vector of;
C0and C1The probability of (d) can be expressed as:
C0and C1The vector of mean values of (a) may be expressed as:
the total average vector of the 2D histogram is:
the probability of leaving the diagonal is ignored, yielding: omega0+ω11, and μT≈ω0μ0+ω1μ1；
The inter-class dispersion matrix is defined as:
the traces of the discrete matrix are represented as:
thus, the optimal threshold is represented as:
the conditions are satisfied: s is less than or equal to f (x, y)*And g (x, y)>t*(ii) a And f (x, y)>s*And g (x, y) is less than or equal to t*Is ignored and set to 0 or 1;
thirdly, taking intersection of the threshold segmentation result and the processed image, and determining the light reflectivity of the image after the intersection is taken;
sequentially selecting function models from a preset function set, and calculating the dissolved oxygen and the permanganate by using the light reflectivity and the function models respectively to obtain the fitting degree of the dissolved oxygen and the fitting degree of the permanganate corresponding to each function model;
the reflectivity and the reference data are subjected to correlation statistics, the correlation degree between the inverted high-manganese acid salt index and the spectral reflectivity is searched, the second waveband reflectivity, the third waveband reflectivity, the product of the second waveband and the third waveband reflectivity, the sum of the second waveband and the third waveband reflectivity are respectively used as independent variables, the corresponding content of the permanganate index is used as a dependent variable, the functional forms such as a power function, an exponential function, a first-order polynomial and a second-order polynomial are used for fitting, and the selected independent variable is (rho)2+ρ3) The dependent variable is the permanganate content, and the function form is a second-order exponential form;
functional form: y is a.ebx+c·edx
Parameter values: a is 13.17; -8.742; c is 1.321; d-4.726
The degree of fit: SSE 0.1098; r-square 0.9596; adjusted R-square 0.9191;
RMSE＝0.1914
the reflectivity and the reference data are subjected to correlation statistics, the correlation degree between the dissolved oxygen content and the spectral reflectivity is searched, the first waveband reflectivity, the second waveband reflectivity, the sum of the first and second waveband reflectivities and the product of the first and second waveband reflectivities are respectively used as independent variables, the corresponding dissolved oxygen content is used as a dependent variable, and the fitting is carried out by using function forms such as a power function, an exponential function, a first-order polynomial and a second-order polynomial, the experimental result shows that the product of the first and second waveband reflectivities is used as the independent variable, the dissolved oxygen content is used as the dependent variable, and the fitting degree of the function type of a third-order polynomial is the best:
functional form: y ═ p1·x3+p2·x2+p3·x+p4
Parameter values: p is a radical of1＝1195000；p2＝-91810；p3＝2183；p4＝-8.265
The degree of fit: SSE 5.294; r-square 0.7923; adjusted R-square 0.5846;
RMSE＝1.328；
fifthly, respectively selecting a function model with the optimal fitting degree from the fitting degree of the dissolved oxygen and the fitting degree of the permanganate according to preset standards to fit;
fifthly, establishing four evaluation indexes, specifically:
the closer the corrected decision coefficient is to 1, the better the fitting degree is;
the closer the absolute value of the determinant coefficient is to 1, the better the fitting degree is;
the closer the sum of squared errors is to 0, the better the degree of fit;
the closer the root mean square error is to 0, the better the fitting degree is;
step two, calculating the absolute value, the error sum of squares and the root mean square error of the corrected decision coefficient, the absolute value, the error sum of squares and the root mean square error which correspond to each function model for the dissolved oxygen and the permanganate content;
step three, counting how many evaluation indexes of each function model are higher than those of other function models; and selecting the function model which meets the conditions and has the most evaluation indexes.
2. The method for fitting water body information based on remote sensing images according to claim 1, characterized by further comprising a multi-temporal remote sensing image registration step before the first step is executed, and specifically comprising:
a, performing Harris angular point feature extraction on a multi-temporal remote sensing image;
b, eliminating mismatching and mismatching vectors by using a normalized cross-correlation matching algorithm to obtain a registered image;
and C, processing the registered image as a remote sensing image in the first step.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710806824.1A CN107590816B (en) | 2017-09-08 | 2017-09-08 | Water body information fitting method based on remote sensing image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710806824.1A CN107590816B (en) | 2017-09-08 | 2017-09-08 | Water body information fitting method based on remote sensing image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107590816A CN107590816A (en) | 2018-01-16 |
| CN107590816B true CN107590816B (en) | 2021-06-15 |
Family
ID=61051032
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710806824.1A Active CN107590816B (en) | 2017-09-08 | 2017-09-08 | Water body information fitting method based on remote sensing image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107590816B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110472697B (en) * | 2019-08-22 | 2021-04-30 | 广州大学 | Water body identification method and device based on iterative classification |
| CN111426806A (en) * | 2020-03-27 | 2020-07-17 | 渤海大学 | Automatic monitoring and early warning method for freshness degree of aquatic product cold-chain logistics based on means of Internet of things |
| CN111949270A (en) * | 2020-07-27 | 2020-11-17 | 中国建设银行股份有限公司 | Method and device for sensing running environment change of process robot |
| CN113326855B (en) * | 2021-06-22 | 2022-07-12 | 长光卫星技术股份有限公司 | Night lamplight remote sensing image feature extraction method based on two-dimensional Gaussian surface fitting |
| CN117649612A (en) * | 2023-10-19 | 2024-03-05 | 成都大学 | Satellite hyperspectral remote sensing data surface water body extraction method based on hybrid algorithm |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101950364A (en) * | 2010-08-30 | 2011-01-19 | 西安电子科技大学 | Remote sensing image change detection method based on neighbourhood similarity and threshold segmentation |
| CN105243367A (en) * | 2015-10-12 | 2016-01-13 | 水利部水利信息中心 | Method and device for monitoring scope of water body based on satellite remote sensing data |
| CN106525762A (en) * | 2016-11-07 | 2017-03-22 | 航天恒星科技有限公司 | Water quality monitoring method and water quality monitoring device based on adaptive model |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101694282B1 (en) * | 2015-09-17 | 2017-01-09 | 계명대학교 산학협력단 | APPARATUS AND METHOD ABOUT PREDICTING CHLOROPHYLL-a FROM RIVER USING SATELLITE SENSOR DATA AND NONLINEAR RANSAC METHOD |
-
2017
- 2017-09-08 CN CN201710806824.1A patent/CN107590816B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101950364A (en) * | 2010-08-30 | 2011-01-19 | 西安电子科技大学 | Remote sensing image change detection method based on neighbourhood similarity and threshold segmentation |
| CN105243367A (en) * | 2015-10-12 | 2016-01-13 | 水利部水利信息中心 | Method and device for monitoring scope of water body based on satellite remote sensing data |
| CN106525762A (en) * | 2016-11-07 | 2017-03-22 | 航天恒星科技有限公司 | Water quality monitoring method and water quality monitoring device based on adaptive model |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107590816A (en) | 2018-01-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107590816B (en) | Water body information fitting method based on remote sensing image | |
| Guo et al. | FSDAF 2.0: Improving the performance of retrieving land cover changes and preserving spatial details | |
| CN112766075B (en) | Hyperspectral remote sensing black and odorous water body grading method based on semi-supervised learning strategy | |
| CN110765941A (en) | Seawater pollution area identification method and equipment based on high-resolution remote sensing image | |
| CN102750701B (en) | Method for detecting spissatus and spissatus shadow based on Landsat thematic mapper (TM) images and Landsat enhanced thematic mapper (ETM) images | |
| CN104867150A (en) | Wave band correction change detection method of remote sensing image fuzzy clustering and system thereof | |
| CN109490306B (en) | Pork freshness detection method based on color and smell data fusion | |
| CN110161233B (en) | Rapid quantitative detection method for immunochromatography test paper card | |
| CN112307901B (en) | SAR and optical image fusion method and system for landslide detection | |
| CN102096824A (en) | Multi-spectral image ship detection method based on selective visual attention mechanism | |
| CN109741322A (en) | A kind of visibility measurement method based on machine learning | |
| CN102360503B (en) | SAR (Specific Absorption Rate) image change detection method based on space approach degree and pixel similarity | |
| CN104102928A (en) | Remote sensing image classification method based on texton | |
| CN114279982A (en) | Water body pollution information acquisition method and device | |
| CN115909256B (en) | Road disease detection method based on road visual image | |
| CN102855485A (en) | Automatic wheat earing detection method | |
| CN105894520A (en) | Satellite image automatic cloud detection method based on Gaussian mixture model | |
| CN107610119A (en) | The accurate detection method of steel strip surface defect decomposed based on histogram | |
| CN102081799A (en) | Method for detecting change of SAR images based on neighborhood similarity and double-window filtering | |
| CN116559111A (en) | Sorghum variety identification method based on hyperspectral imaging technology | |
| CN106204596B (en) | Panchromatic waveband remote sensing image cloud detection method based on Gaussian fitting function and fuzzy mixed estimation | |
| CN111882573A (en) | Cultivated land plot extraction method and system based on high-resolution image data | |
| CN116895019A (en) | Remote sensing image change detection method and system based on dynamic weighted cross entropy loss | |
| Portalés et al. | An image-based system to preliminary assess the quality of grape harvest batches on arrival at the winery | |
| CN112991425B (en) | Water area water level extraction method and system and storage medium |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |