CN109409230A - 一种空地数据互增强的震前热异常自动分级方法 - Google Patents
一种空地数据互增强的震前热异常自动分级方法 Download PDFInfo
- Publication number
- CN109409230A CN109409230A CN201811129842.1A CN201811129842A CN109409230A CN 109409230 A CN109409230 A CN 109409230A CN 201811129842 A CN201811129842 A CN 201811129842A CN 109409230 A CN109409230 A CN 109409230A
- Authority
- CN
- China
- Prior art keywords
- image
- surface temperature
- heat
- extracted
- abnormality
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 49
- 230000005856 abnormality Effects 0.000 claims abstract description 97
- 239000000284 extract Substances 0.000 claims abstract description 45
- 230000005855 radiation Effects 0.000 claims abstract description 36
- 238000010276 construction Methods 0.000 claims description 16
- 230000002159 abnormal effect Effects 0.000 claims description 11
- 238000002372 labelling Methods 0.000 claims description 10
- 238000005457 optimization Methods 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 4
- 230000035699 permeability Effects 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 5
- 238000000605 extraction Methods 0.000 description 20
- 238000004422 calculation algorithm Methods 0.000 description 13
- 238000002834 transmittance Methods 0.000 description 13
- 230000008859 change Effects 0.000 description 7
- 230000009466 transformation Effects 0.000 description 7
- 230000000452 restraining effect Effects 0.000 description 6
- 238000001228 spectrum Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000002708 enhancing effect Effects 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000002968 anti-fracture Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 229940068517 fruit extracts Drugs 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 239000005433 ionosphere Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000009738 saturating Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
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
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/01—Measuring or predicting earthquakes
-
- 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/25—Determination of region of interest [ROI] or a volume of interest [VOI]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Multimedia (AREA)
- Theoretical Computer Science (AREA)
- Remote Sensing (AREA)
- Acoustics & Sound (AREA)
- Astronomy & Astrophysics (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Radiation Pyrometers (AREA)
- Geophysics And Detection Of Objects (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种空地数据互增强的震前热异常自动分级方法，包括以下步骤：对遥感数据进行预处理，获取兴趣区域待热红外异常提取的亮温图像、背景场图像；获取对应日期的地表温度数据，采用克里金插值进行处理；进行地表温度反演；对反演后的地表温度与实际地表温度进行权重配比，得到均衡地表温度，差值计算待提取亮温图像与背景场图像对应的每一个像元值，获取热红外异常提取图像；单位概率化断裂带和震源统计信息，作为权重约束，同时考虑断裂带周围的辐射域，作为权重约束，处理热红外异常提取图像,生成震前热红外异常优化图像，同时消除其他热红外辐射异常对热红外异常提取结果的干扰。通过正态分布将震前热红外异常优化结果分为多个级别。
Description
技术领域
本发明属于地理空间信息系统技术领域，主要涉及一种震前热红外异常的分级方法。
背景技术
随着卫星热红外遥感技术的快速发展，热红外辐射数据获取的精度和可靠性提高，为震前热红外异常提取和地震监测提供了有效可靠的数据支撑，卫星热红外异常提取技术业已成为地震监测预报的一种重要手段。实际的卫星热红外辐射数据往往会受到电离层扰动和气候瞬时变化的影响，然而现有的卫星热红外异常提取方法(功率谱相对变化法，偏移指数法，改进的RST法等)采用多年均值作为背景场，无法充分考虑到气候瞬时变化的影响；并且现有的算法没有注意绝大多数有效的热红外异常也大都发生在断裂带附近，没有充分利用断裂带的约束，卫星热红外异常提取方法提取出大面无效积异常，极大的影响了异常提取的可靠性。因此，如何通过有效的地面数据约束，准确、高效的提取震前热红外异常成为当前亟待解决的难题。
目前国内许多机构都对震前外热红外异常提取展开深入的研究，但是已有的研究主要侧重于背景场的构建和差分提取异常，将多年同期的数据求平均作为背景场以增加算法的可靠性，仅仅利用卫星热红外影像数据进行热红外异常提取，因此，往往存在已下两点不足：
(1)忽略了电离层和气候变化的瞬时性，同时，卫星热红外影像数据不能和地面数据形成有效的验证，异常提取结果的可靠性和准确性无法满足实际应用的需求；
(2)无法解决算法在非地震区域提取出大面积的疑似异常(与地震和构造活动无关的地表热红外异常)，对热红外异常有效提取与分级造成极大干扰。
地震前中期，震块之间相互作用产生大量电荷使得周围低空气体被电离，造成了热红外异常现象，在电离低空气体向外辐射的过程中会被实时变化的大气和地表温度影响。被电离的气体一般情况下都分布在断裂带周围，受到断裂带的约束；同时在在辐射传输过程中会受到大气和地表温度的约束。如果我们能充分考虑这些约束关系，将其引入到热红外异常提取过程中，将大大改善震前热红外提取的效果以及异常提取的准确度。
目前，现有的震前热红外异常提取算法大多数是只考虑卫星热红外图像变化，现有的算法可以看成是某天亮温图像或者是亮温图像处理结果与背景场(历年同期亮温图像或者是亮温图像处理结果求平均)差分，得出的计算结果作为热红外异常提取结果。现有的算法一般可分为两类策略：基于空间域的异常提取策略和基于频率域的异常提取策略：
(1)基于空间域的异常提取策略目前基于空间域的异常提取方法主要有RST法和偏移指数法，RST法是一种基于背景场分析的异常提取方法，结合地表温度时间序列的周期性特征，引入谐波分析，采用傅里叶逼近的方法拟合地表温度离散时序，从中提取其年趋势，建立一个动态的、同时包含局地信息和年际特征的地表温度背景场；将其引入RST 模型，基于“kσ”准则识别地震热异常信息。偏移指数法和RST法的策略类似，只是在构建背景场的时候采用简单的历年同期均值。这些算法只是在卫星热异常图像的基础上进行不同的差分策略，没有考虑到有效的热红外异常往往具有区域约束性，也没考虑到非地震因素造成的热红外异常。
(2)基于频率域的异常提取策略
基于频率域的异常提取策略主要有功率谱相对变化法，功率谱相对变化策略采用7阶小波变换部分减去2阶小波变换部分，保留中间频带部分，达到消除温度场、日变温度场、气流和地形等影响。由于热红外辐射变化呈周期性，短邻地震异常出现时间一般在1个月到3个月之间，对每一个像元建立一个按时间序列排列的亮温数据组，通过傅里叶变换提取地震前后热红外变化的优势频率和幅值，并将对应的幅值做相对处理，将频率幅值变化较大的作为异常提取结果。实际上该方法也属于基于背景场分析的异常提取方法，功率谱相对变化法的数据输入是是亮温数据不同阶小波变换差分结果，小波变换属于通用滤波模型，极其容易平滑震块之间的热红外辐射异常，导致一部分异常无法被提取。加之，缺乏地面数据的支撑验证，无法保证异常提取的可靠性。
发明内容
针对卫星热红外图像提取震前热异常过程中，由于地表温度、日变温度场、气流实时变化，以及缺乏地面数据的约束，导致震前热红外异常提取结果往往存在大范围非地震因素的异常，即使有效异常也容易出现偏差等严重问题，本发明将均衡地表温度引入震前热红外异常提取过程中，使待提取亮温图像与背景场亮温图像对应均衡地表温度的差值结果作为权重约束，提取异常结果；将断裂带和地震震源信息引入，计算单位化的地震概率并权重化，同时将辐射域权重化，在双重权重约束的基础上，减弱大气扰动的影响，克服大面积无关热红外辐射的干扰，并对上一步异常提取结果自动优化分级，建立一套科学合理的震前热红外异常自动分级实施方案。具体而言，本发明是一种多种地面数据约束增强的震前热异常自动分级方法，其各模块关系如图1所示，本发明主要包含两个主要步骤：1)具有地表温度约束的震前热红外异常提取算法；2)顾及断裂带约束的震前热红外异常自动分级方法，算法整体流程如图2所示。
本发明解决其技术问题所采用的技术方案是：一种空地数据互增强的震前热异常自动分级方法，包括以下步骤：
对遥感数据进行预处理，获取兴趣区域待提取的亮温图像，通过前置滑窗选择与待提取亮温图像方差最小的待选背景场亮温图像作为背景场图像；
获取待提取亮温图像和背景场图像对应日期的地表温度数据，采用克里金差值进行处理；对处理后的待提取亮温图像和背景场图像进行地表温度反演；
对反演得到的地表温度Tf与实际地表温度Tp进行权重配比，得到均衡地表温度，通过计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值，获取热红外异常提取图像；
处理热红外异常提取图像，生成震前热红外异常优化图像，进而通过正态分布将震前热红外异常优化图像分为多个级别。
所述对处理后的待提取亮温图像和背景场图像进行地表温度反演如下：
Tf＝A0+A1T31-A1T32
其中，Tf是克里金差值处理后待提取亮温图像Nlst或背景场图像Nf反演后的地表温度，T31和T32分别是MODIS第31和32波段的亮温图像，待提取亮温图像Nlst和背景场图像Nf采用的是MODIS的31波段，故进行地表温度反演时 T31＝Nlst或Nf；A0、A1和A2是劈窗法的参数，分别定义如下：
其中，a31、b31、a32、b32为常数，中间参数Ci和Di通过下式得到：
其中，i是MODIS的第31波段和32波段，分别表示为i＝31或i＝32；τi(θ)是视角为θ的大气透过率；εi是波段为i的地表比辐射率。
所述均衡地表温度T＝(Tp×3/5)+(Tf×2/5)；
当Tp＝Tlst，T31＝Nlst时，Tlst为克里金差值处理后的待提取亮温图像Nlst对应的均衡地表温度，待提取亮温图像得到的均衡地表温度为：Tbalancelst＝(Tlst×3/5)+ (Nlst×2/5)；
当Tp＝Tfl，T31＝Nf时，Tfl为克里金差值处理后背景场亮温图像对应的均衡地表温度，背景场亮温图像得到的均衡地表温度为：Tbalancef＝(Tfl×3/5)+(Nf× 2/5)；
待提取亮温图像与背景场图像对应的均衡地表温度差值ΔT＝Tbalancelst－Tbalancef。
所述通过计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值，获取热红外异常提取图像包括以下步骤：
计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值的差分绝对值ΔN；
当ΔT大于阈值时，异常提取值为N＝0；
当ΔT小于阈值时，异常提取值为N＝((阈值-ΔT)/阈值)×ΔN；
遍历提取亮温图像与背景场图像的所有像元，将每个像元对应的异常提取值组成矩阵构成热红外异常提取图像。
所述处理热红外异常提取图像，生成震前热红外异常优化图像包括以下步骤：
步骤2.1、获取带有地理位置信息的矢量断裂带数据和对应的地震震源信息，得到单位面积下地震次数PER；
步骤2.2、所有断裂带都标上相应的PERj，给每个断裂带Areaj赋相应的权重分级标签Hqj，权重分级标签分为多个等级；j＝1…S，S为矢量断裂带数据个数；
步骤2.3、将每个断裂带的辐射域分为若干层级，每一层级辐射域都分配对应的权重分级标签Hrji；
步骤2.4，对于获取的热红外异常提取图像，逐像元计算震前热红外异常优化值Nconstruction：
如果热红外异常提取图像的像元位于断裂带内，则Nconstruction＝N×Hq；
如果热红外异常提取图像的像元位于辐射域内，则Nconstruction＝N×Hq×Hr；
其他，Nconstruction＝0；
其中，Hq为断裂带权重标签，Hr为辐射域权重标签，N为异常提取值；
步骤2.5，返回步骤2.2直至遍历所有红外异常提取图像，得到震前和热红外异常优化图像；统计震前热红外异常优化图像非零值，求解图像的期望值μ和标准差σ，得到震前热红外异常优化图像非零值正态分布函数；根据期望值μ和标准差σ将震前热红外异常分为若干个级别L，0值归类为正常，其他采用正态分布标记等级。
所述其他采用正态分布标记等级具体如下：
本发明具有以下有益效果及优点：
从卫星热红外图像中提取，如何消除地表温度和大气变化的影响影响到震前热红外异常提取结果的可靠性和稳定性。由于卫星热红外图像提供的信息十分有限，只有加入地表温度约束，才能一定一定程度上克服大气变化带来的影响。然而现有的算法通过小波变换、多期影像傅里叶变换等通用模型，无法达到精准消除大气和温度实时变化带来的影响。因此本发明通过亮温数据均方差最小原则选择背景场，并按权重配比实际地表温度与反演后地表温度，作为权重约束热红外异常提取，从而降低大气变化对地表温度的影响，提高震前热红外异常提取的精度。
在获取震前热红外异常提取结果后，该结果中依然会存在大范围地震无关的热红外异常，现有的“一刀切”方法通常只处理断裂带带区域范围内，而没有充分考虑地震断裂带周围异常情况。本发明单位概率化断裂带和震源统计信息，作为权重约束，同时充分考虑断裂带周围的辐射域，也作为分级权重约束，消除其他热红外辐射异常对热红外异常提取结果的干扰。
附图说明
图1为本发明各模块之间关系图。
图2为本发明算法总体流程图。
图3为具有地表温度约束的震前热红外异常提取算法流程图。
图4为顾及断裂带约束的震前热红外异常自动分级方法流程图。
图5为震前热红外异常分级专题图。
具体实施方式
下面结合实施例对本发明做进一步的详细说明。
具有地表温度约束的震前热红外异常提取方法，具体步骤包括：通过前置滑窗策略确定背景场，并采用克里金内插规则化地表温度数据；利用劈窗法实现地表温度反演，进一步，通过权重化配比反演后的地表温度和实际地表温度，得到均衡地表温度；计算出待提取亮温图像与背景场亮温图像对应均衡地表温度的差值，将该差值结果作为权重约束，得到正则化的热红外异常提取图像，存入内存应用于后续的震前热红外异常自动分级。一种顾及断裂带约束的震前热红外异常自动分级方法，具体步骤包括：依据断裂带和地震震源信息，计算单位化的地震概率，把该概率作为震前热红外异常自动分级的权重约束；采用反距离权的策略，将断裂带辐射域分为三个等级，并作为震前热红外异常自动分级的权重约束；通过两重权重约束，优化上一步获取的热红外异常提取图像, 得到最终的增强结果，最后，采用正态分布策略，将增强后的震前热红外异常分为五个等级L{严重，重度，中度，轻度，正常}。
一种具有地表温度约束的震前热红外异常提取算法，如图3所示，主要步骤包括：
步骤1.1，通过相应的预处理，获取特定区域待热红外异常提取的亮温图像 Nlst，采用前置滑窗的策略(滑动窗口长度设置为w＝3)，选择与待提取亮温图像方差最小的亮温图像作为背景场Nf。
步骤1.2，获取亮温图像(待提取亮温图像和背景场图像)对应日期的地表温度数据，通过克尤金内插实现零散的地表温度站点数据规则化，规则化后的地表温度数据区域范围与空间分辨率同亮温图像保持一致。
步骤1.3，采用劈窗法对待提取亮温图像和背景场图像进行地表温度反演。通过MODIS第2和第19波段估算大气水汽含量，构建水汽含量与大气透过率方程，计算出大气透过率；通过MODIS第31和第32波段估算地表比辐射率，最后利用大气透过率和地表比辐射率实现地表温度高效反演。
步骤1.4，反演后的地表温度Tf与实际地表温度Tp进行2比3的权重配比，可得到：均衡地表温度T＝(Tp×3/5)+(Tf×2/5)，并计算出待提取亮温图像与背景场图像对应的均衡地表温度差值ΔT。
步骤1.5，通过差分计算待提取亮温图像每一个像元值与背景场亮温图像对应的每一个像元值，得到二者差分绝对值ΔN，当ΔT大于阈值1.5时，异常提取值为N＝0；当小于1.5时，异常提取值为N＝((1.5-ΔT)/1.5)×ΔN。
步骤1.6，重复步骤1.5直至遍历所有像元，将所有计算结果正则化后，获取最终的热红外异常提取图像。
顾及断裂带约束的震前热红外异常自动分级方法，如图3所示，主要步骤包括：
步骤2.1，获取带有地理位置信息的矢量断裂带数据和对应的地震震源信息，并统计每一个断裂带中发生地震的次数；将每一个断裂带区域内发生地震的次数除以对应断裂带UTM投影下的面积，得到单位面积下地震次数PER。
步骤2.2，重复步骤2.1直至断裂带遍历结束，所有断裂带都标上相应的编号i，通过相应的阈值约束，给每个断裂带Areaj赋相应的权重分级标签Hqj，权重分级标签分为三个等级。
步骤2.3，考虑到地震发生的可能性同断裂带的距离成反比，将每一个断裂带辐射范围沿法线方向外扩张150千米，并称之为辐射域；辐射域分为三层级，每一层级辐射域的都分配对应的权重分级标签Hrji；如果相邻断裂带的辐射域存在重叠的情况，取权重最大者。
步骤2.4，计算带有多重约束的震前热红外异常提取值Nconstruction，如果热红外异常提取图像的像元位于断裂带内，则Nconstruction＝N×Hq；如果热红外异常提取图像的像元位于辐射域内，则Nconstruction＝N×Hq×Hr；否则，Nconstruction＝0；
步骤2.5，重复步骤2.2至2.4直到处理完所有红外异常提取图像，并生成带有多重约束的震前热红外异常提取图像，进而采用正态分布的策略，将震前热红外异常分为五个级别L{严重，重度，中度，轻度，正常}。
本发明的总体技术流程如图2所示，将首先对MODIS数据进行预处理得到热红外图像，随后采用一套具有多重地表数据权重约束策略进行震前热红外自动分级，具体包括：1)通过预处理，生成待提取亮温数据，并依据前置滑窗策略和均方差最小原则确定背景场；2)引入地表温度数据，权重化配比反演后的地表温度和实际地表温度，计算出均衡地表温度，并将此作为约束条件计算出热红外异常提取结果；3)充分考虑断裂带矢量数据和震源信息，权重化每个断裂带地震次数，同时构建辐射域并权重分级，作为双重约束条件优化热红外异常提取结果；4)采用正态分布策略进行异常分级，得到震前热红外分级结果。为了理解本发明，这里首先整体描述自动化的实现流程，随后将对本发明的两个关键步骤进行详细说明。
1.预处理。MODIS一级数据通常都是不能够直接使用，需要对数据进行相应的波段筛选、拼接、裁剪和重投影，生成栅格数据后方可用于后续计算。我们采用MRT(MODISReprojection Tool)对原始MODIS数据进行预处理，每次参与计算的数据裁剪范围都保持一致，同时将处理好的数据统一投影到统一坐标系下(地理坐标系)，预处理后的栅格数据统一成为亮温图像Nlst。
2.引入地表温度数据。由于站点监测数据是呈离散分布，获取的地表温度站点监测数据，不能和对应的亮温图像一一对应，需要通过克里金内插实现地表温度站点数据的规则化，规则化后的地表温度数据区域范围与空间分辨率同亮温图像保持一致。
3.利用规则化后的地表温度数据，采用本发明中具有地表温度约束的震前热红外异常提取算法，提取震前热红外异常提取结果。
4.引入获取带有地理位置信息的矢量断裂带数据和对应的地震震源信息。
5.实现震前热红外异常自动分级。采用本发明中顾及断裂带约束的震前热红外异常自动分级方法，将步骤4获取的信息作为约束条件，实现震前热红外异常提取结果的优化与增强，同时将震前热红外异常分为五个级别L{严重，重度，中度，轻度，正常}。
6.制作异常分级结果的专题图。通过五种不同的颜色递进标注五个震前热红外异常级别，并制作专题图，如图5所示。
下面将详细介绍本发明中的具有地表温度约束的震前热红外异常提取算法，如图3所示，主要步骤包括：
步骤1.1，通过相应的预处理，获取特定区域待热红外异常提取的亮温图像 Nlst，假设选择特定区域2018年1月11日的亮温图像Nk，下标k为该日(2018 年1月11日)的数值标号，作为待异常提取亮温图像Nlst，采用前置滑窗的策略(滑动窗口长度设置为w＝3)，确定待选背景场图像序列为1月8日-1月10 日的亮温图像集：Nfc＝{Nk-i|i＝(1,2,3)}。下标k为该日(2018年1月11日)的数值标号。依次计算待热红外异常提取的亮温图像Nlst与待选背景场图像序列Nfc的方差：使得背景场图像与待提取图像能量差最小，其中m,n分别为亮温图像的行列数，σ＝{σ1,σ2,σ3}，下标i为待选背景场图像序号，下标j为像元序号，Nkj为待异常提取图像对应序号j的像元值， N(k-i)j为待选背景场图像对应序号j的像元值，最后选择与待提取亮温图像方差最小σmin＝min(σ)的亮温图像作为背景场图像Nf。
步骤1.2，采用克里金内插法(Kriging interpolation method)实现待提取亮温图像和背景场图像对应日期的地表温度数据规则化。具体步骤如下：
步骤1.2.1，获取待提取亮温图像经纬度范围(即左上角经纬度坐标与右下角经纬度坐标)r和空间分辨率res。
步骤1.2.2，采用克里金内插法对离散地表温度站点数据进行规则化内插，内插结果与待提取亮温图像的经纬度范围与空间分辨率保持一致。
步骤1.2.3，得到规则化后的地表温度栅格数据，下文称为地表温度图像，其中规则化后的待提取亮温图像对应地表温度栅格数据为Tlst；规则化后的背景场亮温图像对应地表温度栅格数据为Tfl。
步骤1.2.4，采用步骤1.2.1-步骤1.2.3的方法，完成背景场亮温图像Nf对应的地表温度图像的规则化。
步骤1.3，采用劈窗法对待提取亮温图像和背景场亮温图像进行地表温度反演。通过MODIS第2和第19波段估算大气水汽含量，构建水汽含量与大气透过率方程，计算出大气透过率；通过MODIS第31和第32波段估算地表比辐射率，最后利用大气透过率和地表比辐射率实现地表温度反演,待提取亮温图像 Nlst和背景场图像Nf反演公式相同，对应的反演参数一致，故采用统一的反演公式，具体公式为：
Tf＝A0+A1T31-A1T32
其中，Tf是反演后的地表温度，T31和T32分别是MODIS第31和32波段的亮温值，Nlst和Nf采用的是MODIS的31波段，故进行地表温度反演时T31＝Nlst或 Nf；，T32获取方式与T31相同(预处理同一天的MODIS数据，获得与T31经纬度范围和分辨率一致的亮温图像)，此处T32作为已知值；A0,A1和A2是劈窗法的参数，分别定义如下：
其中常亮的数值为：a31＝-64.603，b31＝0.440817，a32＝-68.72575， b32＝0.473453，中间参数Ci和Di无法直接获得需要通过进一步计算，其公式推导如下：
其中，i是指MODIS的第31波段和32波段，分别i＝31或i＝32；τi(θ)是视角为θ的大气透过率；εi是波段为i的地表比辐射率。从上述推导可得出结论：地表温度反演可以分解两大步，1)大气透过率τi(θ)计算；2)地表比辐射率εi估计。其具体实施步骤如下：
步骤1.3.1，获取亮温数据对应时间的MODIS第2和第19波段数据，并将其裁剪为同待提取亮温图经纬度范围一致。
步骤1.3.2，通过MODTRAN的大气辐射传输模拟，确定大气透过率τi(θ)和大气水分含量w的关系，大气透过率估计方程如表1所示：
表1大气透过率与大气水分含量对应表
其中，SEE为估计标准误差，R2为相关指数，FE为误差自由度，查表求得大气透过率τi(θ)。
步骤1.3.3，大气透过率受传感器视角和大气剖面温度的影响，因此，需要对这两个因素进行校正，其校正模型如表2所示：
表2大气透过率的温度校正模型函数
对比校正模型表格，可以直接求出31波段和32波段大气透过率τi(θ)的准确值。
步骤1.3.4，通过(HanMa,Shun linLiang.Simultaneous inversion of multipleland surface parameters from MODIS optical-thermal observations In:ISPRSJournal of Photogrammetry and Remote Sensing,128(2017)240-254)提供的策略完成地表比辐射率εi估计。
步骤1.3.5，将步骤1.3.3和步骤1.3.4求得的结果地表温度反演公式，解算得到反演后的地表温度栅格数据(下文称为反演后地表温度图像)。
步骤1.3.6，依次执行步骤1.3.1-步骤1.3.5，得到待提取亮温图像和背景场亮温图像对应的地表温度反演后地表温度图像。
步骤1.4，传统的基于空间域的异常提取方法，直接输入的数据通常为未经处理的亮温影像，参与计算亮温图像没有充分考虑大气与地表温度的影像导致异常提取结果不稳定，因此需要通过反演地表温度消除大气影响，权重配比反演后地表温度和实际地表温度提升数据可靠性。计算反演后的地表温度Tf与实际地表温度Tp进行2比3的权重配比，可得到：均衡地表温度T＝(Tp×3/5)+(Tf×2/5)，当Tp＝Tlst，T31＝Nlst时(Tlst为规则化后待提取亮温图像对应的均衡地表温度)，待提取亮温图像得到的均衡地表温度为：Tbalancelst＝(Tlst×3/5)+(Nlst×2/5)；当Tp＝Tfl，T31＝Nf时(Tfl为规则化后背景场亮温图像对应的均衡地表温度)，背景场亮温图像得到的均衡地表温度为：Tbalancef＝(Tfl×3/5)+(Nf×2/5)；并得出待提取亮温图像与背景场图像对应的均衡地表温度差值ΔT＝Tbalancelst－Tbalancef。
步骤1.5，现有热红外异常提取是通过待异常提取图像与背景场图像的差分实现，本发明引入均衡地表温度的概念，同时增加背景场图像和待提取亮温图像对应的均衡地表温度差值ΔT作为热红外异常提取约束条件，通过差分计算待提取亮温图像每一个像元值与背景场亮温图像对应的每一个像元值，得到二者差分绝对值ΔN，当ΔT大于阈值1.5时，异常提取值为N＝0；当小于1.5时，异常提取值为N＝((1.5-ΔT)/1.5)×ΔN，公式如下所示：
步骤1.6，重复步骤1.5直至遍历所有像元，将所有计算结果正则化后，获取热红外异常提取图像：按照热红外异常提取公式N＝((阈值-ΔT)/阈值)×ΔN计算每个一个异常值，实际上亮温图像与背景场图像都可以看做一个矩阵，逐像元计算(N＝((阈值-ΔT)/阈值)×ΔN)就是这两个图像对应行列号矩阵单元值的运算；每个对应像元的计算得到一个异常提取值，当把所有的像元遍历结束，这两个图像对应的每一个像元都会得到一个异常提取结果，这些结果组成矩阵和亮温图像与背景场图像大小一致，称为异常提取图像。
下面将详细介绍本发明中的顾及断裂带约束的震前热红外异常自动分级方法，如图4所示，主要步骤包括：
步骤2.1，获取带有地理位置信息(经纬度坐标)的矢量断裂带数据(每一个矢量断裂带的数据形式为依次存储的点列经纬度坐标)和对应的地震震源信息(100年内有记录大于3.0的地震震源地理坐标和里氏震级)，并统计每一个断裂带中发生地震的次数；将每一个断裂带区域内发生地震的次数除以对应断裂带在UTM投影(Universal TransverseMercator Projection)下的面积，得到单位面积下地震次数PER。具体步骤如下：
步骤2.1.1，统计地震震源(里氏震级大于5)总数M，并将震源信息依次按经度递增原则(第二判别条件为纬度递增)排序{Eq1,Eq2,…,EqM}，震源对应的经纬度坐标为：Eqi＝(Xi,Yi)。i＝1…M。
步骤2.1.2，统计矢量断裂带数据个数S，断裂带按照面积递增原则排序{Fz1,Fz2,…,FzS},每一个矢量断裂带顶点的经纬度坐标均存入数组Fzj＝{(Xj,Yj)|j ≤k}，k为每个断裂带顶点坐标个数。j＝1…S。
步骤2.1.3，断裂带Fzi对应地震震源Eqi的统计数为S组{Num1,Num2,…, NumS}，依次遍历地震震源信息Eqi，对于每一个Eqi，也依次判断震源Eqi是否落入断裂带Fzj内，如果是则Numj增一，否则删除该震源信息，直到地震震源信息遍历结束,其中(i＝1…M；j＝1…S)。
步骤2.1.4，计算每一个断裂带在带UTM投影下的面积{Area1,Area2,…, AreaS}，单位为平方千米。
步骤2.1.5，统计每一个断裂带区域内发生地震的次数除以对应断裂带UTM 投影下的面积，得到单位面积下地震次数PER，公式如下：
步骤2.2，重复步骤2.1.1-步骤2.1.5直至断裂带遍历结束，所有断裂带都标上相应的PERj，给每个断裂带Areaj赋相应的权重分级标签Hqj，权重分级标签分为三个等级1.0，1.5和2.0，公式如下所示
其中，
步骤2.3，地震发生并非局限在断裂带内，经统计发现，离断裂带越远，发生地震的可能逐级递减，考虑到地震发生的可能性同断裂带的距离成反比，将每一个断裂带辐射范围沿法线方向外扩大150千米(不包含断裂带区域)，称该范围为辐射域；辐射域分为三层级，每一层级辐射域的都分配对应的权重分级标签Hrji，公式如下所示：
如果相邻断裂带的辐射域存在重叠的情况，交叠区域会被该区域对应的两个权重标签：Hqj×Hrji的乘积最大者标记，如果某断裂带辐射域与其他断裂带重叠，重叠区域被断裂带标记。
步骤2.4，优化步骤1.6提取的热红外异常结果，逐像元计算带有单位化的地震概率约束和辐射域权重等多重约束的震前热红外异常优化值Nconstruction，如果热红外异常提取图像的像元位于断裂带内，则Nconstruction＝N×Hq；如果热红外异常提取图像的像元位于辐射域内，则Nconstruction＝N×Hq×Hr；否则，Nconstruction＝0，公式如下：
其中，Hq为断裂带权重标签，Hr为辐射域权重标签，N为异常提取值。
步骤2.5，重复步骤2.4直到优化完所有红外异常提取图像，并生成带有多重约束的震前热红外异常优化图像，统计震前热红外异常优化图像非零值，求解图像的期望值μ和标准差σ，得到震前热红外异常优化图像非零值正态分布函数，将震前热红外异常分为五个级别L{严重，重度，中度，轻度，正常}，0值归类为正常，其他均采用正态分布策略标记等级：震前热红外异常优化图像的像元值在(μ-0.3σ,μ+0.3σ)区间内，归类为正常；在(μ-0.7σ,μ-0.3σ)∪ (μ+0.3σ,μ+0.7σ)区间内，归类为轻度；在(μ-1.2σ,μ-0.7σ)∪(μ+0.7σ,μ+ 1.2σ)区间内，归类为中度；在(μ-1.8σ,μ-1.2σ)∪(μ+1.2σ,μ+1.8σ)区间内，归类为重度；在(μ-2.5σ,μ-1.8σ)∪(μ+1.8σ,μ+2.5σ)区间内，归类为严重；具体归类方式公式如下所示：
Claims (6)
1.一种空地数据互增强的震前热异常自动分级方法，其特征在于包括以下步骤：
对遥感数据进行预处理，获取兴趣区域待提取的亮温图像，通过前置滑窗选择与待提取亮温图像方差最小的待选背景场亮温图像作为背景场图像；
获取待提取亮温图像和背景场图像对应日期的地表温度数据，采用克里金差值进行处理；对处理后的待提取亮温图像和背景场图像进行地表温度反演；
对反演得到的地表温度Tf与实际地表温度Tp进行权重配比，得到均衡地表温度，通过计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值，获取热红外异常提取图像；
处理热红外异常提取图像，生成震前热红外异常优化图像，进而通过正态分布将震前热红外异常优化图像分为多个级别。
2.根据权利要求1所述的一种空地数据互增强的震前热异常自动分级方法，其特征在于，所述对处理后的待提取亮温图像和背景场图像进行地表温度反演如下：
Tf＝A0+A1T31-A1T32
其中，Tf是克里金差值处理后待提取亮温图像Nlst或背景场图像Nf反演后的地表温度，T31和T32分别是MODIS第31和32波段的亮温图像，待提取亮温图像Nlst和背景场图像Nf采用的是MODIS的31波段，故进行地表温度反演时T31＝Nlst或Nf；A0、A1和A2是劈窗法的参数，分别定义如下：
其中，a31、b31、a32、b32为常数，中间参数Ci和Di通过下式得到：
其中，i是MODIS的第31波段和32波段，分别表示为i＝31或i＝32；τi(θ)是视角为θ的大气透过率；εi是波段为i的地表比辐射率。
3.根据权利要求1所述的一种空地数据互增强的震前热异常自动分级方法，其特征在于，所述均衡地表温度T＝(Tp×3/5)+(Tf×2/5)；
当Tp＝Tlst，T31＝Nlst时，Tlst为克里金差值处理后的待提取亮温图像Nlst对应的均衡地表温度，待提取亮温图像得到的均衡地表温度为：Tbalancelst＝(Tlst×3/5)+(Nlst×2/5)；
当Tp＝Tfl，T31＝Nf时，Tfl为克里金差值处理后背景场亮温图像对应的均衡地表温度，背景场亮温图像得到的均衡地表温度为：Tbalancef＝(Tfl×3/5)+(Nf×2/5)；
待提取亮温图像与背景场图像对应的均衡地表温度差值ΔT＝Tbalancelst－Tbalancef。
4.根据权利要求1所述的一种空地数据互增强的震前热异常自动分级方法，其特征在于，所述通过计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值，获取热红外异常提取图像包括以下步骤：
计算待提取亮温图像每一个像元值与背景场图像对应的每一个像元值的差分绝对值ΔN；
当ΔT大于阈值时，异常提取值为N＝0；
当ΔT小于阈值时，异常提取值为N＝((阈值-ΔT)/阈值)×ΔN；
遍历提取亮温图像与背景场图像的所有像元，将每个像元对应的异常提取值组成矩阵构成热红外异常提取图像。
5.根据权利要求1所述的一种空地数据互增强的震前热异常自动分级方法，其特征在于，所述处理热红外异常提取图像，生成震前热红外异常优化图像包括以下步骤：
步骤2.1、获取带有地理位置信息的矢量断裂带数据和对应的地震震源信息，得到单位面积下地震次数PER；
步骤2.2、所有断裂带都标上相应的PERj，给每个断裂带Areaj赋相应的权重分级标签Hqj，权重分级标签分为多个等级；j＝1…S，S为矢量断裂带数据个数；
步骤2.3、将每个断裂带的辐射域分为若干层级，每一层级辐射域都分配对应的权重分级标签Hrji；
步骤2.4，对于获取的热红外异常提取图像，逐像元计算震前热红外异常优化值Nconstruction：
如果热红外异常提取图像的像元位于断裂带内，则Nconstruction＝N×Hq；
如果热红外异常提取图像的像元位于辐射域内，则Nconstruction＝N×Hq×Hr；
其他，Nconstruction＝0；
其中，Hq为断裂带权重标签，Hr为辐射域权重标签，N为异常提取值；
步骤2.5，返回步骤2.2直至遍历所有红外异常提取图像，得到震前和热红外异常优化图像；统计震前热红外异常优化图像非零值，求解图像的期望值μ和标准差σ，得到震前热红外异常优化图像非零值正态分布函数；根据期望值μ和标准差σ将震前热红外异常分为若干个级别L，0值归类为正常，其他采用正态分布标记等级。
6.根据权利要求5所述的一种空地数据互增强的震前热异常自动分级方法，其特征在于，所述其他采用正态分布标记等级具体如下：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811129842.1A CN109409230B (zh) | 2018-09-27 | 2018-09-27 | 一种空地数据互增强的震前热异常自动分级方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811129842.1A CN109409230B (zh) | 2018-09-27 | 2018-09-27 | 一种空地数据互增强的震前热异常自动分级方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109409230A true CN109409230A (zh) | 2019-03-01 |
| CN109409230B CN109409230B (zh) | 2021-11-09 |
Family
ID=65465435
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811129842.1A Active CN109409230B (zh) | 2018-09-27 | 2018-09-27 | 一种空地数据互增强的震前热异常自动分级方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109409230B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112665730A (zh) * | 2020-12-30 | 2021-04-16 | 中国地震局地质研究所 | 震前温度异常检测方法、装置、设备和存储介质 |
| CN113834572A (zh) * | 2021-08-26 | 2021-12-24 | 电子科技大学 | 一种无人机非制冷型热像仪测温结果漂移去除方法 |
| CN116701847A (zh) * | 2023-08-07 | 2023-09-05 | 中国石油大学(华东) | 一种基于时空联合的震前热异常提取方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104199084A (zh) * | 2014-08-29 | 2014-12-10 | 中国石油大学(华东) | 一种基于热红外异常信号与小波神经网络的地震预测方法 |
| WO2015116873A2 (en) * | 2014-02-02 | 2015-08-06 | Ertha Space Technologies | Earthquake forecaster |
-
2018
- 2018-09-27 CN CN201811129842.1A patent/CN109409230B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015116873A2 (en) * | 2014-02-02 | 2015-08-06 | Ertha Space Technologies | Earthquake forecaster |
| CN104199084A (zh) * | 2014-08-29 | 2014-12-10 | 中国石油大学(华东) | 一种基于热红外异常信号与小波神经网络的地震预测方法 |
Non-Patent Citations (2)
| Title |
|---|
| YUNTAO MA 等: "《A new method to extract and analyze abnormal phenomenon of earthquake from remote sensing information》", 《2011 IEEE INTERNATIONAL GEOSCIENCE AND REMOTE SENSING SYMPOSIUM》 * |
| 吴立新 等: "《断裂活动及孕震过程遥感热异常分析的研究进展》", 《测绘学报》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112665730A (zh) * | 2020-12-30 | 2021-04-16 | 中国地震局地质研究所 | 震前温度异常检测方法、装置、设备和存储介质 |
| CN112665730B (zh) * | 2020-12-30 | 2023-11-03 | 中国地震局地质研究所 | 震前温度异常检测方法、装置、设备和存储介质 |
| CN113834572A (zh) * | 2021-08-26 | 2021-12-24 | 电子科技大学 | 一种无人机非制冷型热像仪测温结果漂移去除方法 |
| CN113834572B (zh) * | 2021-08-26 | 2023-05-12 | 电子科技大学 | 一种无人机非制冷型热像仪测温结果漂移去除方法 |
| CN116701847A (zh) * | 2023-08-07 | 2023-09-05 | 中国石油大学(华东) | 一种基于时空联合的震前热异常提取方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109409230B (zh) | 2021-11-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110929607B (zh) | 一种城市建筑物施工进度的遥感识别方法和系统 | |
| Zhou et al. | Dynamics of land surface temperature in response to land‐use/cover change | |
| Yan et al. | Comparison of pixel‐based and object‐oriented image classification approaches—a case study in a coal fire area, Wuda, Inner Mongolia, China | |
| CN103413151B (zh) | 基于图正则低秩表示维数约简的高光谱图像分类方法 | |
| Wu et al. | An error-bound-regularized sparse coding for spatiotemporal reflectance fusion | |
| CN103728609B (zh) | 星载多光谱红外传感器交叉辐射定标方法 | |
| CN109409230A (zh) | 一种空地数据互增强的震前热异常自动分级方法 | |
| Ai et al. | Phenology-based Spartina alterniflora mapping in coastal wetland of the Yangtze Estuary using time series of GaoFen satellite no. 1 wide field of view imagery | |
| US20180095959A1 (en) | Large scale processor for satellite data | |
| CN113888386B (zh) | 一种高分辨率的时空无缝地表土壤水分估算方法及系统 | |
| CN111222539A (zh) | 基于多源多时相遥感影像优化和扩充监督分类样本的方法 | |
| Zhu et al. | A robust fixed rank kriging method for improving the spatial completeness and accuracy of satellite SST products | |
| Desjardins et al. | A space–time parallel framework for fine-scale visualization of pollen levels across the Eastern United States | |
| Zhang et al. | Development of a high spatiotemporal resolution cloud-type classification approach using Himawari-8 and CloudSat | |
| CN112668613A (zh) | 基于气象预报和机器学习的卫星红外成像效果预测方法 | |
| Lin | Ionospheric total electron content anomalies due to Typhoon Nakri on 29 May 2008: A nonlinear principal component analysis | |
| Zhao et al. | Automatic extraction of yardangs using Landsat 8 and UAV images: A case study in the Qaidam Basin, China | |
| Bektas Balcik et al. | Determination of magnitude and direction of land use/land cover changes in Terkos Water Basin, Istanbul | |
| Zhao et al. | The distance decay of similarity in climate variation and vegetation dynamics | |
| Yin et al. | Improving LAI spatio-temporal continuity using a combination of MODIS and MERSI data | |
| Tian et al. | Estimation model of global ionospheric irregularities: an artificial intelligence approach | |
| Xu | Application of remote sensing image data scene generation method in smart city | |
| Hutchison et al. | Quantitatively assessing cloud cover fraction in numerical weather prediction and climate models | |
| Shawky et al. | Integrated image processing and GIS-based techniques using knowledge-driven approaches to produce potential radioactivity map for the uraniferous granite of Egypt | |
| Meng et al. | Evaluation of temporal compositing algorithms for annual land cover classification using Landsat time series data |
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 |