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

本发明专利技术提供了一种基于Krylov子空间的深部地层导热系数三维预测方法及装置，包括以下步骤：在均匀半空间研究区域内构建导热系数异常体，设定该研究区域边界条件并开展有限元温度数值模拟，获取地下空间三维温度场d

Three dimensional prediction method and device of thermal conductivity of deep formation based on Krylov subspace

【技术实现步骤摘要】
基于Krylov子空间的深部地层导热系数三维预测方法及装置


[0001]本专利技术涉及导热系数预测领域，尤其涉及一种基于Krylov子空间的深部地层导热系数三维预测方法及装置。

技术介绍

[0002]地热资源作为一种稳定可持续的可再生清洁能源，目前正受到全社会前所未有的高度重视，开发利用地热资源对于助力“双碳”目标具有重要的现实意义及深远影响。深部岩石的热物性如导热系数、比热容以及放射性生热率等直接影响到地球内部各个圈层间岩石的热生成、热储存和热传递，是地表、地球内部温度分布及热传递研究不可缺少的参数。岩石热物性也是定量评价温度随时间变化的先决条件，因为它们直接影响地下传热过程或温度变化。此外，在地源热泵设计以及地热开发工程中，也需要利用热物性参数来定量评价地源热泵系统及换热工质等与围岩地质体之间的热交换。
[0003]导热系数作为岩石矿物固有的物理属性，是各种热物性参数中最重要的参数之一。不同地层中，暴露出不同岩石矿物学类型，岩石导热系数在横向和纵向上都会发生明显变化，从而改变局部和区域的热构造演化、热结构及地热地质特征。目前获取地下岩石导热系数的方法大致可以分为两大类：(1)通过实验室对钻屑或岩心样品测量来确定；(2)通过地球物理测井来刻画。但这些方法仅能刻画点尺度或者线尺度的岩石导热系数，无法预测区域性导热系数展布特征。然而目前对于深部地层导热系数的三维预测尚处于空白阶段，直接制约研究区地热资源热成因机制解译及高温地热靶区选址。

技术实现思路

[0004]针对深部地层导热系数三维预测领域的空缺状态，本专利技术提出了一种基于Krylov子空间的深部地层导热系数三维预测方法及装置。以地下空间温度场作为观测数据，地下热流值及各地层热物性参数等作为先验信息，开展有限元温度数值模拟并构造正则化目标函数，基于Gauss
‑
Newton算法构造Hessian矩阵和梯度向量，利用Jacobian
‑
free Krylov子空间技术求解Jacobian矩阵与任意向量的乘积避免大型稠密的Jacobian矩阵的求解及存储，L
‑
BFGS算法近似求解Hessian矩阵的逆矩阵，用于减少计算量并获取模型修正量，Wolfe准则搜索模型步长更新模型参数，循环预测使得实测数据与模拟数据拟合差小于预设值，输出模型即为最优预测结果，实现深部介质导热系数的三维预测。
[0005]根据本专利技术的第一方面，本专利技术提供了一种基于Krylov子空间的深部地层导热系数三维预测方法，包括以下步骤：
[0006]在均匀半空间研究区域内构建导热系数异常体，对所述均匀半空间研究区域设定边界条件，结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场T(x,y,z)＝d
obs
，x、y、z分别表示x、y、z轴方向；
[0007]根据所述地下空间温度场d
obs
，结合在均匀半空间研究区域先验信息下获取的地下空间温度场T(x,y,z)＝d(m)，构造初始预测模型和正则化目标函数，采用Gauss
‑
Newton
算法并忽略所述正则化目标函数的二阶信息项，构造Hessian矩阵及梯度向量；由Jacobian
‑
free Krylov子空间技术求解Hessian矩阵中Jacobian矩阵及所述Jacobian矩阵的转置与任意向量的乘积；获取所述Hessian矩阵的逆矩阵并求解所述Hessian矩阵及梯度向量组成的线性方程组得到预测模型修正量Δm；
[0008]根据所述预测模型修正量Δm，利用Wolfe准则获得预测模型搜索步长以更新预测模型参数；
[0009]将更新的预测模型参数代入有限元温度数值模拟得到当前预测模型参数下的模拟数据，并计算与实测数据之间的数据拟合差，若小于预设拟合差条件，则输出当前预测模型参数作为最优预测结果；否则，返回更新d(m)并求解新的预测模型修正量Δm。
[0010]优选地，所述边界条件包括：
[0011]上边界温度为恒定温度T
up
，满足：
[0012]T＝T
up
[0013]四周边界为绝热边界，满足绝热边界条件：
[0014][0015]下边界为热流值q
down
，满足热流值条件：
[0016][0017]其中，T为温度，单位℃，n为边界法线方向，无量纲，k为导热系数，单位W/m
·
K，q
down
为热流值，单位W/m2。
[0018]优选地，所述有限元温度数值模拟，采用六面体结构离散方式，六面体单元内温度T、导热系数k均采用线性插值：
[0019][0020][0021]其中，i为节点标号，T
i
和k
i
分别为单元内各节点的温度和导热系数，单位分别为℃和W/m
·
K；N
i
为形函数，满足：
[0022][0023]其中，ξ
i
、η
i
、ζ
i
为节点i在母单元上的坐标，无量纲；母单元中ξ、η、ζ与子单元中坐标x、y、z关系公式：
[0024][0025]其中，x0、y0、z0分别为子单元在x、y、z轴方向上的中点；a、b、c分别为子单元在x、y、z方向上的边长，单位为m。
[0026]优选地，所述结合导热系数分布和所述边界条件开展有限元温度数值模拟，获取地下空间温度场的步骤，包括：
[0027]给定地下空间温度场的微分方程为：
[0028][0029]对所述地下空间温度场的微分方程乘以冲激函数δT并积分：
[0030][0031]根据场论中哈密顿算子运算规则及奥高公式，并代入边界条件，得到地下空间温度场的积分方程为：
[0032][0033]其中，Ω为研究区域，Γ
up
为上边界，Γ
down
为下边界，q
down
为热流值，k表示导热系数，单位W/m
·
K，T为温度，单位℃，为散度，为梯度，δT为冲激函数；
[0034]对所述地下空间温度场的积分方程进行离散剖分、线性插值，求解得到积分方程内各项的单元积分：
[0035][0036][0037][0038]其中，δT
e
为在单元内的冲激函数数组，T
e
为在单元内的温度数组，K
1e
、K
2e
为在单元内的刚度矩阵，P
1e
、P
2e
为在单元内的源向量，求解得到刚度矩阵和源向量内系数：
[0039][0040][0041][0042][0043][0044]其中，k
1ij
、k
2ij
、p
1ij
、p
2ij
分别为刚度矩阵K
1e...

【技术保护点】

【技术特征摘要】
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/m2。3.如权利要求1所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述有限元温度数值模拟，采用六面体结构离散方式，六面体单元内温度T、导热系数k均采用线性插值：k均采用线性插值：其中，i为节点标号，T
i
和k
i
分别为单元内各节点的温度和导热系数，单位分别为℃和W/m
·
K；N
i
为形函数，满足：
其中，ξ
i
、η
i
、ζ
i
为节点i在母单元上的坐标，无量纲；母单元中ξ、η、ζ与子单元中坐标x、y、z关系公式：其中，x0、y0、z0分别为子单元在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/m2；根据每个单元各个节点的编号并扩展至与研究区域相对应的元素位置处，即得合成后的总体刚度矩阵和源向量的线性方程组：∑
Ω
[δT
T
(K1+K2)T
‑
δT
T
(P1+
P2
)]＝0其中，δT为研究区域冲激函数数组，T为研究区域内的温度数组，K1、K2为研究区域内的刚度矩阵，P1、P2为研究区域内的源向量，略去δT
T
，化简后有：KT＝P其中，K＝K1+K2，P＝P1+P2；经过求解该线性方程组KT＝P，实现有限元数值模拟后，得到地下空间温度场T(x,y,z)＝d
obs
，x、y、z分别表示x、y、z轴方向。5.如权利要求4所述的基于Krylov子空间的深部地层导热系数三维预测方法，其特征在于，所述正则化目标函数Φ的表达式为：Φ＝||D[d(m)
‑
d
obs
]||2+λ||W(m
‑
m
apr
)||2其中，D为数据加权矩阵；d
obs
为观测数据，即地下空间温度场T(x,y,z)；m为当前预测模型参数向量；d(m)为在当前预测模型参数下的模拟数据；W为光滑约束矩阵；λ为阻尼因子；m<...

【专利技术属性】
技术研发人员：杨健，胡祥云，黄国疏，
申请(专利权)人：中国地质大学武汉，
类型：发明
国别省市：