CN111898494A - Mining disturbed block boundary identification method - Google Patents

Mining disturbed block boundary identification method Download PDF

Info

Publication number
CN111898494A
CN111898494A CN202010683437.5A CN202010683437A CN111898494A CN 111898494 A CN111898494 A CN 111898494A CN 202010683437 A CN202010683437 A CN 202010683437A CN 111898494 A CN111898494 A CN 111898494A
Authority
CN
China
Prior art keywords
image
ndvi
formula
mining
area
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.)
Granted
Application number
CN202010683437.5A
Other languages
Chinese (zh)
Other versions
CN111898494B (en
Inventor
马占元
杨永均
田利军
朱笑笑
常晓华
侯湖平
王玉明
张绍良
许丽丽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Mining and Technology CUMT
Datong Coal Mine Group Co Ltd
Original Assignee
China University of Mining and Technology CUMT
Datong Coal Mine Group Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Mining and Technology CUMT, Datong Coal Mine Group Co Ltd filed Critical China University of Mining and Technology CUMT
Priority to CN202010683437.5A priority Critical patent/CN111898494B/en
Publication of CN111898494A publication Critical patent/CN111898494A/en
Application granted granted Critical
Publication of CN111898494B publication Critical patent/CN111898494B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Theoretical Computer Science (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Analysis (AREA)

Abstract

A mining disturbed block boundary identification method comprises the following steps of obtaining Landsat images of a research area, performing image preprocessing including radiation correction, atmospheric correction and image cutting, and calculating an NDVI index; iteration is carried out on the images in the same month according to the year sequence for difference making, and all the difference values and the absolute difference values are respectively accumulated to obtain an NDVI difference diagram and an NDVI absolute difference diagram of the research area; segmenting the two images by using an OTSU threshold segmentation method to obtain a threshold segmentation image, and performing post-processing of filtering and denoising on the threshold segmentation image to generate a mining area disturbance plot area and a damaged plot area; extracting boundaries of the disturbed land parcel and the damaged land parcel by using an edge detection algorithm; and performing geometric difference calculation on the disturbed land parcel boundary and the damaged land parcel boundary to obtain a recovered land parcel boundary. The method can accurately identify mining disturbed plots, damaged plots and restored plot boundaries, provide technical support for management of mining area ecological restoration, and provide important data for evaluation of ecological environment restoration quality.

Description

Mining disturbed block boundary identification method
Technical Field
The invention belongs to the field of mine ecological environment monitoring, and particularly relates to a mining disturbed block boundary identification method.
Background
The original landform and landform of the mining area land are changed along with the activities of digging, damage, occupation, collapse, pollution, ecological restoration and the like in the mining production process, and a mining disturbance land block is formed. Wherein the boundaries of the disturbed plot represent the range of influence of the mining activity, the disturbed plot including both types of damaged plots and restored plots. Activities such as assessment of ecological environment influence, law enforcement inspection, ecological restoration and the like in mining areas need to be effectively recognized for boundaries of mining disturbed plots.
In the prior art, the identification of the boundary of the disturbed land parcel is mainly carried out based on field acquisition data or single remote sensing image data. The practical application cost of collecting data or single remote sensing image data on the spot is too high and it is not effectively suitable for mining area. Meanwhile, the recognition accuracy is low, and reliable technical support cannot be provided for the restoration and management of the ecological environment.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a mining disturbed block boundary identification method, which can accurately and effectively identify the boundary of the mining disturbed block and provide reliable technical support for ecological environment restoration and management.
In order to achieve the purpose, the invention provides a mining disturbance plot boundary identification method, which adopts a long-time remote sensing image to identify a mining disturbance boundary and specifically comprises the following steps;
the method comprises the following steps: selecting a mining disturbance area as a research area, and acquiring a Landsat remote sensing image of the area;
a1: carrying out preprocessing of radiometric calibration, atmospheric correction, geometric correction and image cutting on the remote sensing image to obtain a preprocessed image of a research area;
a2: calculating an NDVI index map of each image according to a formula (1), and recording that the NDVI index map of the research area is NDVI;
Figure BDA0002586646680000011
in the formula, NIR is the reflection value of the near infrared band of each preprocessed image, and R is the reflection value of the red light band of each preprocessed image;
step two: carrying out NDVI normalization processing on all the NDVI index diagrams to obtain normalized NDVI index diagrams, and marking the normalized NDVI index diagrams as I;
b1: setting a regression linear equation between the reference image and the image to be normalized as shown in a formula (2);
Figure BDA0002586646680000021
in the formula, x is the NDVI value of the image to be normalized, and y is the NDVI value of the reference image;
b2: adopting a least square method to carry out straight line fitting between the reference image and the image to be normalized, and solving parameters according to a formula (3) and a formula (4) respectively
Figure BDA0002586646680000022
And
Figure BDA0002586646680000023
Figure BDA0002586646680000024
Figure BDA0002586646680000025
in the formula (I), the compound is shown in the specification,
Figure BDA0002586646680000026
and
Figure BDA0002586646680000027
respectively taking the mean values of the m pseudo feature points of the image to be normalized and the reference image, and respectively solving through a formula (5) and a formula (6);
Figure BDA0002586646680000028
Figure BDA0002586646680000029
step three: calculating the normalized NDVI index map to obtain an NDVI difference map and an NDVI absolute difference map of the research area;
c1: respectively carrying out iteration difference on the available images of each month according to the size sequence of the year through a formula (7) and a formula (8), obtaining an NDVI difference diagram and an NDVI absolute difference diagram of each month, and respectively marking as MdAnd AMd；
Figure BDA0002586646680000031
Figure BDA0002586646680000032
Where k denotes the number of available image years for the month, IiSorting the available images for the month by yearThe ith normalized processed NDVI index of (a);
c2: m of each month is expressed by formula (9) and formula (10), respectivelydAnd AMdCarrying out accumulation and summation to obtain an NDVI difference Image of the research areadAnd NDVI absolute difference map AImaged；
Figure BDA0002586646680000033
Figure BDA0002586646680000034
In the formula (I), the compound is shown in the specification,
Figure BDA0002586646680000035
and
Figure BDA0002586646680000036
an NDVI difference map and an NDVI absolute difference map representing an mth month, n representing the number of months of available pictures;
step four: obtaining an optimal segmentation threshold of the NDVI index differential image by using an OTSU threshold segmentation method, and performing threshold segmentation on the NDVI differential image of the research area;
d1: for the image with darker background area, the image size is x × y, the gray histogram is counted, and the proportion ω of the pixel points in the region of interest in the whole image is calculated through the formula (11) and the formula (12) respectively0The proportion omega of the pixel points in the background area in the whole image1；
ω0＝N0/(x×y) (11)；
ω1＝N1/(x×y) (12)；
In the formula, N0Number of pixels whose pixel value gradation is smaller than threshold value T, N1The number of pixels with the pixel value gray scale larger than the threshold value T,
and N is0And N1Satisfies formula (13); omega0And ω1Satisfies the formula (14);
N0+N1＝x×y (13)；
ω0+ω1＝1 (14)；
μ0the average gray scale of the pixel points in the interested area occupying the whole image is obtained;
μ1the average gray scale of the pixel points in the background area occupying the whole image;
g is an inter-class variance score;
μ is the total average gray scale of the image;
μ＝μ0×ω0+μ1×ω1(15)；
g＝ω0(μ0-μ)2+ω1(μ1-μ)2(16)；
d2: obtaining a formula (17) by combining a formula (15) and a formula (16), and solving the inter-class variance g;
g＝ω0×ω1×(μ0-μ1)2(17)；
d3: when the inter-class variance g reaches the maximum, the threshold value T is the optimal threshold value; calculating the NDVI difference diagram and the NDVI absolute difference diagram of the research region according to a formula (18) by taking T as a boundary to respectively obtain a threshold segmentation diagram corresponding to the NDVI difference diagram and the NDVI absolute difference diagram of the research region, and recording the threshold segmentation diagram as ImagetAnd AImaget；
Figure BDA0002586646680000041
D4: dividing graph Image into threshold valuestAnd AImagetImage denoising is carried out, and the Image after mining disturbance area processing is obtainednAnd AImagen；
Step five: automatically identifying boundaries of mining disturbed plots and damaged plots by adopting an edge detection method; detecting mining disturbed plots and damaged plot boundaries by using a Roberts edge detection operator through a formula (19) and a formula (20);
Figure BDA0002586646680000042
R(x,y)＝|f(x,y)-f(x+1,y+1)|+|f(x+1,y)-f(x,y+1)| (20)；
where f (x, y) is a pixel value at the coordinates of the input Image size (x, y), and R (x, y) is an output Image of the input Image size (x, y), and each is referred to as an edge detection ImageRAnd AImageR；
Step six: image of edge detectionRAnd AImageRRespectively converted into vector diagrams, and respectively recorded as vector diagram images of boundaries of mining disturbance plotsVAnd vector image AImage of damaged land block boundaryVAnd performing geometric difference calculation on the two images through a formula (21) to obtain a recovered block boundary vector ImageF；
ImageF＝AImageV-ImageV(21)。
Preferably, in the second step, NDVI normalization is performed by using a pseudo-invariant feature relative radiation normalization method, m pseudo-invariant feature points are selected, and NDVI normalization of the image is completed by using a regression equation.
Preferably, the threshold value is divided into the Image by a mean filtering methodtAnd AImagetAnd carrying out image denoising.
The method adopts the remote sensing image data to identify the boundary of the mining disturbed block, and has the advantages of low cost, high identification precision and wide coverage range. Because of the frequent disturbance of the land in the mining area, some plots may be quickly restored after being damaged. The traditional single-time-phase remote sensing image classification method and the double-time-phase remote sensing image change detection method cannot comprehensively reflect all disturbed regions, namely, regions which are quickly recovered after being damaged may be missed to be detected, the remote sensing images in summer are required to reflect surface changes, and when some regions lack summer images, the disturbed land parcel cannot be identified. The method for extracting the disturbed plaque by using the multi-temporal NDVI difference and the NDVI absolute difference has the advantages that on one hand, the historical accumulated disturbed boundary can be accurately extracted, on the other hand, the method has low requirement on the acquisition time of the remote sensing image, can use the available images in all seasons, and can be suitable for the areas with few available remote sensing images, so that the areas recovered after damage can be effectively identified, and all disturbed areas can be more comprehensively reflected. Therefore, the mining disturbance plot can be identified, and the mining disturbance plot can be further divided into a damaged plot and a recovery plot, so that a data basis can be provided for the ecological environment evaluation work of the mining area, and the mining disturbance plot has important significance for realizing the ecological environment recovery and management of the mining area. The method provided by the invention can accurately identify the boundary of the mining disturbed block, can provide a reliable foundation for mining planning of a mining area, and solves the problem of inaccurate identification of the disturbed block boundary caused by mining in the prior art.
Drawings
FIG. 1 is a schematic flow diagram of the present invention;
FIG. 2 is a NDVI difference diagram according to an embodiment of the invention;
FIG. 3 is an absolute difference diagram of NDVI according to an embodiment of the present invention;
FIG. 4 is a graph of NDVI differential threshold segmentation in accordance with an embodiment of the present invention;
FIG. 5 is a graph of NDVI absolute differential threshold segmentation in accordance with an embodiment of the present invention;
FIG. 6 is an NDVI differential filtering denoising diagram according to an embodiment of the present invention;
FIG. 7 is an NDVI absolute difference filter denoising diagram according to an embodiment of the present invention;
FIG. 8 is an NDVI differential edge detection diagram according to an embodiment of the invention;
FIG. 9 is an NDVI absolute differential edge detection diagram according to an embodiment of the invention;
FIG. 10 is a mining disturbance plot boundary vector diagram according to an embodiment of the present invention;
FIG. 11 is a damaged parcel boundary vector diagram according to an embodiment of the present invention;
fig. 12 is a recovered block boundary vector diagram according to an embodiment of the present invention.
Detailed Description
The invention will be further explained with reference to the drawings.
As shown in fig. 1, a mining disturbance plot boundary identification method adopts a long-time remote sensing image to identify a mining disturbance boundary, and specifically comprises the following steps;
the method comprises the following steps: selecting a mining disturbance area as a research area, and acquiring a Landsat remote sensing image of the area;
a1: carrying out preprocessing of radiometric calibration, atmospheric correction, geometric correction and image cutting on the remote sensing image to obtain a preprocessed image of a research area;
a2: calculating an NDVI index map of each image according to a formula (1), and recording that the NDVI index map of the research area is NDVI;
Figure BDA0002586646680000061
in the formula, NIR is the reflection value of the near infrared band of each preprocessed image, and R is the reflection value of the red light band of each preprocessed image;
step two: carrying out NDVI normalization processing on all the NDVI index diagrams to obtain normalized NDVI index diagrams, and marking the normalized NDVI index diagrams as I; as an optimization, NDVI normalization is carried out by adopting a pseudo-invariant feature relative radiation normalization method, m pseudo-invariant feature points are selected, and NDVI normalization of the image is completed by utilizing a regression equation; as an embodiment, a Landsat 8 remote sensing image of a certain research area may be used as a reference image, and 15 pseudo-invariant feature points are selected, where m is 15; the number of months of available images for this embodiment is 8;
b1: setting a regression linear equation between the reference image and the image to be normalized as shown in a formula (2);
Figure BDA0002586646680000071
in the formula, x is the NDVI value of the image to be normalized, and y is the NDVI value of the reference image;
b2: adopting a least square method to carry out straight line fitting between the reference image and the image to be normalized, and solving parameters according to a formula (3) and a formula (4) respectively
Figure BDA0002586646680000072
And
Figure BDA0002586646680000073
Figure BDA0002586646680000074
Figure BDA0002586646680000075
in the formula (I), the compound is shown in the specification,
Figure BDA0002586646680000076
and
Figure BDA0002586646680000077
respectively taking the mean values of the m pseudo feature points of the image to be normalized and the reference image, and respectively solving through a formula (5) and a formula (6);
Figure BDA0002586646680000078
Figure BDA0002586646680000079
when the m is 15, the reaction time is as short as possible,
Figure BDA00025866466800000710
step three: calculating the normalized NDVI index map to obtain an NDVI difference map and an NDVI absolute difference map of the research area;
c1: respectively carrying out iteration difference on the available images of each month according to the size sequence of the year through a formula (7) and a formula (8), obtaining an NDVI difference diagram and an NDVI absolute difference diagram of each month, and respectively marking as MdAnd AMd；
Figure BDA0002586646680000081
Figure BDA0002586646680000082
Where k denotes the number of available image years for the month, IiAn ith available normalized processed NDVI index sorted by year for the available images for the month;
c2: m of each month is expressed by formula (9) and formula (10), respectivelydAnd AMdCarrying out accumulation and summation to obtain an NDVI difference Image of the research areadAnd NDVI absolute difference map AImagedFor the embodiments in the present application, the corresponding images are shown in fig. 2 and fig. 3, respectively;
Figure BDA0002586646680000083
Figure BDA0002586646680000084
in the formula (I), the compound is shown in the specification,
Figure BDA0002586646680000085
and
Figure BDA0002586646680000086
NDVI difference map and NDVI absolute difference map representing the mth month, n representing the number of months of available pictures, for the embodiment in this application, n is 8;
when n is equal to 8, the reaction is carried out,
Figure BDA0002586646680000087
step four: obtaining an optimal segmentation threshold of the NDVI index differential image by using an OTSU threshold segmentation method, and performing threshold segmentation on the NDVI differential image of the research area;
d1: for the image with darker background area, the image size is x × y, the gray histogram is counted, and the proportion ω of the pixel points in the region of interest in the whole image is calculated through the formula (11) and the formula (12) respectively0And background region pixelsRatio omega of points to whole image1；
ω0＝N0/(x×y) (11)；
ω1＝N1/(x×y) (12)；
In the formula, N0Number of pixels whose pixel value gradation is smaller than threshold value T, N1The number of pixels whose pixel value gray scale is greater than threshold value T, and N0And N1Satisfies formula (13); omega0And ω1Satisfies the formula (14);
N0+N1＝x×y (13)；
ω0+ω1＝1 (14)；
μ0the average gray scale of the pixel points in the interested area occupying the whole image is obtained;
μ1the average gray scale of the pixel points in the background area occupying the whole image;
g is an inter-class variance score;
μ is the total average gray scale of the image;
μ＝μ0×ω0+μ1×ω1(15)；
g＝ω0(μ0-μ)2+ω1(μ1-μ)2(16)；
d2: obtaining a formula (17) by combining a formula (15) and a formula (16), and solving the inter-class variance g;
g＝ω0×ω1×(μ0-μ1)2(17)；
d3: when the inter-class variance g reaches the maximum, the threshold value T is the optimal threshold value; calculating the NDVI difference diagram and the NDVI absolute difference diagram of the research region according to a formula (18) by taking T as a boundary to respectively obtain a threshold segmentation diagram corresponding to the NDVI difference diagram and the NDVI absolute difference diagram of the research region, and recording the threshold segmentation diagram as ImagetAnd AImagetFor the embodiments in the present application, the corresponding structures are shown in fig. 4 and fig. 5, respectively;
Figure BDA0002586646680000091
d4: dividing graph Image into threshold valuestAnd AImagetImage denoising is carried out, and the Image after mining disturbance area processing is obtainednAnd AImagen(ii) a Preferably, the threshold value is divided into the Image by a mean filtering methodtAnd AImagetCarrying out image denoising; for the examples in this application, the corresponding results are shown in fig. 6 and fig. 7, respectively;
step five: automatically identifying boundaries of mining disturbed plots and damaged plots by adopting an edge detection method; detecting mining disturbed plots and damaged plot boundaries by using a Roberts edge detection operator through a formula (19) and a formula (20);
Figure BDA0002586646680000101
R(x,y)＝|f(x,y)-f(x+1,y+1)|+|f(x+1,y)-f(x,y+1)| (20)；
where f (x, y) is a pixel value at the coordinates of the input Image size (x, y), and R (x, y) is an output Image of the input Image size (x, y), and each is referred to as an edge detection ImageRAnd AImageR；
The results are shown in fig. 8 and fig. 9, respectively, for the examples in the present application.
Step six: image of edge detectionRAnd AImageRRespectively converted into vector diagrams, and respectively recorded as vector diagram images of boundaries of mining disturbance plotsVAnd vector image AImage of damaged land block boundaryVThe results are shown in fig. 10 and fig. 11 for the examples in the present application, respectively. The two are geometrically differenced through a formula (21) to obtain a recovered plot boundary vector diagram ImageF(ii) a The results are shown in fig. 12 for the examples in this application.
ImageF＝AImageV-ImageV(21)。

Claims (3)

1. A mining disturbance plot boundary identification method adopts a long-time remote sensing image to carry out mining disturbance boundary identification, and is characterized by comprising the following steps;
the method comprises the following steps: selecting a mining disturbance area as a research area, and acquiring a Landsat remote sensing image of the area;
a1: carrying out preprocessing of radiometric calibration, atmospheric correction, geometric correction and image cutting on the remote sensing image to obtain a preprocessed image of a research area;
a2: calculating an NDVI index map of each image according to a formula (1), and recording that the NDVI index map of the research area is NDVI;
Figure FDA0002586646670000011
in the formula, NIR is the reflection value of the near infrared band of each preprocessed image, and R is the reflection value of the red light band of each preprocessed image;
step two: carrying out NDVI normalization processing on all the NDVI index diagrams to obtain normalized NDVI index diagrams, and marking the normalized NDVI index diagrams as I;
b1: setting a regression linear equation between the reference image and the image to be normalized as shown in a formula (2);
Figure FDA0002586646670000012
in the formula, x is the NDVI value of the image to be normalized, and y is the NDVI value of the reference image;
b2: adopting a least square method to carry out straight line fitting between the reference image and the image to be normalized, and solving parameters according to a formula (3) and a formula (4) respectively
Figure FDA0002586646670000013
And
Figure FDA0002586646670000014
Figure FDA0002586646670000016
in the formula (I), the compound is shown in the specification,
Figure FDA0002586646670000017
and
Figure FDA0002586646670000018
respectively taking the mean values of the m pseudo feature points of the image to be normalized and the reference image, and respectively solving through a formula (5) and a formula (6);
Figure FDA0002586646670000021
Figure FDA0002586646670000022
step three: calculating the normalized NDVI index map to obtain an NDVI difference map and an NDVI absolute difference map of the research area;
c1: respectively carrying out iteration difference on the available images of each month according to the size sequence of the year through a formula (7) and a formula (8), obtaining an NDVI difference diagram and an NDVI absolute difference diagram of each month, and respectively marking as MdAnd AMd；
Figure FDA0002586646670000023
Figure FDA0002586646670000024
Where k denotes the number of available image years for the month, IiAn ith available normalized processed NDVI index sorted by year for the available images for the month;
c2: respectively by the formula (9) And equation (10) will be M for each monthdAnd AMdCarrying out accumulation and summation to obtain an NDVI difference Image of the research areadAnd NDVI absolute difference map AImaged；
Figure FDA0002586646670000025
Figure FDA0002586646670000026
In the formula (I), the compound is shown in the specification,
Figure FDA0002586646670000027
and
Figure FDA0002586646670000028
an NDVI difference map and an NDVI absolute difference map representing an mth month, n representing the number of months of available pictures;
step four: obtaining an optimal segmentation threshold of the NDVI index differential image by using an OTSU threshold segmentation method, and performing threshold segmentation on the NDVI differential image of the research area;
d1: for the image with darker background area, the image size is x × y, the gray histogram is counted, and the proportion ω of the pixel points in the region of interest in the whole image is calculated through the formula (11) and the formula (12) respectively0The proportion omega of the pixel points in the background area in the whole image1；
ω0＝N0/(x×y) (11)；
ω1＝N1/(x×y) (12)；
In the formula, N0Number of pixels whose pixel value gradation is smaller than threshold value T, N1The number of pixels whose pixel value gray scale is greater than threshold value T, and N0And N1Satisfies formula (13); omega0And ω1Satisfies the formula (14);
N0+N1＝x×y (13)；
ω0+ω1＝1 (14)；
μ0the average gray scale of the pixel points in the interested area occupying the whole image is obtained;
μ1the average gray scale of the pixel points in the background area occupying the whole image;
g is an inter-class variance score;
μ is the total average gray scale of the image;
μ＝μ0×ω0+μ1×ω1(15)；
g＝ω0(μ0-μ)2+ω1(μ1-μ)2(16)；
d2: obtaining a formula (17) by combining a formula (15) and a formula (16), and solving the inter-class variance g;
g＝ω0×ω1×(μ0-μ1)2(17)；
d3: when the inter-class variance g reaches the maximum, the threshold value T is the optimal threshold value; calculating the NDVI difference diagram and the NDVI absolute difference diagram of the research region according to a formula (18) by taking T as a boundary to respectively obtain a threshold segmentation diagram corresponding to the NDVI difference diagram and the NDVI absolute difference diagram of the research region, and recording the threshold segmentation diagram as ImagetAnd AImaget；
Figure FDA0002586646670000041
D4: dividing graph Image into threshold valuestAnd AImagetImage denoising is carried out, and the Image after mining disturbance area processing is obtainednAnd AImagen；
Step five: automatically identifying boundaries of mining disturbed plots and damaged plots by adopting an edge detection method; detecting mining disturbed plots and damaged plot boundaries by using a Roberts edge detection operator through a formula (19) and a formula (20);
Figure FDA0002586646670000042
R(x,y)＝|f(x,y)-f(x+1,y+1)|+|f(x+1,y)-f(x,y+1)| (20)；
where f (x, y) is a pixel value at the coordinates of the input Image size (x, y), and R (x, y) is an output Image of the input Image size (x, y), and each is referred to as an edge detection ImageRAnd AImageR；
Step six: image of edge detectionRAnd AImageRRespectively converted into vector diagrams, and respectively recorded as vector diagram images of boundaries of mining disturbance plotsVAnd vector image AImage of damaged land block boundaryVAnd performing geometric difference calculation on the two images through a formula (21) to obtain a recovered block boundary vector ImageF；
ImageF＝AImageV-ImageV(21)。
2. The mining disturbed block boundary identification method according to claim 1, characterized in that in the second step, NDVI normalization is performed by adopting a pseudo-invariant feature relative radiation normalization method, m pseudo-invariant feature points are selected, and NDVI normalization of the image is completed by using a regression equation.
3. A mining disturbance block boundary identification method as claimed in claim 1 or 2, characterized in that in step four, in D4, a mean value filtering method is adopted to threshold segmentation ImagetAnd AImagetAnd carrying out image denoising.
CN202010683437.5A 2020-07-16 2020-07-16 Mining disturbed land boundary identification method Active CN111898494B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010683437.5A CN111898494B (en) 2020-07-16 2020-07-16 Mining disturbed land boundary identification method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010683437.5A CN111898494B (en) 2020-07-16 2020-07-16 Mining disturbed land boundary identification method

Publications (2)

Publication Number Publication Date
CN111898494A true CN111898494A (en) 2020-11-06
CN111898494B CN111898494B (en) 2022-09-27

Family

ID=73192055

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010683437.5A Active CN111898494B (en) 2020-07-16 2020-07-16 Mining disturbed land boundary identification method

Country Status (1)

Country Link
CN (1) CN111898494B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112906648A (en) * 2021-03-24 2021-06-04 深圳前海微众银行股份有限公司 Method and device for classifying objects in land parcel and electronic equipment
CN113486809A (en) * 2021-07-08 2021-10-08 中国矿业大学 Mining area influence boundary identification method
CN113592877A (en) * 2021-03-25 2021-11-02 国网新源控股有限公司 Method and device for identifying red line exceeding of pumped storage power station

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090034868A1 (en) * 2007-07-30 2009-02-05 Rempel Allan G Enhancing dynamic ranges of images
CN103971115A (en) * 2014-05-09 2014-08-06 中国科学院遥感与数字地球研究所 Automatic extraction method for newly-increased construction land image spots in high-resolution remote sensing images based on NDVI and PanTex index
CN104268879A (en) * 2014-09-28 2015-01-07 民政部国家减灾中心 Physical building quantity damage evaluation method based on remote sensing multi-spectral images
CN107103584A (en) * 2017-04-11 2017-08-29 北京师范大学 A kind of production high-spatial and temporal resolution NDVI weighted based on space-time method
CN108592888A (en) * 2018-04-23 2018-09-28 中国科学院地球化学研究所 A kind of Residential area extraction method
CN109598273A (en) * 2018-12-03 2019-04-09 中国矿业大学 A kind of city entity boundary recognition methods of fusion surface temperature and building index
CN110598553A (en) * 2019-08-09 2019-12-20 中国科学院南京地理与湖泊研究所 Original true landform mining damaged area detection method based on remote sensing image and topographic data
CN111199195A (en) * 2019-12-26 2020-05-26 中科禾信遥感科技（苏州）有限公司 Pond state full-automatic monitoring method and device based on remote sensing image

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090034868A1 (en) * 2007-07-30 2009-02-05 Rempel Allan G Enhancing dynamic ranges of images
CN103971115A (en) * 2014-05-09 2014-08-06 中国科学院遥感与数字地球研究所 Automatic extraction method for newly-increased construction land image spots in high-resolution remote sensing images based on NDVI and PanTex index
CN104268879A (en) * 2014-09-28 2015-01-07 民政部国家减灾中心 Physical building quantity damage evaluation method based on remote sensing multi-spectral images
CN107103584A (en) * 2017-04-11 2017-08-29 北京师范大学 A kind of production high-spatial and temporal resolution NDVI weighted based on space-time method
CN108592888A (en) * 2018-04-23 2018-09-28 中国科学院地球化学研究所 A kind of Residential area extraction method
CN109598273A (en) * 2018-12-03 2019-04-09 中国矿业大学 A kind of city entity boundary recognition methods of fusion surface temperature and building index
CN110598553A (en) * 2019-08-09 2019-12-20 中国科学院南京地理与湖泊研究所 Original true landform mining damaged area detection method based on remote sensing image and topographic data
CN111199195A (en) * 2019-12-26 2020-05-26 中科禾信遥感科技（苏州）有限公司 Pond state full-automatic monitoring method and device based on remote sensing image

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
HONG JIANG: "Vegetation monitoring in rugged terrain with one novel topography — Adjusted vegetation index (TAVI)", 《2010 3RD INTERNATIONAL CONGRESS ON IMAGE AND SIGNAL PROCESSING》 *
ZHEN YANG: "Identification of the disturbance and trajectory types in mining areas using multitemporal remote sensing image", 《SCIENCE OF THE TOTAL ENVIRONMENT》 *
周新瑞: "面向自然保护区周边的人为扰动影响遥感探测分析——以永德大雪山国家级自然保护区为例", 《林业资源管理》 *
张晋纶: "黄土高原煤矿区地表采动裂缝扰动范围预计方法研究", 《中国煤炭》 *
李晶等: "基于时序多光谱影像的干旱草原区开采扰动信息提取方法", 《光谱学与光谱分析》 *
杨永均: "矿山土地生态***恢复力及其测度与调控研究", 《中国博士学位论文全文数据库 (工程科技Ⅰ辑)》 *
毕永清等: "基于纹理的高寒地区人为扰动地表信息提取", 《测绘与空间地理信息》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112906648A (en) * 2021-03-24 2021-06-04 深圳前海微众银行股份有限公司 Method and device for classifying objects in land parcel and electronic equipment
CN113592877A (en) * 2021-03-25 2021-11-02 国网新源控股有限公司 Method and device for identifying red line exceeding of pumped storage power station
CN113592877B (en) * 2021-03-25 2024-04-12 国网新源控股有限公司 Method and device for identifying red line exceeding of pumped storage power station
CN113486809A (en) * 2021-07-08 2021-10-08 中国矿业大学 Mining area influence boundary identification method

Also Published As

Publication number Publication date
CN111898494B (en) 2022-09-27

Similar Documents

Publication Publication Date Title
CN111898494B (en) Mining disturbed land boundary identification method
Malambo et al. Automated detection and measurement of individual sorghum panicles using density-based clustering of terrestrial lidar data
CN105784713B (en) Sealing ring detection method of surface flaw based on machine vision
CN106910187B (en) Image data set artificial amplification method for bridge crack detection and positioning
CN108921813B (en) Unmanned aerial vehicle detection bridge structure crack identification method based on machine vision
CN105389586A (en) Method for automatically detecting integrity of shrimp body based on computer vision
CN105719283A (en) Road surface crack image detection method based on Hessian matrix multi-scale filtering
CN111007013B (en) Crop rotation fallow remote sensing monitoring method and device for northeast cold region
CN115115627B (en) Soil saline-alkali soil monitoring method based on data processing
CN105825169A (en) Road-image-based pavement crack identification method
CN111047566A (en) Method for carrying out aquatic vegetation annual change statistics by unmanned aerial vehicle and multispectral satellite image
CN115690081A (en) Tree counting method, system, storage medium, computer equipment and terminal
CN113870202A (en) Far-end chip defect detection system based on deep learning technology
CN107516315B (en) Tunneling machine slag tapping monitoring method based on machine vision
Tao et al. Recognition of ecological vegetation fairy circles in intertidal salt marshes from UAV LiDAR point clouds
CN113076959A (en) Concrete structure surface microcrack feature extraction method based on convolutional neural network
CN116310826B (en) High-resolution remote sensing image forest land secondary classification method based on graphic neural network
CN116363655A (en) Financial bill identification method and system
CN110348443B (en) Multi-feature separation statistical method for tree trunk in complex environment and trunk identification method
David et al. Pathway detection and geometrical description from ALS data in forested mountaneous area
Yao et al. An adaptive and quantitative rubber nutrient status analyzing system by digital foliar images
Foschi A geometric approach to a mixed pixel problem: Detecting subpixel woody vegetation
CN113837996B (en) Automatic subway tunnel defect detection method supporting manual verification
McCarthy et al. Automated internode length measurement of cotton plants under field conditions
CN114049490B (en) Intelligent occupational health early warning method and system

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