CN113192007B - 一种多尺度信息融合的地热异常区提取方法 - Google Patents

Abstract

本发明公开了一种多尺度信息融合的地热异常区提取方法，首先使用全约束最小二乘混合像元分解和亚像元空间引力模型提高卫星影像中热红外波段的空间分辨率，使用单窗算法反演得到高分辨率的温度信息。然后从全局和局部的角度综合提取研究区内的温度异常，并结合地质资料和已知的先验知识圈定地热异常范围，并通过扩大研究区域根据多处已知的温泉位置验证了结果的可靠性。最后在已知的温泉位置等一些典型的地热异常区使用高分辨率的无人机数据做进一步的精细提取。本发明实现了从粗到精的地热异常提取，为后期的人工开采提供参考价值，实现资源的充分利用。

Description

一种多尺度信息融合的地热异常区提取方法

技术领域

本发明涉及遥感领域，特别涉及遥感影像多种地热异常检测方法，尤其涉及一种基于卫星影像和无人机影像从粗到精的多尺度地热异常圈定方法。

背景技术

地热是地球赋予人类的廉价能源，区域性地热调查的有效方式是热红外遥感。因此常使用热红外遥感数据提取地热异常。传统的方法是使用物探、化探进行地热探测，但是这种方法费时费力，效率低下。后来由于对遥感技术的研究程度不断深入，又开发了多种利用热红外遥感数据来获取地表温度的算法，因此采用提取出来的温度异常信息和其他相关因素相结合的方法来圈定地热异常。现有的提取温度异常的方法最常用的是全局阈值，虽然这种方法简单快捷，但是它可能会忽略掉一些局部的高温异常而导致漏掉一些地热异常区域。提取温度异常后常结合地质断裂和已知的先验知识等来确定最终的地热区域，因为地热常常是沿断裂线方向分布的，但是它缺少更高分辨率的热红外数据进行验证，并对一些典型的异常进行精细化提取。

发明内容

发明目的：针对现有使用全局阈值提取温度异常的方法存在的不足，本发明提出一种多尺度信息融合的地热异常区提取方法，旨在从全局和局部的角度分析温度变化的影响，结合地质断裂和已知的资料圈定研究区域地热异常，并采用高分辨率的热红外数据对一些典型的异常区域进一步精细化探测，实现由粗到精的地热异常圈定。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种多尺度信息融合的地热异常区提取方法，包括以下步骤：

步骤1：获取研究区域不同时期的遥感影像并进行预处理，所述预处理包括辐射定标、大气校正、几何校正、镶嵌裁剪；

步骤2：使用全约束最小二乘混合像元分解和亚像元空间引力模型提高遥感影像的分辨率，得到亚像元地物分布图；

步骤3：采用基于地表类型权重叠置法获得高分辨率热红外波段的热辐射值，并在此基础上使用单窗算法反演地表温度；

步骤4：分别使用全局阈值、局部分块和高程分区三种方法提取温度异常范围，计算三种方法所得结果的交集作为最终提取的温度异常区；

步骤5：结合研究区域的地质断裂和步骤4所提取的温度异常区圈定地热异常范围；

步骤6：根据研究区域先验知识采用无人机热红外对地热异常区进一步提取。

作为优选，所述步骤2中，采用的全约束最小二乘混合像元分解包含两个约束条件：一个像元内每种端元的丰度图的像元亮度值范围在0～1之间，并且和为1；首先使用最小噪声分离变换MNF实现数据降维，然后计算纯净像元指数PPI，该指数越大代表像元的纯度越高，接着通过设置阈值选择纯净的像元投影到MNF变换主成分空间，最后结合n维可视化工具确定端元类型。

作为优选，所述步骤2中，使用亚像元空间引力模型SPSAM进行亚像元制图，通过计算像元内各个亚像元和其对应的邻域像元间空间引力大小，确定各亚像元的类别；具体如下：

记p_ij为像元P_ab内的一个亚像元，i,j＝1,2,…,S；a＝1,2,…,L_a；b＝1,2,…,L_b；L_a和L_b分别为低分辨率图像的行数和列数；S为放大比例值；则亚像元p_ij受到其邻域像元中第c类分量的引力之和D_c,ij为：

其中d_k为亚像元p_ij几何中心和像元P_k的几何中心的欧氏距离，c＝1,2,…,C,C是类别数目，N_A是领域像元个数，w_k是各个邻域像元的空间相关性权值，F_c(P_k)是第k个领域像元P_k中第c类的分量值；最后根据D_c,ij的大小来确定像元P_ab内属于类c的亚像元，对应引力值最大的D_c,ij亚像元归于c类。

作为优选，所述步骤3中，计算地表温度时，根据不同的地物对热辐射的贡献差异确定权重；从亚像元地物分布图中选取相应类别的样本，选取同种类别聚集的区域作为样本，以其中一种地物作为参考，采用比值法确定各类型的权重，然后使用叠置法进行像元分解；提高分辨率的热红外波段的热辐射值表示如下：

其中w(i,j)为分割的子像元(i,j)所对应的低空间分辨率像元的权重，I为低分辨率影像中热红外波段像元的辐亮度值，I(i,j)为对应高空间分辨率中子像元(i,j)的辐亮度值；S为放大比例值；

然后使用单窗算法反演地表温度信息，公式如下：

C＝ετ，D＝(1-τ)[1+(1-ε)τ]

其中T_S为地表温度，T_sensor为传感器上的亮度温度，T_a为大气平均温度，ε为比辐射率，τ为大气透射率，a，b为常数。

作为优选，所述步骤4中，使用高程分区方法提取温度异常，确定高温异常区和低温异常区的阈值，公式如下：

T_H＝M^T+F′(1-α_H)σ^T

T_L＝M^T-F′(1-α_L)σ^T

其中T_H为高温异常区阈值，T_L为低温异常区阈值；F是概率论中的F-分布，F′(a)表示分布函数值为a时所对应的标准化温度值；M^T和σ^T分别是研究区域地表温度的均值和方差；α_H和α_L分别是研究区域中高、低温异常区所占比例；

使用全局阈值方法提取温度异常时，全局阈值设置为T_th＝M^T+eσ^T；1≤e≤2；

使用局部分块方法提取温度异常时，将整个研究区域分割成大小一致的子块，然后使用全局阈值T_th分别去提取每个子块中的温度异常。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

(1)由于卫星数据中热红外影像的分辨率较低，本发明在圈定地热异常之前首先进行了亚像元制图提高了热红外影像的分辨率，通过单窗算法获得了高分辨率的温度信息，为后期的地热异常提取奠定了良好的基础。

(2)本发明综合考虑了全局和局部两个角度，分别使用全局阈值、局部分块和高程分区方法提取温度异常，再结合研究区域内的地质资料和先验知识圈定地热异常。

(3)由于研究区域的范围较大，卫星数据的分辨率相对较低，并不能精细地确定出地热存在的位置。本发明使用高分辨率的无人机数据进一步精细化探测，实现了从粗到精的地热异常提取，为后期的人工开采提供参考价值，实现资源的充分利用。

附图说明

图1是本发明的方法流程图；

图2是SPSAM方法距离计算示意图；

图3是局部分块示意图；

图4是两个时期最终提取的温度异常图；

图5是圈定的全域地热异常区和重点地热异常区；

图6是地热异常区结果验证图；

图7是无人机精细化提取图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本实施例选用两个时期获取的Landsat-8卫星数据进行实验，研究区域是整个青岛市，多光谱影像的空间分辨率为30米，热红外波段的空间分辨率为100m，获取时间分别是2017年2月11日和2019年8月28日。

本发明所述的多尺度信息融合的地热异常区提取方法，如图1所示，步骤如下：

步骤1：获取影像和影像预处理。为了减少不同时相对结果的影响，首先分别获取研究区域不同时期(青岛市夏季和冬季两个时期)的遥感影像并进行辐射定标、大气校正、几何校正、镶嵌裁剪。

步骤2：亚像元制图。使用全约束最小二乘混合像元分解和亚像元空间引力模型(SPSAM)提高遥感影像的分辨率，得到亚像元地物分布图。

首先结合研究区域内的具体情况，将区域内的地物分为四类：植被、低反照度地物、高反照度地物和水体，然后使用全约束最小二乘方法进行混合像元分解得到各类别的丰度图像，其误差平均值为0.158268，标准差为0.237763。采用的全约束最小二乘混合像元分解包含两个约束条件：一个像元内每种端元的丰度图的像元亮度值范围在0～1之间，并且和为1；首先使用最小噪声分离变换(MNF)实现数据降维，然后计算纯净像元指数(PPI)，该指数越大代表像元的纯度越高，接着通过设置阈值选择较为纯净的像元投影到MNF变换主成分空间，最后结合n维可视化工具确定端元类型。

使用亚像元空间引力模型(SPSAM)进行亚像元制图，在SPSAM中相关性通过空间引力来描述，通过计算像元内各个亚像元和其对应的邻域像元间空间引力大小，确定各亚像元的类别；具体如下：

记p_ij为像元P_ab内的一个亚像元，i,j＝1,2,…,S；a＝1,2,…,L_a；b＝1,2,…,L_b；L_a和L_b分别为低分辨率图像的行数和列数；S为放大比例值；则亚像元p_ij受到其邻域像元中第c类分量的引力之和D_c,ij为：

其中d_k为亚像元p_ij几何中心和像元P_k的几何中心的欧氏距离，c＝1,2,…,C,C是类别数目，N_A是领域像元个数，w_k是各个邻域像元的空间相关性权值，F_c(P_k)是第k个领域像元P_k中第c类的分量值；最后根据D_c,ij的大小来确定像元P_ab内属于类c的亚像元，对应引力值最大的D_c,ij亚像元归于c类。

本实施例中，放大比例S为4时(即每个像元被分割为4*4即16个亚像元)，可以将100m空间分辨率热红外影像提高到了25m；SPSAM方法的距离计算示意图如图2所示。

步骤3：地表温度反演。采用基于地表类型权重叠置法获得高分辨率热红外波段的热辐射值，并在此基础上使用单窗算法反演地表温度。

计算地表温度时，首先由于不同的地物在温度表现方面存在很大的差异，因此根据不同的地物对热辐射的贡献差异确定权重；从亚像元地物分布图中选取相应类别的一些样本，尽量选取同种类别较为聚集的区域作为样本，选取的样本数量和对应辐亮度的均值和方差如表1所示。以其中一种地物作为参考，通常以较好识别的水体为参考，采用比值法确定各类型的权重，例如植被的权重为9.753999/9.272966＝1.051875。同样的方法依次计算出各类型的权重如表2所示。

表1不同地物的辐亮度值

地物	样本数	平均值	标准差
				水体	1713	9.272966	0.021264
植被	717	9.753999	0.049055
				第反照度地物	721	10.055379	0.070005
高反照度地物	762	10.436804	0.097645

表2各类型的权重(相对于水体)

地物	水体	植被	低反照度地物	高反照度地物
					权重	1.00	1.051875	1.084376	1.125509

使用叠置法进行像元分解；提高分辨率的热红外波段的热辐射值表示如下：

其中w(i,j)为分割的子像元(i,j)所对应的低空间分辨率像元的权重，I为低分辨率影像中热红外波段像元的辐亮度值，I(i,j)为对应高空间分辨率中子像元(i,j)的辐亮度值；S为放大比例值。

使用单窗算法反演地表温度信息，公式如下：

C＝ετ，D＝(1-τ)[1+(1-ε)τ]

其中T_S为地表温度，T_sensor为传感器上的亮度温度，T_a为大气平均温度，ε为比辐射率，τ为大气透射率，a，b为常数，0～70℃时，a＝-67.355351，b＝0.45860。

步骤4：提取温度异常。分别使用全局阈值、局部分块和高程分区三种方法提取温度异常范围，计算三种方法所得结果的交集作为最终提取的温度异常区。

(1)高程分区。对于一些地表起伏较大的区域，提取地热异常时就需要考虑高程对温度变化的影响。通常地表温度的分布基本呈正态分布，对于正态分布来说，选取的阈值T_H表明研究区中地表温度值T高于T_H的区域为高温异常区，而选取的阈值T_L表明地表温度T低于T_L的区域为低温异常区。

式(1)表示确定阈值T_H可得到其对应的高温异常区所占区域比例P为α_H，同样，式(2)表示确定阈值T_L可得到其对应的低温异常区所占比例P为α_L。反之，若假定已知研究区域中高、低温异常区所占比例分别是α_H和α_L，即可确定对应的阈值T_H和T_L。因此，可以由式(1)和式(2)分别推导得到式(3)和式(4)。若用F′(a)表示分布函数值为a时所对应的标准化温度值，通过变换式(3)和式(4)，最终得到用来确定高温异常区和低温异常区的阈值，分别为式(5)和式(6)。当α_H为0.05和0.10时，对应的F′(1-α_H)分别为1.645和1.282。本实施例取α_H为0.10，即以温度信息的“平均值+1.282*标准差”作为阈值进行高温异常提取。

P{T≥T_H}＝α_H (1)

P{T≤T_L}＝α_L (2)

T_H＝M^T+F′(1-α_H)σ^T (5)

T_L＝M^T-F′(1-α_L)σ^T (6)

其中T_H为高温异常区阈值，T_L为低温异常区阈值；F是概率论中的F-分布，F′(a)表示分布函数值为a时所对应的标准化温度值；M^T和σ^T分别是研究区域地表温度的均值和方差；α_H和α_L分别是研究区域中高、低温异常区所占比例。

(2)全局阈值。使用全局阈值方法提取温度异常时，全局阈值设置为T_th＝M^T+eσ^T；1≤e≤2；通常，取反演后温度的“均值+标准差”或“均值+2*标准差”作为全局阈值区提取温度异常，但是由于最后要综合三种方法的结果即计算三种方法的交集，为了保持结果的统一性，使三种方法所取的阈值保持一致，由高程分区提取法中可知地表温度的分布图通常呈正态分布，取高温异常区所占的比例为0.10时，所对应的高温阈值的系数为1.282，即设置阈值为“均值+1.282*标准差”作为全局阈值初次提取温度异常区域。

(3)局部分块。使用局部分块方法提取温度异常时，由于全局阈值可能会忽略掉在一个大范围的研究区域中局部的一些温度异常，因此使用分块的思想，将整个研究区域分割成大小一致的子块，如图3所示，然后使用全局阈值T_th作为阈值分别去提取每个子块中的温度异常。

本实施例在局部分块中，已知研究区域的大小为5742x4488，将整个区域分为了9行8列共72块小区域，每块中包含638x561个像素。在高程分区中，已知研究区域内的高程变化是从-13.1427m～1049.96m，以100m为间隔，整个研究区域被分为了12个高程子区(≤0,(0,100],(100,200],…,(900,1000],>1000)。图4是计算三种方法所得结果的交集获得的青岛市两个时期最终提取的温度异常，图4(a)是2019.8.28的提取结果，图4(b)是2017.2.11的提取结果。

步骤5：圈定地热异常并进行验证。结合研究区域的地质断裂和步骤4所提取的温度异常区圈定地热异常范围；图5是结合青岛市的地质断裂和已知资料圈定地热范围11处。然后对结果进行验证分析，由于青岛市已知的温泉个数只有一个，但是在青岛市附近的烟台等区域已知的温泉数量较多，所以通过扩大范围将烟台市的部分区域包括在内，同样使用步骤4的三种方法提取温度异常，图6显示了结合地质断裂图可知已知的四处温泉位置处都存在延断裂方向分布或在断裂交叉处分布的异常，图6(a)为已知的温度位置示意图，图6(b)～6(e)分别为蓬莱温石汤、栖霞艾山汤、招远汤东泉、牟平于家汤的结果验证图。因此，本实施例基于步骤4的三种方法并结合青岛市的地质断裂圈定的地热范围是相对可靠的。

步骤6：无人机热红外精细化提取。根据研究区域中已知的温泉位置和技术人员的先验知识，对一些典型的或潜在的地热异常区采用高分辨率的无人机热红外数据进一步精细化提取。本实施例中，使用大疆无人机搭载XT2热红外相机分别获得了青岛市东温泉村和西石桥村的高分辨率的热红外数据，图7黑色框中是根据温度异常分布的聚集性和已知的先验知识所圈定的地热异常的区域。

以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围，因此，凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种多尺度信息融合的地热异常区提取方法，其特征在于：包括以下步骤：

步骤1：获取研究区域不同时期的遥感影像并进行预处理，所述预处理包括辐射定标、大气校正、几何校正、镶嵌裁剪；

步骤2：使用全约束最小二乘混合像元分解和亚像元空间引力模型提高遥感影像的分辨率，得到亚像元地物分布图；

步骤3：采用基于地表类型权重叠置法获得高分辨率热红外波段的热辐射值，并在此基础上使用单窗算法反演地表温度；

步骤4：分别使用全局阈值、局部分块和高程分区三种方法提取温度异常范围，计算三种方法所得结果的交集作为最终提取的温度异常区；

其中，使用高程分区方法提取温度异常，确定高温异常区和低温异常区的阈值，公式如下：

T_H＝M^T+F′(1-α_H)σ^T

T_L＝M^T-F′(1-α_L)σ^T

其中T_H为高温异常区阈值，T_L为低温异常区阈值；F是概率论中的F-分布，F′(a)表示分布函数值为a时所对应的标准化温度值；M^T和σ^T分别是研究区域地表温度的均值和方差；α_H和α_L分别是研究区域中高、低温异常区所占比例；

使用全局阈值方法提取温度异常时，全局阈值设置为T_th＝M^T+eσ^T；1≤e≤2；

使用局部分块方法提取温度异常时，将整个研究区域分割成大小一致的子块，然后使用全局阈值T_th分别去提取每个子块中的温度异常；

步骤5：结合研究区域的地质断裂和步骤4所提取的温度异常区圈定地热异常范围；

步骤6：根据研究区域先验知识采用无人机热红外对地热异常区进一步提取。

2.根据权利要求1所述的多尺度信息融合的地热异常区提取方法，其特征在于：所述步骤2中，采用的全约束最小二乘混合像元分解包含两个约束条件：一个像元内每种端元的丰度图的像元亮度值范围在0～1之间，并且和为1；首先使用最小噪声分离变换MNF实现数据降维，然后计算纯净像元指数PPI，该指数越大代表像元的纯度越高，接着通过设置阈值选择纯净的像元投影到MNF变换主成分空间，最后结合n维可视化工具确定端元类型。

3.根据权利要求1所述的多尺度信息融合的地热异常区提取方法，其特征在于：所述步骤2中，使用亚像元空间引力模型SPSAM进行亚像元制图，通过计算像元内各个亚像元和其对应的邻域像元间空间引力大小，确定各亚像元的类别；具体如下：

记p_ij为像元P_ab内的一个亚像元，i,j＝1,2,…,S；a＝1,2,…,L_a；b＝1,2,…,L_b；L_a和L_b分别为低分辨率图像的行数和列数；S为放大比例值；则亚像元p_ij受到其邻域像元中第c类分量的引力之和D_c,ij为：

其中d_k为亚像元p_ij几何中心和像元P_k的几何中心的欧氏距离，c＝1,2,…,C,C是类别数目，N_A是领域像元个数，w_k是各个邻域像元的空间相关性权值，F_c(P_k)是第k个领域像元P_k中第c类的分量值；最后根据D_c,ij的大小来确定像元P_ab内属于类c的亚像元，对应引力值最大的D_c,ij亚像元归于c类。

4.根据权利要求1所述的多尺度信息融合的地热异常区提取方法，其特征在于：所述步骤3中，计算地表温度时，根据不同的地物对热辐射的贡献差异确定权重；从亚像元地物分布图中选取相应类别的样本，选取同种类别聚集的区域作为样本，以其中一种地物作为参考，采用比值法确定各类型的权重，然后使用叠置法进行像元分解；提高分辨率的热红外波段的热辐射值表示如下：

其中w(i,j)为分割的子像元(i,j)所对应的低空间分辨率像元的权重，I为低分辨率影像中热红外波段像元的辐亮度值，I(i,j)为对应高空间分辨率中子像元(i,j)的辐亮度值；S为放大比例值；

然后使用单窗算法反演地表温度信息，公式如下：

C＝ετ，D＝(1-τ)[1+(1-ε)τ]

其中T_S为地表温度，T_sensor为传感器上的亮度温度，T_a为大气平均温度，ε为比辐射率，τ为大气透射率，a，b为常数。

Images