基于改进剥层法的南极冰盖密度反演算法

杨望笑 窦银科* 稂时楠 赵 博 王煜尘 左广宇

①(太原理工大学电气与动力工程学院 太原 030024)

②(北京工业大学信息与通信工程学院 北京 100124)

③(中国科学院空天信息创新研究院 北京 100190)

1 引言

调频连续波(Frequency Modulated Continuous Wave, FMCW)雷达技术目前已在测距、气象监测、极地冰雪探测等领域得到广泛应用。南极冰盖密度对于估算冰盖表面物质平衡是至关重要的[1,2]，文献[3]阐述了冰盖密度与介电常数之间的经验公式，因此，冰盖密度可经雷达反演冰盖垂直剖面上的介电常数转换得到。

目前通过雷达对冰盖密度/介电常数进行反演的主要方法是基于共中心点的雷达反射信号。发射和接收天线之间的距离需在0～80 m范围内不断变化，提取某一反射层在不同天线间距下的传播时间差，再通过Dix反演、相似性分析、干涉测量等方法估算冰盖表面至反射层之间的等效介电常数，经文献[3]中的经验公式转换为平均密度[4,5]。由于这种方法数据采集过程较为繁琐，在进行连续的介电常数剖面反演时使用较少。

当采用单条雷达反射信号进行介电常数反演时，车辆或飞机在采集数据过程中可同时搭载发射和接收天线，高效获取冰盖内部结构和介电常数剖面，这对大范围了解冰盖物质平衡是至关重要的。但是由于少量液态水会导致雷达信号强烈的反射和衰减，因此基于单条反射信号的反演法无法在湿雪带进行介电常数剖面反演，更适用于高海拔、低积聚率地区，且在反演前需分析研究区域中冰的介电特性以确定该类方法的可行性[6]。目前使用单条雷达反射信号反演介电常数的主要方法是高斯牛顿法、退火参数算法或粒子群算法寻找介电常数的最优解[7,8]，但是这些方法的计算过程相对复杂，更适用于层位尺寸较大、数量少的应用条件。实际上冰盖是由无数个薄层组成，每个层之间都会存在微小的介电常数差异，因此并不适用于传统寻找最优解的方法。剥层反演法可提取反射信号峰值的幅度与时延，当已知第1个层位的介电常数时，可根据幅度估算下一层的介电常数，进而递推求出整个介电常数剖面。剥层反演法在应用过程中还存在方法误差，导致反演结果在局部深度范围内存在误差[6]。因此有必要对基于单条反射信号的介电剖面反演算法进行改进。

为了弥补上述算法的不足，本文提出一种基于剥层反演法和密实化公式相结合的优化算法。采用快速傅里叶变换(Fast Fourier Transform, FFT)对差频信号进行频谱估计后，首先使用剥层反演法对回波信号进行初步的密度反演，再根据冰盖密实化机理进行校准，基于校准结果仿真差频信号，通过与原始回波信号进行对比寻找校准方程的最优解。对理论冰盖模型和实测雷达数据的实验结果表明，此算法可有效提高反演精度。

2 FMCW冰雷达工作模型

冰盖是空气与冰的混合物，可将其内部结构假设为L个平行且有介电特性差异的冰层。FMCW雷达差频信号的时域表达式为

以6.25 MS/s的采样率对差频时域信号采样，获取N=25000个采样点。为降低频谱估计误差，在N个采样点之后填充40×N个0，再进行FFT。回波信号由L–1个峰值组成，如图1所示，其中每个峰值表示层位交界面产生的反射。

图1 FMCW雷达工作示意图

回波信号中第k个峰值的位置表示在深度dk处交界面反射信号的时延τk，该峰值幅度AMP的表达式为

3 反演算法

本文算法是基于剥层反演法，并在此基础上进行校准的一种优化算法，下面分别对算法的应用方法和优化过程进行介绍。

3.1 剥层反演法在冰盖上的应用

由式(3)—式(8)可知，每个峰值的幅度由层位间介电常数差异决定。反演时可以提取回波信号中的峰值，当已知某一个层的介电常数εk即可根据差频信号幅度和式(3)、式(4)、式(8)估算出下一层的介电常数εk+1，递推求出整个垂直剖面上的介电常数[6]。但是该方法更适用于层位厚度明显大于雷达分辨率的情况[10]。因此，对传统的剥层反演法进行改进十分必要。

1927年的农村调查，也折射了读书人在村民眼中的尴尬处境：“读书成本太大，出来非但没有官做，即教员位置亦粥少僧多，而况学些空架子，文不象秀才，武不象丁，手不能提篮，肩不能挑担，不事生产，要吃要用。”[21]

文献[11]分别公布了冰芯B31的密度和介电常数测量数据。图2(a)黑线揭示了冰芯5 mm深度分辨率的相对介电常数数据，对30 m附近的数据放大后表明相邻的层位也会表现出介电常数差异，这意味着差频信号中峰值所表示的信号反射并不是由一个冰层的介电特性突变造成的，而是无数个冰层之间微小的介电常数变化引起的许多微弱反射相互干涉形成的[12]。除此之外，当电磁波向由低向高介电常数，或由高向低介电常数的冰层传播时，都会产生电磁波的反射，因此，直接使用差频信号峰值进行剥层反演会产生明显的反演误差。

综上，在进行剥层反演之前，需要通过对整个回波信号的峰值进行光滑样条回归并归一化至[0,1]，可以忽视局部深度范围内冰层介电常数波动产生的反射，从而提取出整个垂直剖面上的幅值变化趋势。根据式(1)可知，第1层的差频信号幅值仅取决于反射系数，所以可以使用冰盖实测介电常数估算冰盖表面反射系数，并将其用于幅值变化趋势的反归一化变换。通过上述步骤可将回波信号转换为反射系数曲线，离散化后再使用剥层反演法，递推求出整个介电常数剖面，进而转化为密度剖面。

文献[3]公布的密度ρ与相对介电常数ε之间的经验公式被广泛应用于冰盖密度反演[4,5]，表达式为

通过式(10)将图2(a)黑线所表示的相对介电常数转化为密度，再分别使用宽度为4 m的窗口对实测密度(图2(b)红线)和由介电常数转化的密度(图2(b)蓝线)进行滑动平均，两者结果吻合，证明了式(10)的有效性。

图2 冰芯B31密度和相对介电常数

实际上差频信号峰值反映的是冰盖内部介电常数标准差而不是密度增长率[13]，后者是反应密度随深度增长过程的重要指标，而剥层反演法将两者视为等效替代。通过对目前已公开的冰芯数据进行统计分析，两者呈现相似的变化特征。以冰芯B31为例，分别使用宽度为0.7 m的窗口计算密度-深度变化率和介电常数标准差，并归一化至[0, 1]后使用光滑样条回归提取变化趋势，结果如图3所示。可以发现两者呈现出相似变化趋势，即在0～30 m内大幅度降低，在30～70 m内升高后，紧接着又逐渐降低。但是在局部区域内又存在差异，介电常数标准差在50 m以后明显高于密度变化率，这会导致50 m以后的剥层反演结果被高估。这属于方法误差，因此，对剥层反演得到的初步反演结果进行校准十分必要。

图3 冰芯B31的介电常数标准差和密度变化率归一化

3.2 密实化公式

初步反演结果中每个冰层的深度和密度均需要校准，其中包含了大量的未知参数，直接使用模拟退火算法、高斯牛顿法等全局优化算法，计算过程过于复杂[14]。本文提出利用冰盖密实化公式来表示密度变化过程，可有效降低未知参数的数量以提高计算速度，同时还可基于冰盖密实化原理约束剥层反演得到的初步反演密度，从而达到数据校准的目的。

目前用于表示密度-深度关系的理论主要是基于密实化原理的半经验公式，并且需要调节参数将理论密度剖面与实际冰芯数据相匹配。本文基于文献[15]的密实化原理，对文献[16]提出的密度-深度公式进行了优化，表达式为其中，Ec,Eg和C是常数，g是重力加速度，ρi是纯冰密度 (917 kg/m3)，Ts是平均冰盖表面温度，R是气体常数，dz是冰层的厚度，dρ是相较于上一冰层的密度增长量。文献[16]根据对南极Jurassic和Dyer 高原统计的密度特征对常数C进行拟合，结果表明在密度-深度变化率的第1个极小值点之前为0.03，而之后为0.07。

根据剥层反演法原理可知，方法误差会随着深度不断的增加，这意味着浅层的反演误差相对较小，因此可以用剥层反演法得到的初步反演结果确定极小值之前的dρ1，再使用优化算法寻找出常数A的最优解来确定dρ2，最终对dρ进行积分求出校准后的密度剖面。该方法在未依托式(11)中环境参数及经验参数的情况下，利用密实化原理完成了数据校准，拓展了本文反演算法的适用范围。

3.3 优化算法

剥层反演法可以得到初步反演密度和密度-深度变化率极小值点所在的深度。通过指数函数对极小值之前的dρ1进行拟合，由于优化过程中的计算量较小，本文通过穷举法列举常数A在[1, 15]范围内可能产生的密度剖面，并在其中寻找校准密度的最优解。

由3.1节可知，回波信号被视为密度增长趋势，通过对回波信号进行积分可以近似的表达密度增长量。在图3中100 m处密度增长率和介电常数标准差的积分结果相近，这表明虽然初步反演密度与冰芯实测密度之间在局部深度范围内有所差异，但是两者在100 m处的密度相近，因此可以将回波信号在0～100 m内的积分差异作为校准密度最优解的判定标准。本文将对每个校准密度的回波信号进行仿真，并将其与雷达原始回波信号的积分相对比，当差异最小时则说明校准密度不仅反应了密度增长趋势，还校正了局部深度范围内的密度误差。

为减少计算时间，本文将每个校准密度固定间隔离散化后基于FMCW雷达原理仿真差频信号。首先根据式(3)、式(4)、式(8)确定每一层的反射系数、传输系数和衰减量，再使用式(9)计算每个层的时延，将上述参数代入式(1)求出每个冰层的反射信号，对其累加后得到差频时域信号，使用FFT进行频谱估计后提取峰值，对模拟的回波信号峰值进行光滑样条回归，归一化至[0, 1]后进行积分，分别将每个校准密度的积分结果与雷达原始回波信号的积分结果相减后取绝对值，取差异最小的密度曲线作为校准密度的最优解。综上所述，冰盖垂直剖面上密度和介电常数曲线的反演和校准流程如图4所示。

图4 冰盖密度反演和校准流程图

4 算法应用

为验证本文所提出算法的有效性，这里给出一个反演实例，将优化算法应用于冰盖模型和雷达实测数据中。

4.1 理论模型实验

本文根据图2(a)冰芯B31在0～100 m内的实测介电常数-深度数据来建立1维水平层状冰盖模型[11]，整个垂直剖面上均使用平均电导率23.16 μS/m，且忽略冰的各向异性[17]。基于时域有限差分方法(Finite Difference Time Domain, FDTD)对差频时域信号进行仿真，频谱估计后使用反演算法估算密度剖面，并与图2(b)中的冰芯实测密度对比。

根据在南极实际使用的FMCW雷达系统参数设置激励源，电磁波发射频率为0.5～2 GHz，扫频周期为4 ms；吸收边界条件为完全匹配层(Perfectly Matched Layer, PML)，吸收层宽度为40；图2(a)介电常数数据的深度分辨率为5 mm，本文设置空间步长dx=2.5 mm；根据Courant稳定性条件计算得到时间步长dt=8.65×10–12s。仿真结束后，将接收信号与发射信号混频，求出差频时域信号，通过FFT获得的回波信号如图5(a)黑线所示。

使用1.5×10-3的表面反射系数和1.61的表面相对介电常数(320 kg/m3的表面密度)，得到初步密度剖面如图5(b)中蓝线所示。密度-深度变化率极小值点位于深度20 m。通过与窗口宽度为4 m的滑动平均密度(图5(b)绿线)进行对比，可以明显观察到50～90 m深度范围内的反演密度被高估，这是由于该深度范围内介电常数标准差大于密度变化率(图3)，导致均方根误差(Root Mean Square Error,RMSE)为2.49%。

使用指数函数对0～20 m内的密度进行拟合，以0.1的间隔分别计算常数A在[1, 15]范围内的密度剖面，并基于式(1)对密度剖面的差频时域信号进行仿真，频谱估计后进行光滑样条回归并归一化后积分，如图5(c)展示了每个回波信号与原始回波信号的积分差异，其中最小值为A=3.7，所代表的密度剖面如图5(b)红线所示。寻找到的最优密度曲线的回波信号如图5(d)红线所示，300～800 ns处的回波信号更低，且50～90 m深度范围内的反演密度被校正，RMSE降低至1.39%，与冰芯实测密度的相关性为0.999。

图5 冰芯B31的密度反演结果

4.2 雷达数据实验

中国第32次南极考察队使用FMCW雷达探测了南极冰盖，本文选择了冰芯LGB69附近的雷达数据来反演平均密度剖面，并与冰芯密度进行对比，结果如图6。

对雷达数据进行反演之前，至少需要使用40道差频信号进行相干平均以减少非相干噪声。使用表面相对介电常数1.82，即415 kg/m3的表面密度，得到初步密度剖面如图6中黑线所示，其中30～60 m内的密度被低估，而70 m之后的密度被明显的高估。当A=1.1时的校准密度为最优解，如图6中红线所示，与冰芯密度(图6绿线)的RMSE从3.86%降低至2.13%。上述验证结果表明，本文的反演方法可以有效提升反演精度，校正了方法误差。尽管密度RMSE仅降低了约1%，但是这对于研究冰盖的密实化率有着至关重要的作用。

图6 冰芯LGB69附近的密度反演结果

表面密度或表面介电常数决定了反演的平均密度，因此本文算法依赖表面密度数据。由于在反演算法中对回波信号进行了归一化处理，当反演不同带宽的FMCW雷达数据时，差频信号之间的绝对电平差异被消除，因此反演结果仍能反映冰盖的平均密度变化趋势。考虑到噪声的存在，对实际雷达数据进行反演前相干平均是必不可少的，对于不同信噪比的雷达系统，可以在相干平均时适当调节雷达数据的道数以保证反演结果的准确性。

5 结束语

本文基于FMCW雷达差频信号在层状介质中的传播特性和冰盖的密实化原理，提出了一种将剥层反演法与密实化公式相结合的冰盖密度反演方法。通过剥层反演法初步得到反演结果，根据冰盖密实化规律对初步反演结果进行校准。理论模型和实际测量结果表明，该方法可有效提高反演精度，并校正了剥层反演法中存在的方法误差，即介电常数标准差与密度变化率之间存在的差异。国内外学者已经获取了大量的探冰雷达数据，基于本文的反演算法建立数据分析系统可以对雷达数据进行更充分的挖掘与利用，进而获取整个南极的密度剖面，极大提升冰盖物质平衡估算精度。