CN117368920B - 基于D-insar的采煤区沉陷监测方法及系统 - Google Patents
基于D-insar的采煤区沉陷监测方法及系统 Download PDFInfo
- Publication number
- CN117368920B CN117368920B CN202311657974.2A CN202311657974A CN117368920B CN 117368920 B CN117368920 B CN 117368920B CN 202311657974 A CN202311657974 A CN 202311657974A CN 117368920 B CN117368920 B CN 117368920B
- Authority
- CN
- China
- Prior art keywords
- obtaining
- value
- fluctuation
- degree
- calculating
- 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
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000012544 monitoring process Methods 0.000 title claims abstract description 35
- 238000005065 mining Methods 0.000 title claims abstract description 30
- 239000003245 coal Substances 0.000 title claims abstract description 29
- 238000004891 communication Methods 0.000 claims abstract description 95
- 230000002159 abnormal effect Effects 0.000 claims abstract description 51
- 238000009499 grossing Methods 0.000 claims abstract description 36
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 26
- 230000005856 abnormality Effects 0.000 claims abstract description 22
- 238000001914 filtration Methods 0.000 claims abstract description 14
- 238000012512 characterization method Methods 0.000 claims description 16
- 238000010586 diagram Methods 0.000 claims description 15
- 230000000694 effects Effects 0.000 claims description 14
- 238000013507 mapping Methods 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000005516 engineering process Methods 0.000 description 5
- 238000006073 displacement reaction Methods 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000010606 normalization Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Processing (AREA)
Abstract
本发明涉及摄影测量技术领域,具体涉及一种基于D‑insar的采煤区沉陷监测方法及系统;根据干涉图中的灰度差异特征获得灰度频率分布图和波动值区间;对干涉图分解获得不同的信息图像;根据波动值区间获得信息图像不同的连通范围;根据连通范围之间的距离特征获得噪声程度;根据连通范围中的灰度差异特征获得像素点的异常指数;根据异常指数和噪声程度获得像素点的异常特征程度;根据所有信息图像不同像素点的异常特征程度和噪声程度,获得连通范围的综合特征程度和平滑系数。本发明根据自适应的平滑系数通过非局部均值滤波算法去噪,提高了去噪和监测地质沉陷的准确性。
Description
技术领域
本发明涉及摄影测量技术领域,具体涉及一种基于D-insar的采煤区沉陷监测方法及系统。
背景技术
在煤矿开采过程中,由于地下空间结构发生改变,可能会引发地表的沉陷,为了保证地质环境安全,需要对地质沉陷情况进行监测。D-insar差分干涉合成孔径雷达技术是利用两个或多个拍摄时间、观测角度和极化状态相同的SAR合成孔径雷达图像,获取地表高程和位移信息的一种雷达遥感技术。
该技术常用于监测地质表面动态变化,在应用过程中需要将多张SAR图像进行处理生成干涉图,位移场的生成效果取决于干涉图的去噪效果的好坏。非局部均值滤波算法是现有一种常用的去噪方法,广泛应用于干涉图的去噪;但由于干涉图中的噪声分布极不均匀,使得难以确定该去噪算法的平滑参数,对去噪效果的影响较大;最终导致监测地址沉陷的准确性降低。
发明内容
为了解决上述非局部均值滤波算法对干涉图的去噪效果较差,导致监测地址沉陷的准确性降低的技术问题,本发明的目的在于提供一种基于D-insar的采煤区沉陷监测方法及系统,所采用的技术方案具体如下:
获取采煤区的干涉图;根据所述干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图;根据所述灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间;
对所述干涉图通过EMD分解算法获得不同的信息图像;根据所述波动值区间获得所述信息图像的不同的连通范围;根据信息图像中所述连通范围之间的距离特征获得连通范围的噪声程度;
根据所述连通范围中的灰度差异特征获得像素点的异常指数;根据所述异常指数和所述噪声程度获得像素点的异常特征程度;根据所有信息图像中相同的连通范围内的不同像素点的异常特征程度的差异特征和噪声程度,获得所述连通范围的综合特征程度;
根据所述综合特征程度获得连通范围的平滑系数;根据所述平滑系数通过非局部均值滤波算法对干涉图去噪获得去噪干涉图,根据去噪干涉图监测地质沉陷。
进一步地,所述根据所述干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图的步骤包括:
将干涉图中任意像素点与预设范围内任意其他像素点构建二元组,并计算所述二元组的灰度差值绝对值,获得任意像素点的邻域灰度差;根据干涉图中所有像素点的邻域灰度差构建横坐标为邻域灰度差、纵坐标为邻域灰度差的频率的直方图,获得干涉图的灰度频率分布图。
进一步地,所述根据所述灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间的步骤包括:
计算所述灰度频率分布图中任意邻域灰度差对应的所述二元组的种类数与所有二元组种类数的比值,获得二元组种类占比,计算所述任意邻域灰度差对应的频率与所述二元组种类占比的乘积并归一化,获得波动权重,计算所述波动权重与所述任意邻域灰度差的乘积,获得波动范围值;
将所述灰度频率分布图中的所述邻域灰度差从小到大排序,计算从正序开始的预设第一数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最小值;计算从倒序开始的预设第二数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最大值;根据所述波动最小值和所述波动最大值构建区间,获得所述干涉图的波动值区间。
进一步地,所述根据所述波动值区间获得所述信息图像的不同的连通范围的步骤包括:
对于任意信息图像中的任意像素点,当任意像素点与预设步长内其他像素点的灰度差值绝对值小于所述波动值区间的任意波动值,则将任意像素点与预设步长内其他像素点构建连通域,遍历所述任意信息图像中所有像素点,获得不同的连通域;当任意两个连通域的平均灰度值的差值小于预设差异值,认为所述任意两个连通域属于同类连通域,获得所述任意信息图像在任意波动值下的不同类别的连通域;
对比所有信息图像中根据不同波动值构建的连通域数量,将连通域数量最多的不同连通域作为所有信息图像中不同的连通范围。
进一步地,所述根据信息图像中所述连通范围之间的距离特征获得连通范围的噪声程度的步骤包括:
对于所述任意信息图像中任意波动值下的任意类别连通域,计算所述任意类别与其他类别的任意两个连通域的中心像素点之间的欧氏距离,获得连通距离差异,计算连通距离差异的最大类间方差值并负相关映射,获得任意类别连通域的规律度;
计算所述任意信息图像中所有包括同一个所述连通范围的任意类别连通域的规律度的平均值并负相关映射,获得所述任意信息图像中所述连通范围的噪声程度。
进一步地,所述根据所述连通范围中的灰度差异特征获得像素点的异常指数的步骤包括:
计算信息图像中任意像素点的预设邻域范围与任意像素点所在连通范围之间的灰度值方差的比值,作为孤立森林算法中的数据特征,对任意像素点所在连通范围通过孤立森林算法计算任意像素点的异常分数值,获得信息图像中任意像素点的异常指数。
进一步地,所述根据所述异常指数和所述噪声程度获得像素点的异常特征程度的步骤包括:
对于任意信息图像的任意连通范围中的任意像素点,计算所述任意像素点的异常指数的归一化值与对应的连通范围的所述噪声程度的乘积,获得异常特征程度。
进一步地,所述获得所述连通范围的综合特征程度的步骤包括:
将所有信息图像的同一位置的像素点的所述异常特征程度构建序列,获得异常表征序列,计算同一连通范围任意两个异常表征序列的皮尔相关系数绝对值的平均值并归一化,获得特征波动程度;计算所有信息图像的同一连通范围的所述噪声程度的平均值并负相关映射,获得区域稳定表征值;计算所述特征波动程度与所述区域稳定表征值的乘积,获得所述连通范围的综合特征程度。
进一步地,所述根据所述综合特征程度获得连通范围的平滑系数的步骤包括:
将所述综合特征程度负相关映射,获得平滑权重,计算预设平滑值与平滑权重的乘积,获得所述连通范围对应的干涉图的平滑系数。
本发明还提出了一种基于D-insar的采煤区沉陷监测系统,包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序,所述处理器执行所述计算机程序实现任意一项一种基于D-insar的采煤区沉陷监测方法的步骤。
本发明具有如下有益效果:
在本发明实施例中,获取灰度频率分布图能够反映干涉图中相邻像素点之间的灰度差异情况的分布特征,根据灰度频率分布图确定波动值区间能够确定干涉图在不同波动值下的不同的连通范围,便于确定不同连通范围的噪声特征。获取干涉图不同的信息图像能够基于不同频率下的噪声特征提高最终去噪效果的准确性;根据信息图像中连通范围之间的距离特征获得连通范围的噪声程度,能够确定不同连通范围的噪声特征,进而为自适应的平滑系数提供基础。获取噪声点的异常指数和异常特征程度能够更准确地分析连通范围内噪声点的分布情况以及对异常情况进行具体化;获得综合特征程度能够准确地基于所有信息图像的异常特征程度和噪声程度确定该连通范围的平滑系数。最终根据不同连通范围的不同平滑系数通过非局部均值滤波进行去噪,提高去噪效果,保证了监测地质沉陷的准确性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案和优点,下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其它附图。
图1为本发明一个实施例所提供的一种基于D-insar的采煤区沉陷监测方法流程图。
具体实施方式
为了更进一步阐述本发明为达成预定发明目的所采取的技术手段及功效,以下结合附图及较佳实施例,对依据本发明提出的一种基于D-insar的采煤区沉陷监测方法及系统,其具体实施方式、结构、特征及其功效,详细说明如下。在下述说明中,不同的“一个实施例”或“另一个实施例”指的不一定是同一实施例。此外,一或多个实施例中的特定特征、结构或特点可由任何合适形式组合。
除非另有定义,本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。
下面结合附图具体的说明本发明所提供的一种基于D-insar的采煤区沉陷监测方法及系统的具体方案。
请参阅图1,其示出了本发明一个实施例提供一种基于D-insar的采煤区沉陷监测方法流程图,该方法包括以下步骤:
步骤S1,获取采煤区的干涉图;根据干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图;根据灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间。
在本发明实施例中,实施场景为对监测地质沉陷所需要的干涉图进行去噪。首先获取采煤区的干涉图,通过D-insar技术获取至少两次的SAR观测图像数据,且需要覆盖相同的采煤地区,通过共轨焦处理将SAR图像对齐得到复数干涉图,需要说明的是,通过D-insar技术获取干涉图属于现有技术,具体步骤不再赘述。
由于采集的干涉图非常容易受到噪声干扰,在后续的相位解缠、影像配准等处理过程中会出现较大的误差,非局部均值滤波算法是干涉图常用的现有去噪方法,该算法中的平滑参数的设定对最终的去噪效果有较大的影响,但干涉图中的噪声分布不均匀,且规模较大,无法很好的确定合适的平滑参数进行滤波处理,因此需要结合需要去噪的干涉图的噪声特征进行自适应选取平滑参数;因非局部均值滤波算法属于现有技术,具体去噪步骤不再赘述。
进一步地,对干涉图去噪则需要分析干涉图中的噪声分布特征,当像素点与邻域内的其他像素点的灰度差异特征越大,意味着噪声越严重;故可根据干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图。
优选地,在本发明一个实施例中,获取灰度频率分布图包括:将干涉图中任意像素点与预设范围内任意其他像素点构建二元组,在本发明实施例中预设范围为像素点的八邻域,实施者可根据实施场景自行确定;计算二元组的灰度差值绝对值,获得任意像素点的邻域灰度差;根据干涉图中所有像素点的邻域灰度差构建横坐标为邻域灰度差、纵坐标为邻域灰度差的频率的直方图,获得干涉图的灰度频率分布图,当二元组内两个像素点相同时只统计一次邻域灰度差。当灰度频率分布图中邻域灰度差的种类越多,二元组的种类越多,意味着该干涉图的相邻像素点之间的灰度差异特征的种类越多,噪声情况越明显。
获得干涉图的灰度频率分布图后,需要分析干涉图中不同区域的噪声特征,故可根据灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间;优选地,获取波动值区间包括:计算灰度频率分布图中任意邻域灰度差对应的二元组的种类数与所有二元组种类数的比值,获得二元组种类占比,横轴上任意一个邻域灰度差都可能是由不同二元组的灰度差值绝对值得到的。计算任意邻域灰度差对应的频率与二元组种类占比的乘积并归一化,获得波动权重,当该二元组种类越多、频率越高,意味着该邻域灰度差在该干涉图中是越明显的,邻域灰度差表征的是相邻像素点之间的波动,该波动程度在干涉图中是越明显的,则计算波动范围中占比权重越高;计算波动权重与对应的任意邻域灰度差的乘积,获得波动范围值;
将灰度频率分布图中的邻域灰度差从小到大排序,计算从正序开始的预设第一数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最小值;在本发明实施例中预设第一数量为5,即前五个最小的邻域灰度差对应的波动范围平均值获得波动最小值,该波动范围最小值表征了干涉图中噪声程度最小的灰度差异的特征。计算从倒序开始的预设第二数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最大值;在本发明实施例中预设第二数量为5,即前五个最大的邻域灰度差对应的波动范围平均值获得波动最大值,该波动范围最大值表征了干涉图中噪声程度最明显的灰度差异的特征。根据波动最小值和波动最大值构建区间,获得干涉图的波动值区间,该波动值区间表征了干涉图中不同区域噪声特征的不同程度。获取波动范围最小值的公式包括:
式中,表示波动范围最小值,/>表示预设第一数量,/>表示第/>个邻域灰度差对应的频率,/>表示第/>个邻域灰度差对应的二元组种类占比,/>表示归一化函数,表示波动权重,/>表示第/>个邻域灰度差。
步骤S2,对干涉图通过EMD分解算法获得不同的信息图像;根据波动值区间获得信息图像的不同的连通范围;根据信息图像中连通范围之间的距离特征获得连通范围的噪声程度。
获得干涉图的波动值区间后,可根据波动值区间确定不同区域的噪声程度,为了提高确定不同区域的噪声程度的准确性,首先将干涉图通过EMD分解算法获得不同的信息图像,需要说明的是EMD分解算法属于现有技术,具体分解步骤不再赘述;不同的信息图像包含着不同频率的信息,因为噪声对应的是高频的信息,故频率越高的信息图像,其噪声特征越明显,后续可根据不同信息图像的相同区域的噪声特征确定非局部均值滤波算法的平滑系数。
获得信息图像后,因为不同区域的噪声特征差异较为明显且地质沉陷区域的噪声更为严重,因此需要对不同噪声程度的区域进行区分,对不同区域进行不同程度的去噪,提高去噪效果;故根据波动值区间获得信息图像的不同的连通范围。
优选地,在本发明一个实施例中,获取连通范围包括:对于任意信息图像中的任意像素点,当任意像素点与预设步长内其他像素点的灰度差值绝对值小于波动值区间的任意波动值,则将任意像素点与预设步长内其他像素点构建连通域,在本发明实施例中预设步长为2,实施者可根据实施场景自行确定,例如任意波动值为5,则将任意像素点与步长2内的其他像素点灰度差小于5的构建为连通域;遍历该任意信息图像中所有像素点,获得不同的连通域;若任意一个像素点不能够构建连通域,则剔除该像素点。当任意两个连通域的平均灰度值的差值小于预设差异值,在本发明实施例中预设差异值为3,认为任意两个连通域属于同类连通域,进而获得该任意信息图像在波动值区间中每个波动值下所划分的不同类别的连通域。
对比所有信息图像中根据不同波动值构建的连通域数量,当波动值越小,划分的连通域的数量越多,越有利于提高去噪效果,将连通域数量最多的不同连通域作为所有信息图像中不同的连通范围,该连通范围在每个信息图像中的位置相同。
进一步地,当干涉图中相同的地理类型在不包含噪声的情况下,根据干涉图的特性,相近的区域其灰度特征相似,连通范围的分布具有一定的规律性,但随着地质沉陷区域以及噪声的出现,其连通范围的分布规律性被破坏,导致不同类别的连通范围的分布随机,因此可根据信息图像中连通范围之间的距离特征获得连通范围的噪声程度。
优选地,在本发明一个实施例中,获取噪声程度包括:对于任意信息图像中任意波动值下的任意类别连通域,计算该任意类别与其他类别的任意两个连通域的中心像素点之间的欧氏距离,获得连通距离差异;例如每个类别连通域有着多个连通域区域,计算该任意类别的任意一个连通域与其他任意类别的任意一个连通域区域的中心像素点的欧氏距离,得到不同的连通距离差异;计算连通距离差异的最大类间方差值并负相关映射,获得任意类别连通域的规律度;最大类间方差法属于现有技术,具体计算步骤不再赘述,当最大类间方差值越大,则连通距离差异的大小越离散,意味着该任意类别连通域的每个连通域与其他类别连通域之间的分布越不规律,与其他类别之间的连通域分布差异越大,该任意类别的连通域的分布的规律度越小,噪声特征越明显。
计算该任意信息图像中所有包括同一个连通范围的任意类别连通域的规律度的平均值并负相关映射,获得该任意信息图像中所述连通范围的噪声程度;计算在不同波动值下获得的包括同一个连通范围的任意类别连通域的规律度的平均值是为了提高噪声程度的准确性,当在所有波动值下对应的该任意类别连通域的规律度都越小,则意味着噪声越明显。
步骤S3,根据连通范围中的灰度差异特征获得像素点的异常指数;根据异常指数和噪声程度获得像素点的异常特征程度;根据所有信息图像中相同的连通范围内的不同像素点的异常特征程度的差异特征和噪声程度,获得连通范围的综合特征程度。
获得不同连通范围的噪声程度后,则需要具体分析该连通范围内不同像素点的异常情况,当该连通范围内异常的像素点越多,也意味着噪声越明显,故根据连通范围中的灰度差异特征获得像素点的异常指数;优选地,获取异常指数包括:计算信息图像中任意像素点的预设邻域范围与任意像素点所在连通范围之间的灰度值方差的比值,在本发明实施例中预设邻域范围为像素点的八邻域,该比值能够反映该任意像素点的八邻域内的灰度差异特征与该区域的灰度差异特征的差异,当比值越大,意味着该任意像素点处的灰度差异越明显,该位置为噪声的可能性越大,故将该比值作为孤立森林算法中的数据特征,对任意像素点所在连通范围通过孤立森林算法计算任意像素点的异常分数值,获得信息图像中任意像素点的异常指数。需要说明的是,孤立森林算法属于现有技术,用来寻找异常数据点,具体计算步骤不再赘述,当异常分数值越大,意味着该像素点越异常,为噪声像素点的可能性越大。
获得连通范围的每个像素点的异常分数值后,则可根据异常指数和噪声程度获得像素点的异常特征程度,具体包括:对于任意信息图像的任意连通范围中的任意像素点,计算该任意像素点的异常指数的归一化值与对应的连通范围的噪声程度的乘积,获得异常特征程度。当该像素点的异常指数越大,且所在连通范围的噪声程度越大,则该像素点的异常特征程度越大,越像噪声点。获取异常特征程度的公式包括:
式中,表示像素点的异常特征程度,/>表示像素点所在连通范围的噪声程度,/>表示像素点的异常指数,/>表示归一化函数。
获得信息图像中像素点的异常特征程度后,则可根据所有信息图像中相同的连通范围内的不同像素点的异常特征程度的差异特征和噪声程度,获得连通范围的综合特征程度;具体包括:将所有信息图像的同一位置的像素点的异常特征程度构建序列,获得异常表征序列,该异常表征序列能够表征同一位置像素点在不同频率的信息图像中的具体数值,更能够明显地反映不同像素点之间的异常特征程度的差异,进而提高判断该连通范围的噪声特征;计算同一连通范围任意两个异常表征序列的皮尔相关系数绝对值的平均值并归一化,获得特征波动程度;当该特征波动程度值越小,意味着该连通范围内的像素点的异常特征程度在不同信息图像中都相似;当该特征波动程度值越大,意味着不同像素点之间的异常特征程度差异越大,则意味着该连通范围内的噪声像素点越多。计算所有信息图像的同一连通范围的噪声程度的平均值并负相关映射,获得区域稳定表征值;当每个信息图像中的该连通范围的噪声程度越大,则意味着区域稳定表征值越小;计算特征波动程度与区域稳定表征值的乘积,获得连通范围的综合特征程度,当该特征波动程度与区域稳定表征值越大,意味着该连通范围的噪声特征越小,则可降低去噪程度。获取综合特征程度的公式包括:
式中,表示连通范围的综合特征程度,/>表示该连通范围任意两个异常表征序列的皮尔相关系数绝对值的平均值,/>表示所有信息图像的该连通范围的噪声程度的平均值,表示以自然常数为底的指数函数。
步骤S4,根据综合特征程度获得连通范围的平滑系数;根据平滑系数通过非局部均值滤波算法对干涉图去噪获得去噪干涉图,根据去噪干涉图监测地质沉陷。
获得综合特征程度值后,可根据综合特征程度值获得连通范围的平滑系数;优选地,获取平滑系数包括:将综合特征程度负相关映射,获得平滑权重,计算预设平滑值与平滑权重的乘积,获得连通范围对应的干涉图的平滑系数;当该连通范围的综合特征程度越大,意味着噪声特征越不明显,则需要适当降低去噪程度,使用更小的平滑系数,避免纹理细节丢失;当该连通范围的综合特征程度越小,意味着噪声特征越明显,需要适当提高去噪程度,使用更大的平滑系数,保证去噪效果。获得不同连通范围的平滑系数后,则可根据不同的平滑系数通过非局部均值滤波算法对干涉图中不同的连通范围进行去噪获得去噪干涉图,需要说明的是,非局部均值滤波算法属于现有技术,具体去噪过程不再赘述。
进一步地,根据去噪干涉图进行相位解包裹处理,并结合现有的地形模型将地形相位去除,进而得到对应的位移场。根据获取的位移场来确定塌陷区域范围、沉陷速度等信息,实现采煤塌陷区地质沉陷监测。至此通过分析干涉图中不同区域的噪声特征进行自适应去噪,提高了去噪效果,保证了地质沉陷监测的准确性。
综上所述,本发明实施例提供了一种基于D-insar的采煤区沉陷监测方法;根据干涉图中的灰度差异特征获得灰度频率分布图和波动值区间;对干涉图分解获得不同的信息图像;根据波动值区间获得信息图像不同的连通范围;根据连通范围之间的距离特征获得噪声程度;根据连通范围中的灰度差异特征获得像素点的异常指数;根据异常指数和噪声程度获得像素点的异常特征程度;根据所有信息图像不同像素点的异常特征程度和噪声程度,获得连通范围的综合特征程度和平滑系数。本发明根据自适应的平滑系数通过非局部均值滤波算法去噪,提高了去噪和监测地质沉陷的准确性。
本发明还提出了一种基于D-insar的采煤区沉陷监测系统,包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序,处理器执行计算机程序实现任意一项一种基于D-insar的采煤区沉陷监测方法的步骤。
需要说明的是:上述本发明实施例先后顺序仅仅为了描述,不代表实施例的优劣。在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中,多任务处理和并行处理也是可以的或者可能是有利的。
本说明书中的各个实施例均采用递进的方式描述,各个实施例之间相同相似的部分互相参见即可,每个实施例重点说明的都是与其他实施例的不同之处。
Claims (8)
1.一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述方法包括以下步骤:
获取采煤区的干涉图;根据所述干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图;根据所述灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间;
对所述干涉图通过EMD分解算法获得不同的信息图像;根据所述波动值区间获得所述信息图像的不同的连通范围;根据信息图像中所述连通范围之间的距离特征获得连通范围的噪声程度;
根据所述连通范围中的灰度差异特征获得像素点的异常指数;根据所述异常指数和所述噪声程度获得像素点的异常特征程度;根据所有信息图像中相同的连通范围内的不同像素点的异常特征程度的差异特征和噪声程度,获得所述连通范围的综合特征程度;
根据所述综合特征程度获得连通范围的平滑系数;根据所述平滑系数通过非局部均值滤波算法对干涉图去噪获得去噪干涉图,根据去噪干涉图监测地质沉陷;
所述根据所述干涉图中像素点与预设范围内其他像素点的灰度差异特征获得灰度频率分布图的步骤包括:
将干涉图中任意像素点与预设范围内任意其他像素点构建二元组,并计算所述二元组的灰度差值绝对值,获得任意像素点的邻域灰度差;根据干涉图中所有像素点的邻域灰度差构建横坐标为邻域灰度差、纵坐标为邻域灰度差的频率的直方图,获得干涉图的灰度频率分布图;
所述根据所述灰度频率分布图中灰度差异的分布特征获得干涉图的波动值区间的步骤包括:
计算所述灰度频率分布图中任意邻域灰度差对应的所述二元组的种类数与所有二元组种类数的比值,获得二元组种类占比,计算所述任意邻域灰度差对应的频率与所述二元组种类占比的乘积并归一化,获得波动权重,计算所述波动权重与所述任意邻域灰度差的乘积,获得波动范围值;
将所述灰度频率分布图中的所述邻域灰度差从小到大排序,计算从正序开始的预设第一数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最小值;计算从倒序开始的预设第二数量个邻域灰度差对应的波动范围值的平均值并向下取整,获得波动最大值;根据所述波动最小值和所述波动最大值构建区间,获得所述干涉图的波动值区间。
2.根据权利要求1所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述根据所述波动值区间获得所述信息图像的不同的连通范围的步骤包括:
对于任意信息图像中的任意像素点,当任意像素点与预设步长内其他像素点的灰度差值绝对值小于所述波动值区间的任意波动值,则将任意像素点与预设步长内其他像素点构建连通域,遍历所述任意信息图像中所有像素点,获得不同的连通域;当任意两个连通域的平均灰度值的差值小于预设差异值,认为所述任意两个连通域属于同类连通域,获得所述任意信息图像在任意波动值下的不同类别的连通域;
对比所有信息图像中根据不同波动值构建的连通域数量,将连通域数量最多的不同连通域作为所有信息图像中不同的连通范围。
3.根据权利要求2所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述根据信息图像中所述连通范围之间的距离特征获得连通范围的噪声程度的步骤包括:
对于所述任意信息图像中任意波动值下的任意类别连通域,计算所述任意类别与其他类别的任意两个连通域的中心像素点之间的欧氏距离,获得连通距离差异,计算连通距离差异的最大类间方差值并负相关映射,获得任意类别连通域的规律度;
计算所述任意信息图像中所有包括同一个所述连通范围的任意类别连通域的规律度的平均值并负相关映射,获得所述任意信息图像中所述连通范围的噪声程度。
4.根据权利要求1所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述根据所述连通范围中的灰度差异特征获得像素点的异常指数的步骤包括:
计算信息图像中任意像素点的预设邻域范围与任意像素点所在连通范围之间的灰度值方差的比值,作为孤立森林算法中的数据特征,对任意像素点所在连通范围通过孤立森林算法计算任意像素点的异常分数值,获得信息图像中任意像素点的异常指数。
5.根据权利要求1所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述根据所述异常指数和所述噪声程度获得像素点的异常特征程度的步骤包括:
对于任意信息图像的任意连通范围中的任意像素点,计算所述任意像素点的异常指数的归一化值与对应的连通范围的所述噪声程度的乘积,获得异常特征程度。
6.根据权利要求1所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述获得所述连通范围的综合特征程度的步骤包括:
将所有信息图像的同一位置的像素点的所述异常特征程度构建序列,获得异常表征序列,计算同一连通范围任意两个异常表征序列的皮尔相关系数绝对值的平均值并归一化,获得特征波动程度;计算所有信息图像的同一连通范围的所述噪声程度的平均值并负相关映射,获得区域稳定表征值;计算所述特征波动程度与所述区域稳定表征值的乘积,获得所述连通范围的综合特征程度。
7.根据权利要求1所述的一种基于D-insar的采煤区沉陷监测方法,其特征在于,所述根据所述综合特征程度获得连通范围的平滑系数的步骤包括:
将所述综合特征程度负相关映射,获得平滑权重,计算预设平滑值与平滑权重的乘积,获得所述连通范围对应的干涉图的平滑系数。
8.一种基于D-insar的采煤区沉陷监测系统,包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序实现如权利要求1-7任意一项所述方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311657974.2A CN117368920B (zh) | 2023-12-06 | 2023-12-06 | 基于D-insar的采煤区沉陷监测方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311657974.2A CN117368920B (zh) | 2023-12-06 | 2023-12-06 | 基于D-insar的采煤区沉陷监测方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN117368920A CN117368920A (zh) | 2024-01-09 |
CN117368920B true CN117368920B (zh) | 2024-02-27 |
Family
ID=89398783
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202311657974.2A Active CN117368920B (zh) | 2023-12-06 | 2023-12-06 | 基于D-insar的采煤区沉陷监测方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN117368920B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117874685B (zh) * | 2024-03-11 | 2024-05-24 | 天津市塘沽第一阀门有限公司 | 一种直流变频电动半球阀的压力数据异常预警方法 |
CN118154456B (zh) * | 2024-05-07 | 2024-08-20 | 宝鸡拓普达钛业有限公司 | 用于钛合金加工刀具模型构建的图像智能去噪方法 |
CN118364235B (zh) * | 2024-06-14 | 2024-09-17 | 陕西意联电气设备有限公司 | 一种基于大数据的吸湿器状态监测方法及系统 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103760619A (zh) * | 2014-01-07 | 2014-04-30 | 中国神华能源股份有限公司 | 煤田火区的监测方法和装置 |
CN106226764A (zh) * | 2016-07-29 | 2016-12-14 | 安徽理工大学 | 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法 |
CN114442094A (zh) * | 2022-02-10 | 2022-05-06 | 中国地质科学院岩溶地质研究所 | 一种地表形变预测方法及系统 |
CN115980745A (zh) * | 2022-11-21 | 2023-04-18 | 中国地质大学(北京) | 地下开采沉陷区智能识别方法、装置、电子设备及存储介质 |
CN116026283A (zh) * | 2023-01-05 | 2023-04-28 | 中国矿业大学 | 一种煤矿充填开采地表减沉效果监测与评价方法 |
CN116092015A (zh) * | 2023-04-06 | 2023-05-09 | 安徽乾劲企业管理有限公司 | 一种道路施工状态监测方法 |
CN116609781A (zh) * | 2023-05-25 | 2023-08-18 | 北京理工大学 | 多星数据联合的北斗InSAR DEM误差补偿方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
AUPR187100A0 (en) * | 2000-12-04 | 2001-01-04 | Cea Technologies Inc. | Slope monitoring system |
-
2023
- 2023-12-06 CN CN202311657974.2A patent/CN117368920B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103760619A (zh) * | 2014-01-07 | 2014-04-30 | 中国神华能源股份有限公司 | 煤田火区的监测方法和装置 |
CN106226764A (zh) * | 2016-07-29 | 2016-12-14 | 安徽理工大学 | 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法 |
CN114442094A (zh) * | 2022-02-10 | 2022-05-06 | 中国地质科学院岩溶地质研究所 | 一种地表形变预测方法及系统 |
CN115980745A (zh) * | 2022-11-21 | 2023-04-18 | 中国地质大学(北京) | 地下开采沉陷区智能识别方法、装置、电子设备及存储介质 |
CN116026283A (zh) * | 2023-01-05 | 2023-04-28 | 中国矿业大学 | 一种煤矿充填开采地表减沉效果监测与评价方法 |
CN116092015A (zh) * | 2023-04-06 | 2023-05-09 | 安徽乾劲企业管理有限公司 | 一种道路施工状态监测方法 |
CN116609781A (zh) * | 2023-05-25 | 2023-08-18 | 北京理工大学 | 多星数据联合的北斗InSAR DEM误差补偿方法 |
Non-Patent Citations (4)
Title |
---|
Monitoring Extractive Activity-Induced Surface Subsidence in Highland and Alpine Opencast Coal Mining Areas with Multi-Source Data;Shuqing Wang et al.;Remote Sensing;第14卷(第14期);第1-12页 * |
基于D-InSAR的矿区地面沉降监测;刘付刚等;黑龙江科技大学学报;第27卷(第3期);第265-269页 * |
姚薇 ; 钱玲玲 ; .矿山遥感图像自适应加权改进中值滤波算法.金属矿山.2016,(第04期),第101-105页. * |
矿区塌陷区遥感影像改进自适应维纳滤波算法;冯丽慧;;金属矿山(第07期);第151-154页 * |
Also Published As
Publication number | Publication date |
---|---|
CN117368920A (zh) | 2024-01-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN117368920B (zh) | 基于D-insar的采煤区沉陷监测方法及系统 | |
CN115829883B (zh) | 一种异性金属结构件表面图像去噪方法 | |
CN106296655B (zh) | 基于自适应权值和高频阈值的sar图像变化检测方法 | |
CN109919870B (zh) | 一种基于bm3d的sar图像相干斑抑制方法 | |
CN103226832B (zh) | 基于光谱反射率变化分析的多光谱遥感影像变化检测方法 | |
CN104200471A (zh) | 基于自适应权值图像融合的sar图像变化检测方法 | |
CN115512231B (zh) | 适用于国土空间生态修复的遥感解译方法 | |
CN103871039B (zh) | 一种sar图像变化检测差异图生成方法 | |
CN111008585A (zh) | 基于自适应分层高分辨sar图像的舰船目标检测方法 | |
CN102360503B (zh) | 基于空间贴近度和像素相似性的sar图像变化检测方法 | |
CN106940782B (zh) | 基于变差函数的高分sar新增建设用地提取软件 | |
CN117272210A (zh) | 一种建筑施工异常隐患数据检测方法及系统 | |
CN116309757A (zh) | 基于机器视觉的双目立体匹配方法 | |
CN113628234B (zh) | 基于综合邻域信息的显著性极化sar图像变化检测方法 | |
CN104680536B (zh) | 利用改进的非局部均值算法对sar图像变化的检测方法 | |
CN111290021A (zh) | 基于梯度结构张量谱分解的碳酸盐岩缝洞增强识别方法 | |
CN115937231B (zh) | 一种频谱结构约束的红外图像迭代去噪方法及系统 | |
CN118035664B (zh) | 基于多维数据的地质信息数据分析决策方法及系统 | |
CN110929598B (zh) | 基于轮廓特征的无人机载sar图像匹配方法 | |
CN115840205A (zh) | 一种基于激光雷达技术的地形面积计量方法和系统 | |
CN113837074B (zh) | 结合后验概率和空间邻域信息的遥感影像变化检测方法 | |
CN112669332B (zh) | 一种基于双向局部极大值和峰值局部奇异性判断海天条件和检测红外目标的方法 | |
CN117830415A (zh) | 一种车载大尺寸光电玻璃盖板加工定位方法 | |
CN117475154A (zh) | 一种基于优化Canny边缘检测的仪表图像识别方法及系统 | |
CN109697474B (zh) | 基于迭代贝叶斯的合成孔径雷达影像变化检测方法 |
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 |