一种用于浅层温度分析的反演方法，用于估算加拿大不列颠哥伦比亚省梅格山（Mount Meager）南部地热场中的热异常现象

《Applied Computing and Geosciences》：An inversion method for shallow temperature analysis to estimate thermal anomalies in the southern Mount Meager geothermal field, British Columbia, Canada

　　刘小军 | 斯蒂芬·E·格拉斯比 | 陈卓恒
加拿大自然资源部，加拿大地质调查局，卡尔加里，AB T2L 2A7，加拿大

**摘要**
部署浅层温度探头以记录温度变化是识别高地面热流异常区域的一种方法，随后对地下温度时间序列进行的数据处理和分析对于解释结果至关重要。本研究提出了一种不同的处理工作流程，以更准确和可靠地估计热异常。为了减少日太阳辐射的影响，计算随深度变化的热量，并估计近表面的相对热异常，采用了带有蒙特卡洛（MC）参数生成器的正则化共轭梯度反演算法来解决目标函数，并克服与复杂参数变化相关的不确定性。通过假设热导率和热扩散率随土壤含水量和密度的季节性变化来创建一个物理模型，然后使用经验方程来估计热参数的变化。使用Python语言开发了用于热数据集预处理和分析的工作流程和代码。基于处理后的瞬态温度变化，利用不同深度的长期温度趋势来推断更深层的热梯度。在一项使用真实数据的研究中，所提出的反演和模拟工作流程被应用于计算2.0米深度处的温度变化。结果在梅格山（Mount Meager）南部区域识别出一个相对热流较高的异常区域。通过与之前从浅层温度井获得的结果进行比较，所提出的估计异常热场的方法在识别区域热梯度变化方面具有显著优势，可以成为不同案例研究的实用工具。

**1. 引言**
地热能作为一种清洁且可再生的能源资源具有许多优势，包括相对较高的容量系数、低运营成本和最小的环境影响（Gupta和Roy，2007；Grasby等，2012）。地下温度分布是评估地热资源潜力的重要因素。定量评估热梯度不仅有助于确定热资源量，还有助于定位可能将热液引导到更浅层的水力断裂带（Grasby等，2020；Emekwi等，2024）。自20世纪70年代末以来，一些研究小组（Fairbank，1979；GeothermExInc.，2004；Grasby等，2020）基于地质学、地球物理学、地球化学和热梯度井钻探等综合解释方法，研究了加拿大不列颠哥伦比亚省梅格山地区的地热潜力。为了识别热流路径和渗透带，包括阐明区域和局部水文效应对温度的影响，我们研究了利用地面温度探头记录地下热流的潜力，并开发了相关的数据处理和分析工具（Liu等，2022；Chen等，2023）。分析近表面热流垂直分布的数值方法已经发展了几十年。一些算法已被提出（例如Bennett等，2008），用于估计地表热流，其中总结了地热勘探中常用的热流技术。在高土壤深度研究中，Gao等（2010，2017）提出了一种评估土壤热流（SHF）的新方法，该方法使用量热法以及基于热扩散方程解析解的方法来计算地表以下的热流，而不是直接测量地表热流。尽管浅层地表（<1.0米）的岩石和土壤温度主要受日太阳加热的影响，但这些热数据仍然包含了一些关于深层热流的信息，这也是本研究的重点。使用数值建模可以帮助消除日太阳辐射的影响，并提取仅与深层热源相关的信息。例如，Assouline等（2019）提出了一种傅里叶建模策略来估计表观热扩散率，并将其应用于瑞士的几个地点。Campbell（1985）和Sharratt等（1992）提出了一种使用瞬态热流方程计算热导率的有限差分方法，这比谐波方法更好。

尽管以前的方法能够确定热流，但土壤热属性的不确定性和环境变化可能会使计算结果产生偏差。例如，具有不同土壤组成的土地覆盖会影响地下物理属性，这是估计土壤/岩石热属性的一个影响因素。这是量化随时间变化的地下温度的先决条件。在模型分析中需要考虑热导率和热扩散率的不确定性，特别是对于具有异质岩性的组合模型。如果不考虑这些因素，计算出的热导率/扩散率可能无法代表真实的热流。可以外部管理的影响土壤热容量的因素包括含水量和土壤密度（Xie等，2018），这些因素在土壤热容量中起重要作用，但非常不确定。

此前，我们展示了沿梅格山火山复合体南侧收集的时间序列数据集的初步分析（Liu等，2022；Chen等，2023）。为了减少热流属性的影响并可靠地估计异常热流区域，我们修订了数据分析方法，专注于确定不同深度的近表面温度和热流异常。由于近表面的日变化和季节性温度变化包含了来自上层地热流的信息，观察到的温度变化可以转化为热流的定量估计和热异常的识别。

**2. 方法论**
**2.1. 热属性模型**
热属性影响单位面积单位时间内在单位温度梯度下的热传递量（热流密度）。影响土壤热属性的因素可以归类为组成属性、结构、环境条件等（Chen等，2025），并可以估计创建这些属性与因素之间关系的关键经验系数。当干燥土壤变湿时，热属性会急剧增加。对于这个特定研究区域，影响因素可以简化，例如土壤热导率（λ）主要取决于土壤的体积密度和含水量。对于低频向上热异常研究，我们假设热属性是均匀的，以简化方程。根据Lu等（2014）提出的经验关系（公式（1）：
λ = λdry × exp(β?θ?α)，
其中α和β是与体积密度相关的形状因子，λdry是干燥土壤的热导率。根据这个经验方程，λ（Wm?1K?1）与含水量θ之间的关系，以及土壤体积含水量与公式（1）中相应参数之间的幂关系可以描述为：
λ = 0.2 + exp(1.46?θ?0.34)，
其中0.20 Wm?1K?1是干燥土壤的热导率，1.46和0.34是拟合曲线的形状因子，使用全局优化算法进行拟合（Vrugt等，2003；Gupta和Roy，2007；Fuchs等，2015；Xie等，2018）。利用类似的经验和基于物理的关系，可以基于热导率和热容量模型以及土壤颗粒的体积分数和体积密度来估计热扩散率κs：
κs = (0.69 + exp(3.06?θ?0.26))。

表1列出了一些热属性作为参考。土壤中的水分可以改善土壤颗粒之间的热接触，其热导率比空气低20倍。

**表1. 土壤成分的热属性（Hillel，1980；Brutsaert，1982）**
成分 λ [Wm?1K?1] 比热容 [Jkg?1K?1] 密度 [kgm?3] 体积热容 [MJm?3K?1]
--- --- --- --- ---
土壤矿物（平均） 2.97 33 2650 1.94
土壤有机物 0.25 19 2613 0.25
水（液态） 0.57 41 1821 0.41
冰（0°C） 2.22 10 89 1.93
空气 0.02 5 10 0.12

根据热属性的经验方程（公式（1）、（2）、（3），我们创建了一组表观热扩散率值，代表用于热异常分析的年度变化。由于土壤湿度是一个重要因素，因此将研究区域的月降雨量考虑到了热属性模型中。图1(a)显示了2016年至2021年彭伯顿（Pemberton）的降水量（数据来源：https://climate.weather.gc.ca/climate_data/）。这5年平均值的降水量模式可用于估计土壤含水量。为了简化，忽略了春季融雪的影响。考虑到季节性变化，创建了一个简化的经验模型，通过结合其他效应因素（如土壤矿物组成、颗粒大小、干密度和自然含水量）来模拟热扩散率。而不是使用土壤孔隙度，而是使用土壤密度与矿物密度的比率作为近似因素。使用了一个简化的变热扩散率公式来计算热扩散率（Abu-Hamdeh，2003；Fuchs等，2015；Xie等，2018），并随后用于模拟和数据分析：
κs = 0.25 + exp(b?Sr?a)
其中，a、b、c是曲线的形状参数。我们将这些形状参数固定为0.39、5.25、3.11（Xie等，2018），Sr是根据公式（5）中θ和相对密度rd（ρs，ρm）的比率计算得出的，θ是自然含水量的参数，ρs是土壤密度，ρm是矿物矩阵的密度，本研究中使用了一个典型的值2.67 g/cc（Grasby等，2020）。

**图1.** (a) 彭伯顿的降水量；(b) 土壤含水量与热扩散率的关系；(c) 土壤密度与热扩散率的关系；(d) 使用MC生成器模拟的热扩散率。

图1(b)显示了不同土壤密度下土壤含水量与热扩散率的关系；图1(c)显示了不同含水量下土壤密度与热扩散率的变化。含水量和土壤密度都与热扩散率成正比。由于在记录温度序列的期间没有测量相对湿度，因此只能通过季节性水位/土壤湿度来估计含水量（Duff等，1997）。为了理解不同环境中热属性不确定性的影响，创建了一个带有蒙特卡洛方法的样本生成器来生成含水量和土壤密度的季节性变化（表2）（Lu等，2014；Xie等，2018；Grasby等，2020）。创建了一个基于谐波的模型来模拟变化：
κs(t) = 0.25 + exp(5.25?S(θ,rd)?0.39) × 4.18 × S(θ,rd) + 3.11 × sin(ωt+φ)，
其中ω是径向频率，φ是相位角。图1(d)显示了不同参考参数范围内未冻结土壤的热扩散率概率分布。种子数量设置为100，对于这个简化的热扩散率模型来说是足够的。灰色线条显示了使用蒙特卡洛生成器模拟的热扩散率值。根据热参数的复杂性，未冻结土壤的热扩散率在3到10（10?7m2s?1）之间变化，建模结果被用作计算更深层温度变化和估计异常的输入。

**表2. MC生成器的热扩散率方程参数**
参数 θ ρs α φ
--- --- ------ ---
未冻结土壤 0.35 1.81 π/2 ---
范围 (0.01, 0.9) (0.85, 2.15) (0.8, 1.0) (1e-5, 5e-1)

**2.2. 建模方法**
在本小节中，描述了热流模拟方法。地下能量传递过程的基本理论与之前的出版物类似（Hurley和Wiltshire，1993；Sauer等，2007；Gao等，2010；Xie等，2018）。我们使用谐波方法来模拟测量深度的温度变化，并预测目标深度的地下温度。使用谐波分析模拟热流的详细信息在附录A中有所说明。方程(Eq. (A5))描述了随时间和深度变化的温度，其中包含了使用MC生成器计算的热扩散值。图2展示了在理想条件下48小时内的模拟温度变化。初始表面温度为25°C，热导率为0.9 Wm?1K?1，热扩散率为6(10?7m2s?1)。温度剖面图展示了在不同深度下稳态热流通过介质时的热量传输情况，这表明了太阳热效应的日变化。在较深的区域（>1.0米），这种效应可以忽略不计。由于假设的地下热异常会影响太阳辐射的向下渗透，因此这种用于计算深层温度变化的建模方法可以与季节性温度变化分析结合使用，以估计向下的地下热通量。

下载：下载高分辨率图片（343KB）
下载：下载全尺寸图片
图2. 在理想条件下，由于日太阳效应导致的不同深度的模拟温度变化。

2.3. 反演工作流程
为了减少与复杂地下物理属性相关的不确定性，在模拟中我们假设热属性会随季节变化，并且受到土壤含水量以及土壤密度与矿物密度比例的影响。基于这一假设，通过参数化模型计算热导率和热扩散率的年度变化。我们开发了一个工作流程，结合地球物理反演算法和蒙特卡洛（MC）参数生成器来估计近地表的热分布，然后模拟不同深度的时间序列温度变化。建模结果用于消除太阳辐射的日变化效应，并预测更深层次的温度变化。通过统计方法分析计算出的深层温度序列，以定位可能的热异常。

对于热通量建模，计算方程(Eq. (A5))中的An和Cn的一种简单方法是使用快速傅里叶变换（FFT）（Fourier, 1822; Sauer et al., 2007; Assouline et al., 2019）。FFT的优点是计算速度快，但其输出包括了来自高频噪声的所有频谱，并且不适用于热属性动态变化的模型。为了使用特定的角频率范围，例如日周期和年周期，我们假设研究区域内岩石/土壤的地下表观热扩散率是土壤湿度和密度的函数，并采用经典反演算法来计算这一点。这种处理方法抑制了高频噪声信息，并提高了谐波方程(Eq. (A5))中热扩散率动态变化的灵活性，从而允许估计不同深度的温度。

基于之前的地下热通量计算方法（Gao et al., 2010, 2017），我们创建了一个目标函数Pα(m)，并采用地球物理反演工作流程来最小化失配函数。目标函数中包含了日周期和年周期。然后可以将计算出的幅值和相位用于模拟。我们假设离散的反问题d=Ψ(m)是从模型空间M到数据空间D的连续算子。对于任何d∈D和任何参数α≥0，都存在一个模型mα∈M，使得参数函数为：
(6)Pα(m)=φ(m)+αs(m)=∑i=1n(Ψi(m)?di)2+α∑j=1m(mj?mapr)2=min，其中(Eq. (6))由失配函数φ(m)和稳定函数s(m)组成。我们使用了梯度型技术——正则化共轭梯度（RCG）算法来求解方程(Eq. (6))（Zhdanov, 2002）。详细的迭代步骤在附录A中有所描述。在反演过程中，权重被设置为1，参考模型mapr被设置为0。图3展示了使用迭代反演预测谐波模型、随后进行向下模拟和热异常聚类分析的工作流程。迭代优化的总计算时间约为5分钟，以达到最小失配水平。

下载：下载高分辨率图片（455KB）
下载：下载全尺寸图片
图3. 用于预测谐波模型和热异常分析的RCG优化工作流程。

在获得模型参数后，可以使用方程(Eq. (A5))计算每日温度值。方程(Eq. (7))以等效形式表示，用于计算每个时间间隔i的热流密度，然后可以使用傅里叶第一定律（Fourier, 1822; Sauer et al., 2007）来确定深度处的热流：
(7)Gi=?kiΔTiΔzi, i=1,2?
为了测试所提出的反演程序，收集了合成和公共的温度数据集。原始数据来自不列颠哥伦比亚省彭伯顿的天气数据库（红色虚线）。数据集的时间跨度为两年。然后使用图2中的相同物理参数进行反演和建模计算。图4展示了使用计算出的谐波参数模拟的年度温度分布。结果显示了温度随深度增加的趋势（0.5, 1.0, 2.0米）。根据模型测试结果，所提出的方法可用于实际数据处理和模拟向下的温度分布及其年度变化。

下载：下载高分辨率图片（464KB）
下载：下载全尺寸图片
图4. 彭伯顿（2019-2021年）反演测试结果和模拟的年度向下热变化（日均值）的图表。

3. 实际数据和分析
加里巴尔迪火山带的研究重点在于开发新的方法，以帮助降低对盲目的地热资源的勘探风险（Grasby et al., 2020）。除了地质、地球物理和地球化学研究外，还在梅格尔山火山区的南侧部署了温度探针（Chen et al., 2023）。这些实际数据被用来评估这里提出的方法。背景、温度记录和热通量异常如下所述。

3.1. 背景和先前的工作
梅格尔山位于温哥华以北约170公里处的一片未开发的山区（图5(b)）。它位于加里巴尔迪火山带的北端，该火山带是美国喀斯喀特山脉在加拿大的延伸部分。梅格尔火山复合体是一个相对较年轻的火山，大约在200万年前至几千年前处于活跃状态（图5(a)中的粉红色区域），至今仍存在温泉和喷气孔（Read, 1979; GeothermExInc., 2004）。火山序列包括安山岩熔岩流、火山碎屑单元以及英安岩和流纹英安岩穹丘和流。火山复合体下方的基底由 intrusive quartz diorite、granodiorite 和变质岩组成。在梅格尔山南侧的地质热勘探区钻探的井穿透了较老的侵入岩，以及少量岩脉和热液角砾岩（GeothermExInc., 1995）。岩脉和角砾岩之间的密切关联表明宿主岩石经历了显著的加热和随后的断裂。这种断裂对于热液系统的发展是必要的，因为它在结晶岩石中创造了必要的渗透性，使热水能够从更深、更热的地层（如冷却的岩浆附近）传导到可以通过钻探到达的较浅环境（GeothermExInc., 2004）。

下载：下载高分辨率图片（1MB）
下载：下载全尺寸图片
图5. 不列颠哥伦比亚省梅格尔山南侧研究区域的位置和温度探针（红色三角形）的分布（基础地图来源：http://wms.ess-ws.nrcan.gc.ca/wms）。关于此图中颜色参考的解释，请参阅文章的网络版本。梅格尔山地区近50年来一直是地热勘探的目标。20世纪70年代，加拿大启动了地热能计划，首次了解了加拿大地热状况（Fairbank et al., 1980; B.C. Hydro, 1983; Jessop, 2008）。在梅格尔山南部钻探的勘探井显示了高温地热资源（>250°C）和高热梯度。该地区最近的研究包括结合地质线状构造、地球物理地震及大地电（MT）和岩石样本分析等数据集进行地热资源评估（Grasby et al., 2020; Chen et al., 2023; Hormozzade Ghalati et al., 2023; Liu et al., 2024）。

3.2. 野外测量
本研究中使用的温度探针是由Onset制造的Hobo Water Temp Pro v2型号。该探针的可记录温度范围是从-40到70°C，精度为±0.2°C，分辨率为0.02°C，年漂移量为0.1°C。探针防水，以半小时为间隔记录数据，并配备可支持约六年数据记录的电池。探针的地理位置如图5(a)所示。数据收集时间为2020年9月至2021年8月（Chen et al., 2023）。探针的埋设深度在1至7厘米之间，以避免直接接收太阳辐射，站点海拔高度在550至1400米之间（表3）。

表3. 所有温度探针站点的一般信息。
站点ID A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22
遮荫指数 0.2 0.05 0.1 0.7 0.4 0.5 0.3 0.5 0.6 0.8 0.6 0.5 0.6 0.8 0.6 0.4 0.9 0.95
埋设深度（厘米） 21 25 45 42 53 33 36 24 36 42 37 42 33 0 25 25 25 25 33 30 25 25 71 57 15 70 59 35
海平面高度（米） 138 61 28 31 117 81 110 51 117 81 44 98 29 07 88 28 39 80 97 84 77 67 47 74 74 73 96 93 56 25 71 57 05

在数据预处理期间，开发的代码会扫描所有数据集，确保时间序列记录的连续性，例如检查温度值以统一单位，并应用dropna(.)函数来移除异常点。图6(a)显示了9月份的所有测量数据，图6(b)显示了从12月到3月选定的冬季时段的数据。在遮荫度低的区域部署的探针记录的日温度变化较大。当表面被导热性低的雪覆盖时，会阻挡太阳热通量，导致温度没有日变化（Chen et al., 2023）。为了在同一尺度上分析时间序列数据，选择了夏季或冬季相同时间段的温度数据集进行模拟和统计分析，在RCG迭代处理期间输入所有间隔为30分钟的记录点。

下载：下载高分辨率图片（731KB）
下载：下载全尺寸图片
图6. 经过预处理后，所有站点在此期间记录的所有温度数据，(a) 全年；(b) 选定的冬季时段。

3.3. 反演和解释
利用创建的热扩散率变化模型，采用所提出的RCG迭代反演方法来计算梅格尔山数据集中的100个变量，用于谐波分析。为了避免积雪和雪融化对热属性模型的影响，我们认为只有5月15日至8月15日期间的数据对向下模拟和解释是有意义的。图7显示了选定样本站点和时间段的变化图。蓝色虚线代表记录的温度，其他线条代表在固定热扩散率6(10?7m2s?1)下，不同深度（0.2, 0.5, 1.0米）的模拟温度。观察到太阳热通量的日变化和模拟温度变化。不同深度的计算温度被绘制出来，以显示太阳辐射随深度增加而减少的效果。即使在直接暴露于阳光下的站点，温度变化也随着深度的增加而稳定。1.0米深度的温度显示了最小的日阳光效应（红色线条）。

下载：下载高分辨率图片（1MB）
下载：下载全尺寸图片
图7. 2021年5月15日三个月期间，在不同深度（0.2至1.0米）的计算温度变化，固定热扩散率的情况下，以及在不同热扩散率模型下（2.0米在站点4A(a)、16A(b)和21A(c)）。由于我们假设来自下方的表观热通量是恒定的，而近地表的热扩散率随季节变化且模式不确定，我们使用MC生成器生成的100个热参数，在2.0米深度计算了所有可能的温度分布。图7中的灰色区域显示了2.0米深度可能的温度范围。整个时期的模型温度变化的不确定性小于5°C。与向上热通量相关的预测变化被保存下来，用于后续的统计分析，以检测相对异常和解释。

为了更容易解释温度随深度的变化，计算了四个组成部分指数的统计结果。根据能量平衡方程（Eq. (1)），温度梯度将受到地下热通量的影响。这意味着温度变化包括了向上的热信息，因此可以计算温度变化梯度和其相对异常指数。这四个组成部分指数包括正（向上）和负（向下）方向的梯度值；2.0米处的温差；以及冬季的平均温度。公式表示具有正向和负向梯度分量的标准化相对差异，Gp=1/n?1∑i=2^n(Ti?Ti?1)；2.0米处峰值和最低温度之间的温差，Gm=max(T2)?min(T2)；以及冬季（12月至3月）的平均温度，Gw=mean(T2)。所有100个MC生成的热扩散率值都用于所有地点的计算，模拟的热分布则用于计算不同的组分指数。图8显示了每个地点的四个组分统计量的误差条，这些指数用于识别高热流量的潜在位置。从异常指数分布来看，地点10A、11A、15A到18A以及21A显示出较高的数值范围，因此被识别为热异常区。图9(a)和(b)展示了冬季的温度统计结果，以及夏季的温度梯度统计结果。使用scikit-learn中的无监督机器学习聚类算法——基于密度的带噪声空间聚类（DBSCAN）对地点进行分组，该算法将群集视为不同密度区域（图9(c)）。图9(d)绘制了结合了四个组分指数的标准化加权汇总的指数，我们可以看到潜在增强热流量的高相对差异区域围绕Meager Creek的地点15A和21A。

图8. 基于夏季预测热变化的计算组分指数：(a)表示正梯度；(b)表示负梯度；(c)表示夏季的温差；(d)表示包含所有加权指数的合并特征。
图9. 异常指数的统计分布：(a)基于选定冬季的平均温度；(b)基于正梯度值；(c)通过提取不同组分进行分类；(d)基于图8中四个指数的组合分布，线M-M'是一个示意性剖面（点的大小代表高热异常的相对差异值，等高线表示海拔高度）。

4. 讨论
2021年夏季收集的温度数据显示了每日太阳热流量的影响。因此，近地表处的计算温度不足以独立反映向上的热流量效应。在本研究中，对所有探针地点运用了变化的表观热导率/扩散率进行反演处理，这提供了合理的相对热异常范围。由于我们仅将降水作为水分含量的因素，年度热扩散率变化的不确定性可能会实际受到影响，从而影响预测结果。为了进行综合比较，还需要测量土壤样本的热导率。将模拟结果与MC生成的热特性进行统计分析，只能显示热异常的高概率区域。在模拟过程中可以考虑不同地点的热导率变化，这将缩小综合解释的偏差范围。使用反演算法模拟的时间序列温度变化可以生成一年的热变化趋势，但在Mount Meager地区的每个地点仅有一个温度探针，无法完全消除太阳热辐射的影响，以计算地热流量值。根据当前的数据集和反演结果，我们只能计算温度的相对变化，以显示基于差异指数的深层地热热流的异常分布。结合多年温度序列的综合分析将有助于减少解释的不确定性。在积雪覆盖期间，记录的温度不能准确反映季节变化，因为雪的热导率远低于地表岩石或土壤（图6(b)）。

作为与太阳辐射效应相关的两个重要因素，探针地点的遮阴指数（SI）和土壤湿度（SM）在部署探针时进行了测量。为了评估这两个因素对记录温度的影响，使用了一种简单的统计方法。为简化计算，在解释过程中未考虑地形的影响。图10显示了遮阴指数和土壤湿度与温度指数之间的关系（对于>20°C的情况），我们可以看到随着SI的增加，太阳辐射效应减小，而随着SM水平的增加，太阳辐射效应略有增加。由于其他因素（如地形）的影响，一些地点的统计数据显示出较大的偏差。在未来的研究中，应考虑埋藏深度、遮阴区和地形效应，以减少数据分析和解释的不确定性。

图10. 遮阴指数（a）和土壤湿度（b）与温度指数的交叉图。
基于之前的研究（Fairbank, 1979, Fairbank et al., 1980; Lewis and Souther, 1978），我们更新了一个剖面图（图11），以展示研究区域内的横截面视图（见图9中的M-M'部分）。图11(a)显示了M-M'剖面上不同钻孔的温度记录剖面。相对温差较大的地点（15A、21A）与附近的钻孔温度记录（M1-74D和301-2）一致（Lewis and Souther, 1978）。这两个地点在浅层区域（<50米）显示出较高的平均温度。这一结果验证了Meager Creek沿线的浅层热梯度区以及统计分析的可靠性。尽管其他井（例如M5-78D井）在深层显示出高温，但在近地表区域的温度相对较低。一个可能的原因是上覆英安岩和安山岩的低孔隙度和渗透性。在这些井旁边增加更多的探针站点，并进行多年观测，将有助于提高信噪比，并在未来的研究中更准确地解释异常区域。在当前的模拟框架中，我们专注于通过简化热特性模型和排除流体流动动力学的地下效应来进行方法学研究，这可能限制了输出的应用性。将本研究与其他地球科学数据结合进行综合解释，将提高识别热异常的准确性。

图11. 研究区域内的垂直剖面M-M'：(a) M-M'剖面上附近井的温度记录剖面；(b) M-M'剖面显示了储层温度分布和地质情况（根据Fairbank, 1979; Liu et al., 2022进行修订）。

5. 结论
时间序列温度记录数据的预处理、反演和解释工作流程将有助于识别潜在的热异常区域。MC生成器可以创建动态热特性结果，可用于预测温度变化，并协助进行热异常的统计分析，同时减少数据噪声带来的不确定性。本研究侧重于模拟和反演算法的方法学发展，以减少太阳辐射的影响，并通过温度时间序列的统计分析识别浅层热异常，其中包含了灵活选择角频率分布的选择。尽管采用了RCG迭代反演和MC生成器，计算仍然高效。输出可以提供支持性信息和证据，用于区域地表热分布，这对于地热资源评估和勘探非常有用。尽管仅使用了单深度水平的温度时间序列，但通过与以往研究的比较，仍可以在Mount Meager地区识别出可能的过高热流量的异常区域。这种方法可用于未来研究中的多深度测量的数据反演分析。将其他地球科学数据结合进行综合解释，将减少不确定性并产生更准确的建模结果。

作者贡献声明：
刘晓军：撰写——原始草稿、软件、方法学、形式分析
Stephen E. Grasby：撰写——审核与编辑、验证、监督、项目管理、调查、资金获取
陈卓恒：撰写——审核与编辑、验证、监督、调查、数据管理、概念化

资金声明：
本项工作完全由British Columbia省的Geoscience BC项目和Natural Resources Canada的GeoEnergy计划资助。