CN114707385A - 基于Krylov子空间的深部地层导热系数三维预测方法及装置 - Google Patents

Abstract

本发明提供了一种基于Krylov子空间的深部地层导热系数三维预测方法及装置，包括以下步骤：在均匀半空间研究区域内构建导热系数异常体，设定该研究区域边界条件并开展有限元温度数值模拟，获取地下空间三维温度场d_obs；构造初始预测模型和正则化目标函数，在预测过程中采用Jacobian‑freeKrylov子空间技术求解Jacobian矩阵与任意向量的乘积避免大型稠密的Jacobian矩阵的求解及存储；利用Gauss‑Newton算法及L‑BFGS算法分别构造Hessian矩阵并近似求解Hessian矩阵的逆矩阵，用于减少存储需求与计算量，并获取模型修正量Δm，基于Wolfe准则搜索模型步长更新模型参数，循环预测使得实测数据与模拟数据拟合差小于预设值，输出最优预测结果。本发明能够量化表征深部介质导热系数分布特征，预测精度高范围广，实用性强。

Description

基于Krylov子空间的深部地层导热系数三维预测方法及装置

技术领域

本发明涉及导热系数预测领域，尤其涉及一种基于Krylov子空间的深部地层导热系数三维预测方法及装置。

背景技术

地热资源作为一种稳定可持续的可再生清洁能源，目前正受到全社会前所未有的高度重视，开发利用地热资源对于助力“双碳”目标具有重要的现实意义及深远影响。深部岩石的热物性如导热系数、比热容以及放射性生热率等直接影响到地球内部各个圈层间岩石的热生成、热储存和热传递，是地表、地球内部温度分布及热传递研究不可缺少的参数。岩石热物性也是定量评价温度随时间变化的先决条件，因为它们直接影响地下传热过程或温度变化。此外，在地源热泵设计以及地热开发工程中，也需要利用热物性参数来定量评价地源热泵系统及换热工质等与围岩地质体之间的热交换。

导热系数作为岩石矿物固有的物理属性，是各种热物性参数中最重要的参数之一。不同地层中，暴露出不同岩石矿物学类型，岩石导热系数在横向和纵向上都会发生明显变化，从而改变局部和区域的热构造演化、热结构及地热地质特征。目前获取地下岩石导热系数的方法大致可以分为两大类：(1)通过实验室对钻屑或岩心样品测量来确定；(2)通过地球物理测井来刻画。但这些方法仅能刻画点尺度或者线尺度的岩石导热系数，无法预测区域性导热系数展布特征。然而目前对于深部地层导热系数的三维预测尚处于空白阶段，直接制约研究区地热资源热成因机制解译及高温地热靶区选址。

发明内容

针对深部地层导热系数三维预测领域的空缺状态，本发明提出了一种基于Krylov子空间的深部地层导热系数三维预测方法及装置。以地下空间温度场作为观测数据，地下热流值及各地层热物性参数等作为先验信息，开展有限元温度数值模拟并构造正则化目标函数，基于Gauss-Newton算法构造Hessian矩阵和梯度向量，利用Jacobian-free Krylov子空间技术求解Jacobian矩阵与任意向量的乘积避免大型稠密的Jacobian矩阵的求解及存储，L-BFGS算法近似求解Hessian矩阵的逆矩阵，用于减少计算量并获取模型修正量，Wolfe准则搜索模型步长更新模型参数，循环预测使得实测数据与模拟数据拟合差小于预设值，输出模型即为最优预测结果，实现深部介质导热系数的三维预测。

根据本发明的第一方面，本发明提供了一种基于Krylov子空间的深部地层导热系数三维预测方法，包括以下步骤：

在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向；

根据所述地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；

根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；

将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。

优选地，所述边界条件包括：

上边界温度为恒定温度T_up，满足：

T＝T_up

四周边界为绝热边界，满足绝热边界条件：

下边界为热流值q_down，满足热流值条件：

其中，T为温度，单位℃，n为边界法线方向，无量纲，k为导热系数，单位W/m·K，q_down为热流值，单位W/m²。

优选地，所述有限元温度数值模拟，采用六面体结构离散方式，六面体单元内温度T、导热系数k均采用线性插值：

其中，i为节点标号，T_i和k_i分别为单元内各节点的温度和导热系数，单位分别为℃和W/m·K；N_i为形函数，满足：

其中，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲；母单元中ξ、η、ζ与子单元中坐标x、y、z关系公式：

其中，x₀、y₀、z₀分别为子单元在x、y、z轴方向上的中点；a、b、c分别为子单元在x、y、z方向上的边长，单位为m。

优选地，所述结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场的步骤，包括：

给定地下空间温度场的微分方程为：

对所述地下空间温度场的微分方程乘以冲激函数δT并积分：

根据场论中哈密顿算子运算规则及奥高公式，并代入边界条件，得到地下空间温度场的积分方程为：

其中，Ω为研究区域，Γ_up为上边界，Γ_down为下边界，q_down为热流值，k表示导热系数，单位W/m·K，T为温度，单位℃，

为散度，

为梯度，δT为冲激函数；

对所述地下空间温度场的积分方程进行离散剖分、线性插值，求解得到积分方程内各项的单元积分：

其中，δT_e为在单元内的冲激函数数组，T_e为在单元内的温度数组，K_1e、K_2e为在单元内的刚度矩阵，P_1e、P_2e为在单元内的源向量，求解得到刚度矩阵和源向量内系数：

其中，k_1ij、k_2ij、p_1ij、p_2ij分别为刚度矩阵K_1e、K_2e、源向量P_1e、P_2e内的元素，i、j、l为节点标号，a、b、c分别为子单元在x、y、z方向上的边长，单位为m，k表示导热系数，单位W/m·K，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲，T_up为上边界温度，单位℃，q_down为热流值，单位W/m²；

根据每个单元各个节点的编号并扩展至与研究区域相对应的元素位置处，即得合成后的总体刚度矩阵和源向量的线性方程组：

∑_Ω[δT^T(K₁+K₂)T-δT^T(P₁+P₂)]＝0

其中，δT为研究区域冲激函数数组，T为研究区域内的温度数组，K₁、K₂为研究区域内的刚度矩阵，P₁、P₂为研究区域内的源向量，略去δT^T，化简后有：

KT＝P

其中，K＝K₁+K₂，P＝P₁+P₂；

经过求解该线性方程组KT＝P，实现有限元数值模拟后，得到地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向。

优选地，所述正则化目标函数Φ的表达式为：

Φ＝||D[d(m)-d_obs]||²+λ||W(m-m_apr)||²

其中，D为数据加权矩阵；d_obs为观测数据，即地下空间温度场T(x,y,z)；m为当前预测模型参数向量；d(m)为在当前预测模型参数下的模拟数据；W为光滑约束矩阵；λ为阻尼因子；m_apr为预测模型先验信息；

基于有限元中总体刚度矩阵和源向量的线性方程组KT＝P，由于刚度矩阵K中含有导热系数k，而源向量P与导热系数k无关，因此对该线性方程组中导热系数k求导得到：

进行转换，得到转换公式：

为保证温度值及导热系数非负，预测时设定数据参数和模型参数为温度的对数和导热系数的对数：d＝lnT，m＝lnk；根据链式求导法则，Jacobian矩阵中元素为：

其中，i＝1,2,…,N，N为数据参数个数，j＝1,2,…,M，M为模型参数个数，d_i为研究区域内第i个节点的数据参数，m_j表示研究区域内第j个节点的模型参数，T_i表示研究区域内第i个节点的温度，k_j表示研究区域内第j个节点的导热系数，则Jacobian矩阵写为：

其中，Q和

为对角矩阵，

G为元素G(i,j)组成的矩阵；

对于由无量纲元素组成的任意向量v＝(v₁,v₂,…,v_M)^T，并基于有限元刚度矩阵所得的转换公式有：

其中，e₁,e₂,…,e_M均为单位向量，K(v)为导热系数取v时的刚度矩阵；因此，Jacobian矩阵J与任意向量v的乘积为：

同理，计算得到Jacobian矩阵的转置与由无量纲元素组成的任意向量y＝(y₁,y₂,…,y_N)^T的乘积。

根据Gauss-Newton算法，忽略二阶信息项，得到如下线性方程组求解模型修正量Δm：

HΔm＝-g

其中，H为Hessian矩阵，H＝(J^TD^TDJ+λW^TW)；g为梯度向量，g＝[J^TD^TD(d-d_obs)+λW^TW(m-m_apr)]。

优选地，所述获取所述Hessian矩阵的逆矩阵的步骤，包括：

采用L-BFGS算法近似求得所述Hessian矩阵的逆矩阵。

优选地，所述采用L-BFGS算法近似求得所述Hessian矩阵的逆矩阵时，只利用最近n次迭代的信息：

其中，

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k，

I为单位矩阵，k-n,k-n+1,…,k+1为迭代次数，

为Hessian矩阵H在第k+1次的近似逆矩阵，m_k、m_k+1分别为第k次和第k+1次迭代时的模型参数，g_k、g_k+1分别为第k次和第k+1次迭代时的梯度向量。

优选地，所述根据所述预测模型修正量，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数，需满足如下不等式：

其中，α为搜索步长，无量纲，c为无量纲的数；在求取搜索步长过程中，令α初值为1，若不满足该不等式，则令

α＝ρα

其中，ρ为大于0且小于1的小数，循环直至满足Wolfe准则；因此，模型参数的更新公式为：

m_k+1＝m_k+αΔm_k

其中，m_k、m_k+1分别为第k次和第k+1次迭代时的模型参数，Δm_k为第k次迭代所得的模型修正量。

优选地，所述数据拟合差的计算公式如下式：

其中，RMS为均方根，用来衡量实测值与预测值之间的偏差，值越小表明预测结果越可靠；d(m)ⁱ、d_obs ⁱ分别为第i个预测值和实测值，n为观测数据个数，无量纲。

若数据拟合差小于预设值，则停止迭代并输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量。

根据本发明的另一方面，本发明提供了一种基于Krylov子空间的深部地层导热系数三维预测装置，包括以下模块：

有限元温度数值模拟模块，用于在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向；

模型修正量计算模块，用于根据所述地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；

模型参数更新模块，用于根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；

导热系数预测模块，用于将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。

本发明提供的技术方案具有以下有益效果：

本发明提出了基于Krylov子空间的深部地层导热系数三维预测方法及装置，采用Jacobian-free Krylov子空间技术减少在预测过程中所需的计算内存及计算量，实现深部介质三维导热系数准确快速预测，实用性强、应用范围广，结果准确率达到了85.4％。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明的一种基于Krylov子空间的深部地层导热系数三维预测方法的流程图；

图2为本发明的真实导热系数模型分布图，其中图2(a)为XZ方向截面(y＝2.5km)导热系数分布图，图2(b)为YZ方向截面(x＝10km)导热系数分布图；

图3为本发明的研究区域地下空间(y＝2.5km)温度场特征图；

图4为本发明的在先验信息下(y＝2.5km)的温度场特征图；

图5为本发明的在最优预测结果下(y＝2.5km)的温度场特征图；

图6为本发明的研究区域地下最优导热系数分布图，其中图6(a)为XZ方向截面(y＝2.5km)导热系数分布图，图6(b)为YZ方向截面(x＝10km)导热系数分布图；

图7为本发明的数据拟合差随迭代次数变化曲线图；

图8为本发明的一种基于Krylov子空间的深部地层导热系数三维预测装置的结构图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，图1为本发明的一种基于Krylov子空间的深部地层导热系数三维预测方法的流程图；该方法包括以下步骤：

S1：在均匀半空间研究区域内构建导热系数异常体，对该均匀半空间研究区域设定边界条件，上边界为恒定温度T_up，四周边界为绝热边界，下边界为热流值q_down；结合导热系数分布及边界条件开展有限元温度数值模拟，获取地下空间x,y,z方向温度场T(x,y,z)＝d_obs(参考图3)；

在本实施例中，步骤S1的具体为：在均匀半空间研究区域内构建导热系数异常体(参考图2)，上边界温度为恒定值T_up＝15.0℃：

T＝T_up

其中，T为温度，℃。

四周绝热边界条件：

其中，n为边界法线方向，无量纲。

下边界热流值条件：

其中，k为导热系数数据，W/m·K；q_down为热流值，W/m²。

在上述边界条件下开展有限元温度数值模拟，获取地下空间x,y,z方向温度场T(x,y,z)＝d_obs；

作为可选地实施方式，有限元温度数值模拟，采用六面体结构离散方式，六面体单元内温度T、导热系数k均采用线性插值：

其中，i为节点标号，T_i和k_i分别为单元内各节点的温度和导热系数，单位分别为℃和W/m·K；N_i为形函数，满足：

其中，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲；母单元中ξ、η、ζ与子单元中坐标x、y、z关系公式：

其中，x₀、y₀、z₀分别为子单元在x、y、z轴方向上的中点；a、b、c分别为子单元在x、y、z方向上的边长，单位为m。

给定地下空间温度场的微分方程为：

对所述地下空间温度场的微分方程乘以冲激函数δT并积分：

根据场论中哈密顿算子运算规则及奥高公式，并代入边界条件，得到地下空间温度场的积分方程为：

其中，Ω为研究区域，Γ_up为上边界，Γ_down为下边界，q_down为热流值，k表示导热系数，单位W/m·K，T为温度，单位℃，

为散度，

为梯度，δT为冲激函数；

对所述地下空间温度场的积分方程进行离散剖分、线性插值，求解得到积分方程内各项的单元积分：

其中，δT_e为在单元内的冲激函数数组，T_e为在单元内的温度数组，K_1e、K_2e为在单元内的刚度矩阵，P_1e、P_2e为在单元内的源向量，求解得到刚度矩阵和源向量内系数：

其中，k_1ij、k_2ij、p_1ij、p_2ij分别为刚度矩阵K_1e、K_2e、源向量P_1e、P_2e内的元素，i、j、l为节点标号，a、b、c分别为子单元在x、y、z方向上的边长，单位为m，k表示导热系数，单位W/m·K，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲，T_up为上边界温度，单位℃，q_down为热流值，单位W/m²；

根据每个单元各个节点的编号并扩展至与研究区域相对应的元素位置处，即得合成后的总体刚度矩阵和源向量的线性方程组：

∑_Ω[δT^T(K₁+K₂)T-δT^T(P₁+P₂)]＝0

其中，δT为研究区域冲激函数数组，T为研究区域内的温度数组，K₁、K₂为研究区域内的刚度矩阵，P₁、P₂为研究区域内的源向量，略去δT^T，化简后有：

KT＝P

其中，K＝K₁+K₂，P＝P₁+P₂。

经过求解该线性方程组，实现有限元数值模拟后，得到地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向。

S2：对步骤S1获取的地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略正则化目标函数二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及其转置与任意向量的乘积，避免对于大型稠密的Jacobian矩阵的直接求解及存储；获取Hessian矩阵的逆矩阵并求解Hessian矩阵及梯度向量组成的线性方程组得到模型修正量Δm；

在本实施例中，步骤S2的具体为：设定研究区域的先验信息为均匀半空间(m_ref＝3.0W/m·K)，代入边界条件进行有限元温度数值模拟，得到先验信息下研究区域的地下空间温度场d(m)(参考图4)；

正则化预测中，其正则化目标函数的表达式为：

Φ＝||D[d(m)-d_obs]||²+λ||W(m-m_apr)||²

其中，D为数据加权矩阵；d_obs为观测数据，即地下空间温度场T(x,y,z)；m为当前预测模型参数向量；d(m)为在当前预测模型参数下的模拟数据；W为光滑约束矩阵；λ为阻尼因子；m_apr为预测模型先验信息；

基于有限元中总体刚度矩阵和源向量的线性方程组KT＝P，由于刚度矩阵K中含有导热系数k，而源向量P与导热系数k无关，因此对该线性方程组中导热系数k求导得到：

进行转换，得到转换公式：

为保证温度值及导热系数非负，预测时设定数据参数和模型参数为温度的对数和导热系数的对数：d＝lnT，m＝lnk；根据链式求导法则，Jacobian矩阵中元素为：

其中，i＝1,2,…,N，N为数据参数个数，j＝1,2,…,M，M为模型参数个数，d_i为研究区域内第i个节点的数据参数，m_j表示研究区域内第j个节点的模型参数，T_i表示研究区域内第i个节点的温度，k_j表示研究区域内第j个节点的导热系数，

则Jacobian矩阵可写为：

其中，Q和

为对角矩阵，

G为元素G(i,j)组成的矩阵

对于由无量纲元素组成的任意向量v＝(v₁,v₂,…,v_M)^T，并基于有限元刚度矩阵所得的转换公式有：

其中，e₁,e₂,…,e_M均为单位向量，K(v)为导热系数取v时的刚度矩阵；因此，Jacobian矩阵J与任意向量v的乘积为：

同理，计算出Jacobian矩阵的转置与由无量纲元素组成的任意向量y＝(y₁,y₂,…,y_N)^T的乘积。

根据Gauss-Newton算法，忽略二阶信息项，得到如下线性方程组求解模型修正量Δm：

HΔm＝-g

其中，H为Hessian矩阵，H＝(J^TD^TDJ+λW^TW)；g为梯度向量，g＝[J^TD^TD(d-d_obs)+λW^TW(m-m_apr)]。

作为可选的实施方式，对步骤S2中的Hessian矩阵，利用L-BFGS算法近似求得Hessian矩阵的逆矩阵，减少存储需求及计算量。

具体包括：近似Hessian矩阵逆的更新时只利用最近n次迭代的信息：

其中，

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k，

I为单位矩阵，k-n,k-n+1,…,k+1为迭代次数，

为Hessian矩阵H在第k+1次的近似逆矩阵，m_k、m_k+1分别为第k次和第k+1次迭代时的模型参数，g_k、g_k+1分别为第k次和第k+1次迭代时的梯度向量。

S3：对步骤S2所求得的模型修正量Δm，利用Wolfe准则获得模型搜索步长以更新模型参数。

在本实施例中，步骤S3的具体为：Wolfe准则即需满足如下不等式：

其中，α为搜索步长，无量纲；c通常为一较小的数，如10^-4，无量纲；在求取搜索步长过程中，令α初值为1，若不满足该不等式，则令

α＝ρα

其中，ρ为0.5，可根据实际情况调整在0-1之间调整，循环直至满足Wolfe准则；因此，模型参数的更新公式为：

m^(k+1)＝m^(k)+αΔm^(k)

S4：对步骤S3更新的模型参数代入正演模型得出当前模型参数下的模拟数据，计算与实测数据之间的数据拟合差，若数据拟合差小于预设值，则停止迭代并输出当前预测模型参数作为最优预测结果，并计算当前模型参数下的模拟数据作为最优温场(参考图5)，否则，返回更新d(m)并求解新的预测模型修正量。

在本实施例中，步骤S5中：数据拟合差计算公式如下式：

其中，RMS为均方根，用来衡量实测值与预测值之间的偏差，值越小表明预测结果越可靠；n为观测数据个数，无量纲。

设置拟合差预设值为0.02，若数据拟合差小于预设值，则停止迭代并输出当前模型参数作为最优预测结果(参考图6)；否则进入下一次循环直至数据拟合差小于预设值。参考图7，图7为本发明的数据拟合差随迭代次数变化曲线图；在本实施例中，计算每次迭代中实测数据和模拟数据间的拟合差。实验结果表明，模拟数据与实测数据间的拟合差由0.892收敛至0.018，拟合度较好。

参考图2，图2为本发明的真实导热系数模型分布图；在本实施例中，定义模型评价标准异常体积重叠率(AVOR)，以评估预测结果的准确率：

其中，V_o为预测异常体与实际异常体相交部分的体积，m³；V_p为预测异常体的体积，m³；V_t为实际异常体体积，m³。实验结果表明，小异常体导热系数范围0.88～1.50W/m·K，小异常体积重叠率达到83.3％；大异常体导热系数范围4.50～5.54W/m·K，大异常体积重叠率达到85.4％。由此表明，导热系数预测结果与真实结果拟合程度高。

参考图8，图8为本发明的一种基于Krylov子空间的深部地层导热系数三维预测装置的结构图，该装置包括以下模块：

有限元温度数值模拟模块1，用于在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向；

模型修正量计算模块2，用于根据所述地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；

模型参数更新模块3，用于根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；

导热系数预测模块4，用于将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (10)

1.一种基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，包括以下步骤：

在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向；

根据所述地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；

根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；

将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。

2.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述边界条件包括：

上边界温度为恒定温度T_up，满足：

T＝T_up

四周边界为绝热边界，满足绝热边界条件：

下边界为热流值q_down，满足热流值条件：

其中，T为温度，单位℃，n为边界法线方向，无量纲，k为导热系数，单位W/m·K，q_down为热流值，单位W/m²。

3.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述有限元温度数值模拟，采用六面体结构离散方式，六面体单元内温度T、导热系数k均采用线性插值：

其中，i为节点标号，T_i和k_i分别为单元内各节点的温度和导热系数，单位分别为℃和W/m·K；N_i为形函数，满足：

其中，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲；母单元中ξ、η、ζ与子单元中坐标x、y、z关系公式：

其中，x₀、y₀、z₀分别为子单元在x、y、z轴方向上的中点；a、b、c分别为子单元在x、y、z方向上的边长，单位为m。

4.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场的步骤，包括：

给定地下空间温度场的微分方程为：

对所述地下空间温度场的微分方程乘以冲激函数δT并积分：

根据场论中哈密顿算子运算规则及奥高公式，并代入边界条件，得到地下空间温度场的积分方程为：

其中，Ω为研究区域，Γ_up为上边界，Γ_down为下边界，q_down为热流值，k表示导热系数，单位W/m·K，T为温度，单位℃，

为散度，

为梯度，δT为冲激函数；

对所述地下空间温度场的积分方程进行离散剖分、线性插值，求解得到积分方程内各项的单元积分：

其中，δT_e为在单元内的冲激函数数组，T_e为在单元内的温度数组，K_1e、K_2e为在单元内的刚度矩阵，P_1e、P_2e为在单元内的源向量，求解得到刚度矩阵和源向量内系数：

其中，k_1ij、k_2ij、p_1ij、p_2ij分别为刚度矩阵K_1e、K_2e、源向量P_1e、P_2e内的元素，i、j、l为节点标号，a、b、c分别为子单元在x、y、z方向上的边长，单位为m，k表示导热系数，单位W/m·K，ξ_i、η_i、ζ_i为节点i在母单元上的坐标，无量纲，T_up为上边界温度，单位℃，q_down为热流值，单位W/m²；

根据每个单元各个节点的编号并扩展至与研究区域相对应的元素位置处，即得合成后的总体刚度矩阵和源向量的线性方程组：

∑_Ω[δT^T(K₁+K₂)T-δT^T(P₁+^P ₂)]＝0

其中，δT为研究区域冲激函数数组，T为研究区域内的温度数组，K₁、K₂为研究区域内的刚度矩阵，P₁、P₂为研究区域内的源向量，略去δT^T，化简后有：

KT＝P

其中，K＝K₁+K₂，P＝P₁+P₂；

经过求解该线性方程组KT＝P，实现有限元数值模拟后，得到地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向。

5.如权利要求4所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述正则化目标函数Φ的表达式为：

Φ＝||D[d(m)-d_obs]||²+λ||W(m-m_apr)||²

其中，D为数据加权矩阵；d_obs为观测数据，即地下空间温度场T(x,y,z)；m为当前预测模型参数向量；d(m)为在当前预测模型参数下的模拟数据；W为光滑约束矩阵；λ为阻尼因子；m_apr为预测模型先验信息；

基于有限元中总体刚度矩阵和源向量的线性方程组KT＝P，由于刚度矩阵K中含有导热系数k，而源向量P与导热系数k无关，因此对该线性方程组中导热系数k求导得到：

进行转换，得到转换公式：

为保证温度值及导热系数非负，预测时设定数据参数和模型参数为温度的对数和导热系数的对数：d＝lnT，m＝lnk；根据链式求导法则，Jacobian矩阵中元素为：

其中，i＝1,2,…,N，N为数据参数个数，j＝1,2,…,M，M为模型参数个数，d_i为研究区域内第i个节点的数据参数，m_j表示研究区域内第j个节点的模型参数，T_i表示研究区域内第i个节点的温度，k_j表示研究区域内第j个节点的导热系数，则Jacobian矩阵写为：

其中，Q和

为对角矩阵，

G为元素G(i,j)组成的矩阵；

对于由无量纲元素组成的任意向量v＝(v₁,v₂,…,v_M)^T，并基于有限元刚度矩阵所得的转换公式有：

其中，e₁,e₂,…,e_M均为单位向量，K(v)为导热系数取v时的刚度矩阵；因此，Jacobian矩阵J与任意向量v的乘积为：

同理，计算得到Jacobian矩阵的转置与由无量纲元素组成的任意向量y＝(y₁,y₂,…,y_N)^T的乘积；

根据Gauss-Newton算法，忽略二阶信息项，得到如下线性方程组求解模型修正量Δm：

HΔm＝-g

其中，H为Hessian矩阵，H＝(J^TD^TDJ+λW^TW)；g为梯度向量，g＝[J^TD^TD(d-d_obs)+λW^TW(m-m_apr)]。

6.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述获取所述Hessian矩阵的逆矩阵的步骤，包括：

采用L-BFGS算法近似求得所述Hessian矩阵的逆矩阵。

7.如权利要求6所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述采用L-BFGS算法近似求得所述Hessian矩阵的逆矩阵时，只利用最近n次迭代的信息：

其中，

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k，

I为单位矩阵，k-n,k-n+1,…,k+1为迭代次数，

为Hessian矩阵H在第k+1次的近似逆矩阵，m_k、m_k+1分别为第k次和第k+1次迭代时的模型参数，g_k、g_k+1分别为第k次和第k+1次迭代时的梯度向量。

8.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述根据所述预测模型修正量，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数，需满足如下不等式：

其中，α为搜索步长，无量纲，c为无量纲的数；在求取搜索步长过程中，令α初值为1，若不满足该不等式，则令

α＝ρα

其中，ρ为大于0且小于1的小数，循环直至满足Wolfe准则；因此，模型参数的更新公式为：

m_k+1＝m_k+αΔm_k

其中，m_k、m_k+1分别为第k次和第k+1次迭代时的模型参数，Δm_k为第k次迭代所得的模型修正量。

9.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述数据拟合差的计算公式如下式：

其中，RMS为均方根，用来衡量实测值与预测值之间的偏差，值越小表明预测结果越可靠；d(m)ⁱ、d_obs ⁱ分别为第i个预测值和实测值，n为观测数据个数，无量纲；

若数据拟合差小于预设值，则停止迭代并输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量。

10.一种基于Krylov子空间的深部地层导热系数三维预测装置，其特征在于，包括以下模块：

有限元温度数值模拟模块，用于在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d_obs，x、y、z分别表示x、y、z轴方向；

模型修正量计算模块，用于根据所述地下空间温度场d_obs，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss-Newton算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian-free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；

模型参数更新模块，用于根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；

导热系数预测模块，用于将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。

Images