冻土水热耦合介质步进频率探地雷达响应特征

张世榕; 孙怀凤; 葛峻恺; 刘锐; 董英凯

doi:10.12404/j.issn.1671-1815.2600092

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (25) :10717-10728 DOI: 10.12404/j.issn.1671-1815.2600092

天文学、地球科学

冻土水热耦合介质步进频率探地雷达响应特征

张世榕 ¹^,²^,³ ,
孙怀凤 ¹^,²^,³ ,
葛峻恺 ¹^,²^,³ ,
刘锐 ¹^,²^,³ ,
董英凯 ¹^,²^,³^,^*

作者信息 +

Response Characteristics of Stepped Frequency Ground Penetrating Radar in Hydro-thermal Coupled Frozen Soil Media

ZHANG Shi-rong ¹^,²^,³ ,
SUN Huai-feng ¹^,²^,³ ,
GE Jun-kai ¹^,²^,³ ,
LIU Rui ¹^,²^,³ ,
DONG Ying-kai ¹^,²^,³^,^*

Author information +

文章历史 +

摘要

冻土内部的冰水分布与电性参数密切相关,是探地雷达进行冻土探测的关键因素。然而传统探地雷达正演建模方法缺乏从热力学角度刻画冻土冻融分层演化与电磁响应之间关系的方法。基于此,构建了基于冻土水热耦合理论与有效介质理论的综合探地雷达正演模型框架,揭示了冻土在不同热力状态下的电性结构及其步进频率探地雷达响应特征。首先利用 COMSOL 建立温度场与水分场耦合模型,获得冻土活动层在不同冻融阶段的温度、饱和度、含冰率分布;随后基于有效介质理论计算冻土模型的等效介电常数与电导率;最后将电性参数导入 gprMax 软件,并采用步进频率连续波(stepped frequency continuous wave,SFCW)与Ricker子波进行正演模拟对比分析。结果表明,基于水热耦合理论的冻土建模能够重构冻土内部冻结与融化锋面的动态演化过程;对比分析显示,SFCW信号与Ricker子波在界面位置和走势上保持一致且未出现明显波形畸变,表明了SFCW能够有效识别冻土中由温度和水分剧烈变化引起的复杂介电分层结构。该研究为冻土区步进频率探地雷达冻土活动层探测提供了正演参考。

Abstract

The distribution of ice and water within frozen soil is closely correlated with electrical parameters and is a key factor for frozen soil detection by ground penetrating radar (GPR). However, the relationship between the evolution of freeze-thaw layering and electromagnetic response fails to be characterized by traditional GPR forward modeling methods from a thermodynamic perspective. To address this issue, a comprehensive GPR forward modeling framework based on frozen soil hydro-thermal coupling theory and effective medium theory was proposed. The electrical structure of frozen soil under different thermal states and its stepped frequency continuous wave (SFCW) GPR response characteristics were revealed. First, a coupled model of temperature and moisture fields was established using COMSOL to obtain the distribution of temperature, saturation, and ice content of the active layer of frozen soil at different freeze-thaw stages. Subsequently, the equivalent dielectric constant and electrical conductivity of the frozen soil model were calculated based on effective medium theory. Finally, the electrical parameters were imported into gprMax software, and forward modeling comparative analysis was conducted using SFCW signals and Ricker wavelets. It is indicated by the results that the dynamic evolution process of freezing and thawing fronts within frozen soil can be reconstructed by the proposed modeling method. It is demonstrated that SFCW signals and Ricker wavelets remain consistent in interface position and trend without obvious waveform distortion. It is confirmed that complex dielectric layered structures in frozen soil caused by drastic changes in temperature and moisture can be effectively identified by SFCW. A forward modeling reference is provided for the detection of active layers in permafrost regions using stepped frequency GPR.

Graphical abstract

关键词

探地雷达 / 冻土水热耦合 / 有效介质理论 / 步进频率连续波 / 正演分析

Key words

ground penetrating radar / frozen soil hydro-thermal coupling / effective medium theory / stepped frequency continuous wave / forward modeling analysis

引用本文

引用格式 ▾

张世榕, 孙怀凤, 葛峻恺, 刘锐, 董英凯. 冻土水热耦合介质步进频率探地雷达响应特征[J]. 科学技术与工程, 2026, 26(25): 10717-10728 DOI:10.12404/j.issn.1671-1815.2600092

ZHANG Shi-rong, SUN Huai-feng, GE Jun-kai, LIU Rui, DONG Ying-kai. Response Characteristics of Stepped Frequency Ground Penetrating Radar in Hydro-thermal Coupled Frozen Soil Media[J]. Science Technology and Engineering, 2026, 26(25): 10717-10728 DOI:10.12404/j.issn.1671-1815.2600092

登录浏览全文

4963

注册一个新账户忘记密码

冻土是指0 ℃以下,含有冰的各种岩石和土壤^[1-2]。一般可分为季节冻土以及多年冻土^[3]。多年冻土广泛分布于高纬度和高海拔地区,其物理性质对气候变化和工程建设具有重要影响^[4]。多年冻土区的活动层是指位于多年冻土层之上、地表下一定深度内暖季融化、冷季完全回冻的近地表层^[5]。冻土为四相体系,相较于普通土,冻土中冰的存在影响着冻土的工程性质,气候的季节变化引起冻土中温度的变换,而这一过程导致了水的液相与固相的转变,引起冻土性质的变化^[6],从而导致多年冻土具有热稳定性差、地下冰含量丰富、对气候变暖反应极为敏感以及水热活动强烈等特点。地温升高导致多年冻土活动层退化的现象日益显著,由此带来冻融、地表沉降、热融沉陷、滑塌等地质灾害^[7]。因此,开展冻土精细化探测,揭示其内部结构及动态演化规律,对寒区工程建设、灾害防控及环境保护具有重要的科学价值与实践意义^[8-9]。

目前,探地雷达已被广泛应用于多年冻土上限深度识别、活动层厚度变化监测以及地下冰层结构探测等领域^[10-12]。然而,脉冲探地雷达受限于中心频率带宽,往往难以在探测深度与分辨率之间取得平衡。相比之下,步进频率连续波(stepped frequency continuous wave,SFCW)探地雷达通过发射步进频率信号并在频率域接收反射回波,具有更高的信号带宽,因此具有广泛的发展前景^[13]。尽管如此,受现场条件、天线布设、环境干扰及反演精度的限制,实测数据解释往往存在不确定性^[14]。基于探地雷达的正演仿真研究,可以通过构建理想化或半真实模型,分析SFCW在被测介质中传播规律,从而为数据解释与反演提供参考,减少实测试验成本,提高探测效率^[15-16]。

冻土是一种受环境温度、水热条件及外部荷载长期影响而不断演化的介质,其内部含冰率与含水率随时间动态变化^[17]。目前关于探地雷达正演模型的研究大多集中在介质的随机分布,少有针对介质随时间变化特征的研究^[18-19]。冻土的动态变化是复杂的物理过程。目前,针对冻土的研究主要是水热耦合理论^[20-22]。冻土动态变化过程中存在热传导、水分迁移及冰水相变过程^[23-24]。水热耦合理论能描述冻土在不同温度场和水分场条件下的相态分布及演变规律,为构建真实的冻土动态物理模型提供了理论基础^[25-26]。

冻土内部的冰水分布与电性参数密切相关,而常规冻土介质探地雷达正演建模方法缺乏由热力学机理出发的多场耦合分析,难以真实刻画冻土在气候交替作用下内部水热状态的动态演化过程及其对电磁波传播的时变影响。针对这一局限,现构建基于冻土水热耦合理论与有效介质理论的综合探地雷达正演模型框架。该方法突破传统静态建模的单一性,从热力学机理出发重构冻土冻融锋面的动态演化过程。通过对比分析SFCW与Ricker子波的回波信号,揭示宽带雷达信号在复杂冻土介质中的响应特征。以期为寒区工程步进频率探地雷达探测提供理论参考与技术支持,对寒区工程的长期监测与灾害预警具有一定的实用价值。

1 基本理论

1.1 冻土水热耦合理论

在冻土冻胀研究中,冻土多场耦合模型应用广泛。该理论认为,冻土变化是由温度场、水分场、应力场共同作用的结果。其中温度场通过控制水分子的相变与迁移,直接影响含水率、含冰率、应力-应变的变化,并起决定性作用;水分场控制水分子的渗透、迁移,直接影响含水率与应力应变;应力场通过应力-应变的变化改变土体结构,从而改变土体导热率和渗透性,间接影响温度场和水分场。为简化计算,对冻土模型进行假设,认为土颗粒和冰分子不可压缩且土体形变为弹性形变。这样温度场通过温度变化来改变含水率;水分场通过含水率变化来改变传热系数;在温度场和水分场共同控制下通过含冰率变化来改变土体应变。简化后的冻土多场耦合关系如图1所示。

在确定水分场和温度场后,由于假设土体变形是弹性的,因此可忽略应力场改变对水分场、温度场的影响,不研究冻土的冻胀与应力-应变关系,故仅考虑水热耦合作用。

1.1.1 温度场

在前述假设的基础上推导温度场控制方程。根据能量守恒定理,冻土内部任意一处的内能变化等于外部传递热量与土体内部热量变化(冰水相变)之和。则冻土温度场控制方程如式(1)所示^[27]。

(1)$\mathit{\rho }\mathit{C}\left(\mathit{\theta }\right)\frac{\partial \mathit{T}}{\partial \mathit{t}}=\mathit{\lambda }\left(\mathit{\theta }\right){\nabla }^{2}\mathit{T}+\mathit{L}{\mathit{\rho }}_{\mathrm{I}}\mathrm{ }\frac{\partial {\mathit{\theta }}_{\mathrm{I}}}{\partial \mathit{t}}$

式(1)中:ρ为冻土的密度,kg/m³;C为冻土等效热容,J/(kg·K);θ为体积含水率;T为温度,℃;t为时间;λ为导热系数,W/(m·℃);L为水分子相变潜热,kJ/kg,取334.56;θ_I为体积含冰率;ρ_I为冰的密度,kg/m³。

温度场控制方程中等号左侧为内能变化量,等号右侧为热传递量与相变热量差。冻土等效热容由土体的热容、水的热容、冰的热容共同控制,采用加权平均的方式计算,冻土体积含水率计算如式(2)所示,等效热容计算如式(3)所示。

(2)$\mathit{\theta }={\mathit{\theta }}_{\mathrm{w}}+\frac{{\mathit{\rho }}_{\mathrm{I}}}{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\theta }}_{\mathrm{I}}$

(3)$\mathit{C}\left(\mathit{\theta }\right)=\frac{{\mathit{C}}_{\mathrm{S}}+{\mathit{\theta }}_{\mathrm{w}}{\mathit{C}}_{\mathrm{w}}+{\mathit{\theta }}_{\mathrm{I}}{\mathit{C}}_{\mathrm{I}}}{1+\mathit{\theta }}$

式中:θ_w为未冻水的体积含水率;C_S为土体的比热;C_w为水的比热,取4 200 J/(kg·K);C_I为冰的比热,取2 100 J/(kg·K)。

土体导热系数表示土体传递热量的效率,土体导热系数的计算方法可以分为直接加权平均和指数加权平均两种形式。直接加权平均法相较于指数加权平均法误差较大^[28],如式(4)所示。

(4)$\mathit{\lambda }=\mathit{\lambda }{\mathit{C}}_{\mathrm{S}}^{{\mathit{\theta }}_{\mathrm{S}}}{\mathit{\lambda }}_{\mathrm{w}}^{{\mathit{\theta }}_{\mathrm{w}}}{\mathit{\lambda }}_{\mathrm{I}}^{{\mathit{\theta }}_{\mathrm{I}}}{\mathit{\lambda }}_{\mathrm{a}}^{{\mathit{\theta }}_{\mathrm{a}}}$

式(4)中:θ_S为土颗粒的体积分数;θ_a为空气的体积分数;λ_S为土体的导热系数;λ_w为水的导热系数;λ_I为冰的导热系数;λ_a为空气的导热系数。

1.1.2 温度场

在前述假设的基础上推导水分场控制方程。冻土中的未冻水的迁移变化遵循达西定律,根据质量守恒定律,冻土内部任意一处的水分的迁移量等于水分渗透量与冰水相变的损耗之和,考虑重力作用,水分场控制方程及相关参数如式(5)~式(7)所示。

(5)$\frac{\partial {\mathit{\theta }}_{\mathrm{w}}}{\partial \mathit{t}}+\frac{{\mathit{\rho }}_{\mathrm{I}}}{{\mathit{\rho }}_{\mathrm{w}}}\frac{\partial {\mathit{\theta }}_{\mathrm{I}}}{\partial \mathit{t}}=\nabla \left[\mathit{D}\right({\mathit{\theta }}_{\mathrm{w}})\nabla {\mathit{\theta }}_{\mathrm{w}}]+\frac{\partial {\mathit{k}}_{\mathrm{g}}\left({\mathit{\theta }}_{\mathrm{w}}\right)}{\partial \mathit{z}}$

(6)$\mathit{D}\left({\mathit{\theta }}_{\mathrm{w}}\right)=\frac{\mathit{k}\left({\mathit{\theta }}_{\mathrm{w}}\right)}{\mathit{c}\left({\mathit{\theta }}_{\mathrm{w}}\right)}\mathit{I}$

(7)$\mathit{I}={10}^{-10{\mathit{\theta }}_{\mathrm{I}}}$

式中:D(θ_w)为冻土中水的扩散率;k_g(θ_w)为重力方向水的渗透率^[29],m/s;k(θ_w)为冻土的渗透率,m/s;c(θ_w)为比水容量,m^-1;I为阻抗因子,为冻土中冰对未冻水渗透的作用。

采用VG模型确定冻土的渗透率与比水容量,各项参数表达式如式(8)~式(10)所示。

(8)$\mathit{S}=\frac{{\mathit{\theta }}_{\mathrm{w}}-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$

(9)$\mathit{c}\left({\mathit{\theta }}_{\mathrm{w}}\right)=\frac{{\mathit{a}}_{0}\mathit{m}}{1-\mathit{m}}{\mathit{S}}^{\frac{1}{\mathit{m}}}(1-{\mathit{S}}^{\frac{1}{\mathit{m}}}{)}^{\mathit{m}}$

(10)$\mathit{k}\left({\mathit{\theta }}_{\mathrm{w}}\right)={\mathit{k}}_{\mathrm{s}}\sqrt{{\mathit{S}}_{\mathrm{e}}}[1-(1-{\mathit{S}}^{\frac{1}{\mathit{m}}}{)}^{\mathit{m}}{]}^{2}$

式中:S为冻土的饱和度;θ_s为冻土饱和含水率;θ_r为冻土残余含水率;m、a₀为由土体性质确定的参数;k_s为饱和土渗透系数^[30-31]。

在温度场控制方程和水分场控制方程之间增加冰与水的体积之比,固液比B(T)作为耦合项。固液比计算如式(11)所示。

(11)$\mathit{B}\left(\mathit{T}\right)=\frac{{\mathit{\theta }}_{\mathrm{I}}}{{\mathit{\theta }}_{\mathrm{w}}}=\left\{\begin{array}{ll}1.1{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{f}}}\right)}^{\mathit{B}}-1.1,& \mathit{T}<{\mathit{T}}_{\mathrm{f}}\\ 0,& \mathit{T}\ge {\mathit{T}}_{\mathrm{f}}\end{array}\right.$

式(11)中:T_f为土体冻结温度;θ_w为未冻水的体积含水率;θ_I为体积含冰率;B为土盐常数,按土质取砂土0.61,粉土0.47,黏土0.5。

1.2 有效介质理论

有效介质理论是一种用于简化复杂多相介质问题的理论方法,其核心思想是将非均质或非连续的多相介质等效为一个均匀的连续介质,从而便于分析和计算。这种方法适用于介质的宏观特性远大于微观结构尺度的情况。在宏观尺度上,冻土为各向同性的混合物,因此适用于有效介质理论。本文采用Bruggeman方程,介电常数计算方程如式(12)所示,电导率计算方程如式(13)所示。

(12)$\stackrel{\mathit{N}}{\sum _{\mathit{i}=1}}{\mathit{\theta }}_{\mathit{i}}\mathrm{ }\frac{{\mathit{\varepsilon }}_{\mathit{i}}-{\mathit{\varepsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}}}{{\mathit{\varepsilon }}_{\mathit{i}}+2{\mathit{\varepsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}}}=0$

(13)$\stackrel{\mathit{N}}{\sum _{\mathit{i}=1}}{\mathit{\theta }}_{\mathit{i}}\frac{{\mathit{\sigma }}_{\mathit{i}}-{\mathit{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}}{{\mathit{\sigma }}_{\mathit{i}}+2{\mathit{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}}=0$

式中:θ_i为复合材料中第i相的体积分数;ε_i为第i相的相对介电常数;ε_eff为N相复合材料的等效相对介电常数;σ_i为第i相的电导率;σ_eff为N相复合材料的等效电导率。

1.3 步进频率连续波正演原理

探地雷达利用高频电磁波进行探测,这一过程主要用时域有限差分(finite-difference time-domain,FDTD)来求解麦克斯韦方程组进行正演模拟。FDTD通过离散麦克斯韦方程组中的旋度方程,采用时间和空间在电场和磁场交替间隔取样的方式,并沿时间逐步推进,进行方程求解,从而获取空间域电磁场的分布特征。目前,常用探地雷达正演软件 gprMax 完成^[15-16]。

对于发射SFCW的雷达系统, gprMax 无法对一段完整的 SFCW 信号直接正演,而正演过程可视为一个线性时不变系统,因此可以采用冲激响应法将 SFCW 信号与模型系统的冲激响应h(t)卷积,即可得到对应 SFCW 信号的回波信号^[32-33]。对于一段有N个频率点的 SFCW 信号,第n个频点的发射频率为

(14)f_n=f₀+(n-1)Δf

式(14)中:f_n为第n个频点的频率;f₀为起始频率;Δf为频率步长。第n个频点的发射信号为

(15)T_x(t)=sin(2πf_nt)

将式(15)与冲激响应进行卷积,并考虑电磁波在介质中的传播损耗与延时,得到第n个频点的目标体回波信号为

(16)$\begin{array}{c}{\mathit{R}}_{\mathit{x}}\left(\mathit{t}\right)={\mathit{T}}_{\mathit{x}}\left(\mathit{t}\right)\mathit{h}\left(\mathit{t}\right)\approx {\mathit{A}}_{\mathit{n}}\mathit{s}\mathit{i}\mathit{n}\left[2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\left(\mathit{t}-\frac{2\mathit{R}}{\mathit{v}}\right)\right]\end{array}$

式(16)中:A_n为回波幅度;R为目标距离;v为电磁波在介质中的波速。

为了提取回波信号的幅度与相位信息,需要对接收信号进行正交解调。将回波信号分别与同相参考信号和正交参考信号进行混频处理。

(17)Ref_I(t)=sin(2πf_nt)

(18)Ref_Q(t)=cos(2πf_nt)

(19)${\mathit{I}}_{\mathrm{m}\mathrm{i}\mathrm{x}}\left(\mathit{t}\right)=\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{ }\left(2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)-\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{ }\left(4\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\mathit{t}-2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)$

(20)${\mathit{Q}}_{\mathrm{m}\mathrm{i}\mathrm{x}}\left(\mathit{t}\right)=\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{ }\left(4\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\mathit{t}-2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)-\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{ }\left(2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)$

式中:Ref_I(t)为同向参考信号;Ref_Q(t)为正交参考信号;I_mix(t)为混频后的同相分量;Q_mix(t)为混频后的正交分量。

由式(1)~式(20)可知,混频后的信号中同时包含高频倍频分量和携带目标信息的直流差频分量,通过设置归一化截止频率的低通滤波器滤除高频项后,得到最终的直流分量如下。

(21)${\mathrm{I}}_{\mathit{n}}=\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{ }\left(2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)$

(22)${\mathrm{Q}}_{\mathit{n}}=-\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{ }\left(2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)$

式中:I_n和Q_n分别为回波信号在复平面上的实部和虚部;第n个频点的频域复信号为

(23)$\mathit{K}\left(\mathit{n}\right)={\mathrm{I}}_{\mathit{n}}+\mathrm{j}{\mathrm{Q}}_{\mathit{n}}=\frac{{\mathit{A}}_{\mathit{n}}}{2}\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{ }\left(-\mathrm{j}2\mathrm{\pi }{\mathit{f}}_{\mathit{n}}\frac{2\mathit{R}}{\mathit{v}}\right)$

该复信号K(n)包含了第n个频率点的振幅衰减与相位延迟信息。对所有N个频点的K(n)进行逆傅里叶变换,即可得到SFCW的时间域回波信号。在IFFT的过程中,由于SFCW本身的频域截断导致信号产生旁瓣,因此在处理过程中使用汉宁窗作为加权函数对旁瓣进行抑制^[34]。

2 数值模拟

2.1 COMSOL 冻土建模

为进行冻土冻融过程模拟,采用COMSOL Multiphysics 软件中的偏微分方程(partial differential equation,PDE)模块建立温度场和水分场。

2.1.1 模型建立

为研究冻土在冻融过程中的温度与含水变化规律,基于 COMSOL Multiphysics 建立二维矩形区域的水热耦合模型。如图2所示,模型长10 m,深5 m,上边界代表地表,下边界代表地下基岩层,上边界3~7 m处为道路表面。几何模型网格如图3所示,模型网格大小为0.05 m×0.05 m。设置的材料温度场参数如表1所示,水分场参数如表2所示。

2.1.2 初始条件和边界条件

初始温度场设定为垂向分段线性分布。地表温度5 ℃,1 m深温度为0 ℃,底部温度-2 ℃;上边界采用狄利克雷边界条件模拟温差及季节变温。根据青藏高原气象资料年平均气温为-1.15 ℃,并考虑未来气候变暖导致温度每50年升高2.6 ℃,以正弦函数的形式模拟季节温差16 ℃,给出上边界的边界条件。地表气温变化函数如式(24)所示,车道温度变化函数如式(25)所示^[35-36]。

(24)${\mathit{T}}_{\mathrm{d}}=-1.15+16\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{2\mathrm{\pi }}{365}\mathit{d}\right)+\frac{2.6\mathit{d}}{365\times 50}$

(25)${\mathit{T}}_{\mathrm{c}}=3.15+16\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{2\mathrm{\pi }}{365}\mathit{d}\right)+\frac{2.6\mathit{d}}{365\times 50}$

式中:T_d为气温变化函数,℃;T_c为道路温度变化函数,℃;d为时间,d。

考虑地热影响,下边界采用诺伊曼边界条件,初始温度梯度为0.03 ℃/m,热流密度为-0.06 W/m²。水分场的边界条件定义为无通量,仅设置水分场的初始饱和度S₀=0.6,初始含冰率为0。

2.2 电性模型计算

冻土的电磁特性主要由介电常数与电导率决定。在得到冻土的物性特征分布后,可根据物性分布计算整体的电性参数。通过 COMSOL Multiphysics 中建立二维矩形区域的水热耦合模型,计算各个点各物质的体积分数。其中温度T、饱和度S通过PDE方程直接求出。各项体积分数如式(26)~式(29)计算。

(26)θ_w=S(θ_s-θ_r)+θ_r

(27)θ_I=B(T)θ_w

(28)θ_S=1-θ_s

(29)θ_a=θ_s-(θ_w+θ_I)

由于空气、冰的电导率数量级远小于土中的矿物组分和水,为简化计算不计空气和冰的电导率。不同土的电性参数差异主要由它们的持水能力不同导致。为简化计算,土中的矿物组分取同一个电性参数的电性参数。各物质的电性参数如表3所示。

2.3 探地雷达正演建模

通过gprMax软件进行FDTD正演,在电性模型的基础上编写可被gprMax识别的HDF5模型和对应的材料文件。HDF5模型中的值为材料代码,对应的材料参数储存在材料文件中。为了保证数值精度与文件体积的合理性,在编写前对电性模型的数据进行预处理,为了兼顾计算精度与模型效率,对连续变化的电性参数进行了离散化处理。电性模型为数据矩阵,包括介电常数矩阵和电导率矩阵。该文件以逐点扫描的方式,对每一位置对应的介电常数与电导率组合进行判定。若当前位置的电性参数在材料文件中尚未存在,则在材料文件中新增并在HDF5模型中将该材料的标识符作为值赋到对应位置上。经过上述步骤即可完成从电性参数到可被gprMax识别的模型的转化。

3 结果与分析

3.1 冻土分布

为揭示冻土路基在年周期内的动态演化规律及其非均匀冻融特征,选取一年内具有代表性的典型时刻进行温度场与水分场结果分析,结果如图4~图7所示。由图4地表气温变化函数可见,d=0 d时设定为夏季最高温,随着气温降低,土体经历降温、冻结、深冻、回暖及融化过程。其中,d=182 d时土体达到年度最低温,路基处于深冻状态;d=365 d时土体温度重新达到最高,浅层冻土完全融化。为体现冻结起始和融化起始阶段模型的剧烈相变过程,选取16个关键时间节点进行全周期分析。

冻结过程从地表开始。在降温初期,路基土体热量持续散失,地表温度逐渐跌破冻结温度,表层水分发生相变,形成初始含冰区,如图7中d=95 d所示。在冻结期(d=120~210 d),随着地表温度进一步降低,地表及浅层土体温度低于冻结温度,冻土活动层中含冰量明显增加。同时由于深部永冻土层产生的聚水效应,导致水分向永冻土层上表面聚集。此时活动层内部含水率较低,导致在温度最低时(d=210 d)的含冰率小于水分聚集前(d=150 d)的含冰率。

在融化期(d=255~365 d),气温回升,路基中心的表层土首先开始融化。但此时浅层土体仍处于冻结温度之下,如图5中d=270 d所示,融化产生的水由于低温聚水效应向下移动,当移动到温度小于冻结温度的活动层时产生原位冻结,最终导致表层融区和活动层之间产生饱冰带。随着温度的不断升高,融化产生的水越来越多并不断向下富集,如图7中d=300~365 d所示,饱冰带下移且增厚。这个过程中融区中的土体处于饱和状态,融区厚度不断增大,如图6中d=300~365 d所示。

3.2 电性参数变化

根据有效介质理论,由式(26)~式(29)计算得到的各相体积分数,结合表3中的电性参数计算模型的等效介电常数和电导率分布。模型等效介电常数变化如图8所示,模型等效电导率变化如图9所示。

在冻土介质中,水相决定了冻土整体的电性参数分布。冻结过程中,水相逐渐转化为冰相,冰的介电常数为3.15,而水为81,导致等效介电常数显著下降;融化过程中冰相转化为水相,使介电常数升高。不同土质中介电常数变化幅度不同,这与其保水能力相关。对于电导率,由于各个组分电导率存在数量级的差异,其中水的电导率最高,因此电导率分布直接反映了含水率的分布区间。

3.3 探地雷达响应特征

为了分析冻土在不同阶段的电性差异对探地雷达信号的影响,采用gprMax对不同时刻的等效电性参数模型进行了二维电磁波场B-scan正演模拟。模型网格大小为0.01 m×0.01 m×0.01 m,天线间隔1 m,时窗设置为 200 ns,测线起点为x=1 m,测线终点为x=9 m。设置 SFCW 起始频率为100 MHz,步进点数为1 001,频率步进值为1 MHz,模型时域响应波形如图10所示。为进行对比分析,采用中心频率为 200 MHz 的 Ricker 子波进行正演,其他参数相同,结果如图11所示。

通过图10和图11对比分析可见,SFCW信号与Ricker子波信号在界面位置和走势上保持高度一致,且均能清晰地反映出冻土冻融过程中的主要界面特征。在d=95 d时,随着冻结锋面的下移,两幅B-scan图像均准确勾勒出了未冻土与冻土之间的介电常数差异界面;从d=270 d到d=365 d的融化阶段,冻融锋面(饱冰带)下降过程也有明显的反射,并未出现明显的波形畸变,验证了 SFCW 方法在冻土探测中的有效性。

相比于Ricker子波,经过加权函数修正后的SFCW信号在时域上表现出良好的脉冲压缩特性,其反射同相轴更加细腻、连续。在d=95 d和d=120 d的冻结阶段,SFCW图像对浅层冻结界面的刻画展示出更高的垂直分辨率,这说明 SFCW 信号能够有效识别冻土中由温度和水分剧烈变化引起的复杂介电常数分层,具有优异的层界识别能力。

注意到在d=365 d时,SFCW结果和Ricker子波图像均出现双曲线,结合图8中d=365 d分析,这是由于冻融锋面处存在的局部介电常数非均匀突变点电磁波衍射现象,进一步表明了SFCW和Ricker子波对微小非均匀体均具有敏锐的捕捉能力,能够为寒区工程监测中的实际数据解释提供可靠的理论依据。

3.4 工程实例

为进一步验证本文模型在实际寒区工程中的实用性与正确性,选取青藏高原 G109 国道五道梁段K3017+528 处进行探地雷达冻土探测,测线布置方向垂直于 G109 走向,横跨路基,车道位于测线中央。探地雷达采用美国 GSSI SIR-4000 型号主机,采用的天线中心频率为100 MHz,以时间采集的方式进行拖曳式测量。设置探测时间窗为 200 ns,由前期 CMP 数据计算得地下波速为 0.08 m/s。

图12展示了探地雷达探测结果与地质解释的对比。其中,图12(a)为现场工作照片,图12(b)为探地雷达探测结果,图12(c)为层位追踪结果与钻孔数据对比。

由图12(b)可见,在深度2~3 m 处存在一条清晰、连续的强反射同相轴,结合图12(c)中的层位追踪结果与现场钻孔勘探资料对比分析可知。该强反射同相轴对应于冻土的冻融界面。在路基两侧的天然场地,冻融界面埋深为 2~3 m,这与钻孔揭示的活动层的分层结构吻合良好。在车道正下方区域,同相轴出现明显的中断与不连续现象。这是由于路面吸热及车辆荷载产生的热扰动导致路基下伏多年冻土上限发生改变,局部融化深度增加,破坏了原有的层状结构。

从物理机制分析,探测处于暖季,浅层活动层发生融化,而深部多年冻土层仍保持冻结状态。在重力势与温度梯度的双重驱动下,土壤水分向冻结锋面迁移并聚集,导致界面处未冻水含量显著升高。由于液态水与冰、土骨架之间存在巨大的介电常数差异,该富水层在雷达图像上形成了显著的强反射界面。

实测探地雷达图像中冻融界面表现为清晰、连续的强反射并可追踪到同相轴,这一特征与3.3节正演结果相符。由图12(c)可见,在车道两侧,地面以下 2~3 m 范围内存在强反射,推断地面以下 2~3 m 为融区,融区下限在地下 3 m 处;在车道位置,由于车道温度荷载等因素,同相轴于此处中断。同时,通过图12(c)的探地雷达层位追踪结果与钻孔数据高度对比可见,探地雷达结果能较好地反映地下冻融界面位置。

综上所述,实测数据的波形特征、界面深度及反射强度规律,与3.1 节构建的水热分布模型及3.3节的正演模拟结果基本一致,验证了本文模型在寒区工程探测中的可靠性与实用性。

4 结论

基于冻土水热耦合机理构建数值模型,结合有效介质理论构建冻土在不同阶段的电性模型,并采用冲激响应法通过gprMax探地雷达正演软件计算步进频率探地雷达响应特征。主要结论如下。

(1)冻土水热耦合模型结果表明,冻土中温度梯度驱动的冻结与融化过程导致土体内部矿物组分、水和空气在时间与深度方向上发生变化,从而导致冻土介电常数、电导率的变化。模拟结果展现了冻结与融化阶段的聚水效应及冻融锋面特征,且雷达响应特征与电性参数的变化规律保持一致,验证了水热耦合建模的有效性。

(2)在计算电性参数的过程中使用有效介质理论,有效介质理论中的Maxwell-Garnett(MG)方程与Bruggeman方程应用广泛。相较于MG方程,Bruggeman方程具有对称性,无需确定某一组分为基质。在冻结阶段存在冰的体积分数大于矿物组分的体积分数的情况,此时MG方程的基质假设失效。因此Bruggeman模型更适于描述冻土这类多相并存的混合物体系。

(3)提出一种新型的探地雷达电性参数建模方法,利用冻土水热耦合结果计算电性模型。相较于目前探地雷达正演模型的静态介质或随机介质建模方法,本研究可模拟冻土模型随时间变化情况,使建模更加贴合自然界真实的冻土环境,为深入分析冻土区探地雷达响应特征提供参考。

(4)对比了SFCW和Ricker子波两种不同体制的探地雷达在冻土探测的响应特征,并结合实际工程案例进行了验证。结果表明,SFCW结果与Ricker子波结果在界面位置和走势上保持高度一致,且未出现明显波形畸变,证明了SFCW在探测复杂冻土分层中的适用性。同时数值模拟结果与实测探地雷达数据在冻融界面位置及反射特征上吻合良好,这一实例验证不仅证实了本文模型的有效性,也为步进频率雷达在寒区工程监测中的数据解释提供了一定的理论与实践依据。

(5)主要考虑水热耦合,未考虑冻土实际环境中的盐分场与应力场。其中盐分场不仅影响冰点变化和水动力过程,也直接改变液相电导率,从而影响电性反演精度;应力场则与冻土冻胀有关,对于冻土区工程分析必不可少。同时本文模型采用一维垂向渗流假设,未显示包含融雪径流、表层径流及侧向水分迁移等真实水动力过程。未来的研究可从完善冻土多场耦合机制和多维度水分迁移方向入手改进,进一步提高探地雷达模型与实际冻土环境的吻合度。此外,本文研究暂未考虑介电常数的频变特性,未来研究将引入Debye模型或Cole-Cole模型来更精确描述宽带响应。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	彭小清, 田伟伟, 李璇佳, 等. 青藏高原和环北极冻土变化研究进展[J]. 冰川冻土, 2023, 45(2): 521-534.

[2]	Peng Xiaoqing, Tian Weiwei, Li Xuanjia, et al. Research progress of permafrost changes in the Qinghai-Tibet Plateau and the Circum-Arctic[J]. Journal of Glaciology and Geocryology, 2023, 45(2): 521-534.

[3]	高英, 马艳霞, 杨丰华, 等. 季节冻土区黄土湿陷变形及微观机理[J]. 科学技术与工程, 2024, 24(34): 14778-14786.

[4]	Gao Ying, Ma Yanxia, Yang Fenghua, et al. Collapsible deformation and microscopic mechanism of loess in seasonal frozen soil regions[J]. Science Technology and Engineering, 2024, 24(34): 14778-14786.

[5]	Shen J J, Wang Q, Chen Y T, et al. Evolution process of the microstructure of saline soil with different compaction degrees during freeze-thaw cycles[J]. Engineering Geology, 2022, 304: 106699.

[6]	Ma Y M, Hu Z Y, Xie Z P, et al. A long-term (2005-2016) dataset of hourly integrated land-atmosphere interaction observations on the Tibetan Plateau[J]. Earth System Science Data, 2020, 12(4): 2937-2957.

[7]	罗栋梁, 金会军, 吴青柏, 等. 天然状态下多年冻土区活动层厚度研究进展与展望[J]. 冰川冻土, 2023, 45(2): 558-574.

[8]	Luo Dongliang, Jin Huijun, Wu Qingbai, et al. Progress and prospect of active layer thickness research in permafrost regions under natural state[J]. Journal of Glaciology and Geocryology, 2023, 45(2): 558-574.

[9]	杜发, 李宗省, 贡觉扎西, 等. 多年冻土地下水及其变化研究进展[J]. 冰川冻土, 2023, 45(6): 1789-1803.

[10]	Du Fa, Li Zongsheng, Gongjue Zhaxi, et al. Research progress of groundwater and its changes in permafrost regions[J]. Journal of Glaciology and Geocryology, 2023, 45(6): 1789-1803.

[11]	陈珊珊, 臧淑英, 孙丽. 东北多年冻土退化及环境效应研究现状与展望[J]. 冰川冻土, 2018, 40(2): 298-306.

[12]	Chen Shanshan, Zang Shuying, Sun Li. Status and prospect of permafrost degradation and its environmental effects in Northeast China[J]. Journal of Glaciology and Geocryology, 2018, 40(2): 298-306.

[13]	杨沙, 王琪洁, 李佳晨, 等. 顾及土壤热传导-对流效应的InSAR冻土活动层厚度估计[J]. 地球物理学报, 2024, 67(7): 2568-2581.

[14]	Yang Sha, Wang Qijie, Li Jiachen, et al. Estimation of active layer thickness of permafrost by InSAR considering soil thermal conduction-convection effect[J]. Chinese Journal of Geophysics, 2024, 67(7): 2568-2581.

[15]	孙超, 宋韬, 郭浩天, 等. 温度对季冻区粉质黏土强度及变形特性的影响[J]. 科学技术与工程, 2024, 24(4): 1609-1618.

[16]	Sun Chao, Song Tao, Guo Haotian, et al. Effect of temperature on strength and deformation characteristics of silty clay in seasonal frozen regions[J]. Science Technology and Engineering, 2024, 24(4): 1609-1618.

[17]	Huang Y H, Niu F J, He P F, et al. Effectiveness evaluation of cooling measures for express highway construction in permafrost regions based on GPR and ERT[J]. Cold Regions Science and Technology, 2024, 228: 104339.

[18]	Tregubov O D, Uyagansky K K. Substantiation of the monitoring network of talik zones in urbanized permafrost areas based on gpr profiling data (Anadyr, Chukotka)[J]. Urban Science, 2024, 8(3): 94.

[19]	Ahmad F, Bulik A, Muto A. Dual-scan dual-frequency GPR imaging for mapping the permafrost table[C]// Radar Sensor Technology XXVIII. Washington D C: SPIE, 2024: 3017833.

[20]	熊思宇, 晋良念. 步进频率探地雷达快速超分辨成像方法[J]. 雷达科学与技术, 2025, 23(3): 253-261.

[21]	Xiong Siyu, Jin Liangnian. Fast super-resolution imaging method for stepped frequency ground penetrating radar[J]. Radar Science and Technology, 2025, 23(3): 253-261.

[22]	Luo T X H, Lai W W L, Giannopoulos A. Forward modelling on GPR responses of subsurface air voids[J]. Tunnelling and Underground Space Technology, 2020, 103: 103521.

[23]	Warren C, Giannopoulos A, Giannakis I. gprMax: open source software to simulate electromagnetic wave propagation for ground penetrating radar[J]. Computer Physics Communications, 2016, 209: 163-170.

[24]	Warren C, Giannopoulos A, Gray A, et al. A CUDA-based GPU engine for gprMax: open source FDTD electromagnetic simulation software[J]. Computer Physics Communications, 2019, 237: 208-218.

[25]	Smith S L, O’Neill H B, Isaksen K, et al. The changing thermal state of permafrost[J]. Nature Reviews Earth & Environment, 2022, 3(1): 10-23.

[26]	郭昰龙, 鹿琪, 刘四新, 等. 基于土壤分形的探地雷达三维随机介质建模[J]. 地球物理学进展, 2025, 40(5): 2201-2210.

[27]	Guo Xialong, Lu Qi, Liu Sixin, et al. 3D random medium modeling of ground penetrating radar based on soil fractal[J]. Progress in Geophysics, 2025, 40(5): 2201-2210.

[28]	张清河, 吴欣悦, 刘含, 等. 一种新颖的探地雷达快速正演模拟及埋地目标探测机器学习方法[J]. 地球物理学报, 2023, 66(8): 3482-3492.

[29]	Zhang Qinghe, Wu Xinyue, Liu Han, et al. A novel fast forward modeling simulation of ground penetrating radar and machine learning method for buried target detection[J]. Chinese Journal of Geophysics, 2023, 66(8): 3482-3492.

[30]	Qu D X, Luo Y, Li X P, et al. Study on the stability of rock slope under the coupling of stress field, seepage field, temperature field and chemical field[J]. Arabian Journal for Science and Engineering, 2020, 45(10): 8315-8329.

[31]	Li G F, Li N, Bai Y, et al. A novel simple practical thermal-hydraulic-mechanical (THM) coupling model with water-ice phase change[J]. Computers and Geotechnics, 2020, 118: 103357.

[32]	Han Z H, Li R R, Zhang L Q, et al. Multi-field coupling simulation of impact of temperature and density of heat injection well on carbon budget during hydrate exploitation in Qilian Mountain Permafrost Region[J]. Journal of Earth Science, 2024, 35(6): 1934-1943.

[33]	白青波, 李旭, 田亚护, 等. 冻土水热耦合方程及数值模拟研究[J]. 岩土工程学报, 2015, 37(S2): 131-136.

[34]	Bai Qingbo, Li Xu, Tian Yahu, et al. Hydro-thermal coupling equation of frozen soil and its numerical simulation[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(S2): 131-136.

[35]	Zhang S, Teng J D, He Z Y, et al. Importance of vapor flow in unsaturated freezing soil: a numerical study[J]. Cold Regions Science and Technology, 2016, 126: 1-9.

[36]	Li X K, Li X, Liu J K. A dynamic soil freezing characteristic curve model for frozen soil[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2024, 16(8): 3339-3352.

[37]	Zhan Y X, Zhao M Y, Lu Z, et al. Study on water-heat coupling migration law and frost heave effect of soil slope in seasonal frozen regions during groundwater recharge[J]. Environmental Earth Sciences, 2023, 82(17): 401.

[38]	王东源, 张科, 陈熹, 等. 一种冻土温度场方程的解析解[J]. 地震工程学报, 2024, 46(6): 1373-1379.

[39]	Wang Dongyuan, Zhang Ke, Chen Xi, et al. An analytical solution for the temperature field equation of frozen soil[J]. China Earthquake Engineering Journal, 2024, 46(6): 1373-1379.

[40]	谭贤君, 褚以惇, 陈卫忠, 等. 考虑冻融影响的岩土类材料导热系数计算新方法[J]. 岩土力学, 2010, 31(S2): 70-74.

[41]	Tan Xianjun, Chu Yidun, Chen Weizhong, et al. A new method for calculating thermal conductivity of geotechnical materials considering freeze-thaw effect[J]. Rock and Soil Mechanics, 2010, 31(S2): 70-74.

[42]	郝君明, 贾佩钱, 李旺平. 多年冻土斜坡稳定性的模拟分析[J]. 自然灾害学报, 2024, 33(2): 152-163.

[43]	Hao Junming, Jia Peiqian, Li Wangping. Simulation analysis of stability of permafrost slopes[J]. Journal of Natural Disasters, 2024, 33(2): 152-163.

[44]	祁子寒, 王云琦, 李通, 等. 缙云山典型林地根系分布对土壤土水特征和渗透性的影响[J]. 中国水土保持科学(中英文), 2024, 22(4): 162-170.

[45]	Qi Zihan, Wang Yunqi, Li Tong, et al. Effects of root distribution on soil-water characteristics and permeability of typical forest lands in Jinyun Mountain[J]. Science of Soil and Water Conservation, 2024, 22(4): 162-170.

[46]	闫相君, 王云琦, 祁子寒, 等. 火烧迹地根系分布对浅层含根土土水特征及渗透特性的影响[J]. 中国水土保持科学(中英文), 2023, 21(3): 42-51.

[47]	Yan Xiangjun, Wang Yunqi, Qi Zihan, et al. Effects of root distribution on soil-water characteristics and permeability of shallow root-bearing soil in burned areas[J]. Science of Soil and Water Conservation, 2023, 21(3): 42-51.

[48]	张经纬, 刘新, 叶盛波, 等. 一种用于准单站配置的步进频探地雷达建模方法[J]. 雷达学报, 2019, 8(4): 519-526.

[49]	Zhang Jingwei, Liu Xin, Ye Shengbo, et al. A modeling method of stepped-frequency ground penetrating radar for quasi-monostatic configuration[J]. Journal of Radars, 2019, 8(4): 519-526.

[50]	陶春康, 汪海勇, 瞿诗华. 步进扫频宽带探地雷达信号处理[J]. 微波学报, 2020, 36(2): 67-70, 78.

[51]	Tao Chunkang, Wang Haiyong, Qu Shihua. Signal processing of stepped frequency continuous wave ground penetrating radar[J]. Journal of Microwaves, 2020, 36(2): 67-70, 78.

[52]	Liu S, Cao Y H, Yeo T S, et al. Range sidelobe suppression for randomized stepped-frequency chirp radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2021, 57(6): 3874-3885.

[53]	Zhao L, Zou D F, Hu G J, et al. A synthesis dataset of permafrost thermal state for the Qinghai-Tibet (Xizang) Plateau, China[J]. Earth System Science Data, 2021, 13(8): 4207-4218.

[54]	赵丹丹, 刘家福, 刘吉平, 等. 东北寒区多年冻土退化的主要研究综述[J]. 科学技术与工程, 2022, 22(30): 13151-13161.

[55]	Zhao Dandan, Liu Jiafu, Liu Jiping, et al. Review on degradation of permafrost in cold regions of Northeast China[J]. Science Technology and Engineering, 2022, 22(30): 13151-13161.

基金资助

国家重点研发计划(2023YFB2604801)

泰山学者工程项目

PDF (12799KB)

访问

被引

详细

导航

Received	Published	Revised
2026-01-06	2026-09-08	2026-02-03
Online Date
2026-09-18

摘要