一种基于短时傅里叶变换的极地冰下水体探测方法

摘要

本发明涉及一种基于短时傅里叶变换的极地冰下水体探测方法，包括以下步骤：1)获取无线电回波测深数据；2)获取航向模糊处理后的二维雷达图像；3)对航向模糊处理后的二维雷达图像进行裁切拼接处理，得到裁切后的二维雷达图像；4)对裁切后的二维雷达图像中的分界面信号进行波形重整；5)对波形重整后的二维雷达图像进行短时傅里叶变换处理，获取每一道测线在分界面处的频率响应及幅值；6)结合归一化后的频率响应、短时傅里叶变换后频率响应幅度以及地形坡度计算冰底物质探测值，进行冰下水体的自动识别。与现有技术相比，本发明以冰底物质探测值作为判断冰底物质是否为冰下湖和冰下水体的依据，实现单道雷达信号冰底物质的高精度分辨。

说明书

技术领域

本发明涉及地球气候和地球动力学技术领域，尤其是涉及一种基于短时傅里叶变换的极地冰下水体探测方法。

背景技术

极地冰川是地球气候和地球动力学系统的重要组成部分，长久以来仍然缺乏来自南极冰盖底部足够丰富的地质和地球物理数据，难以支撑对南极历史变迁的认识和对全球气候变化发展的预估。为解决这一问题，世界各国科研团队通过航空地球物理测绘，收集了包括重力，磁力，测高，冰雷达等南极地区的地球物理数据，对于了解冰盖平衡的潜在变化起到了非常重要的作用，其中，冰下湖是南极冰底环境的重要组成部分之一，它对冰盖物质平衡、冰盖稳定性、全球气候变化和海平面变化有着重要的影响，截止到2022年2月，全球共发现了773个冰下湖，其中，南极洲有675个冰下湖，格陵兰岛有64个。

目前探测冰盖及冰下的技术有很多，然而无线电回波测深技术对冰盖内部及底部的探测是最为敏感的，弥补了深冰芯钻探在恶劣的南极地区无法大规模实施的缺陷，同时也克服了地震、重力勘探等地球物理探测技术对于冰盖不敏感的不足，为南极冰下水体探测提供了重要的数据源。在寻找冰下水体的过程中，最常用的方法是目视解译，然而它的局限性也很明显，例如，它通常很费时，这对未来大量无线电回波测深数据的解译极为不利，而且解译标准随着解译人员而变化，难以统一，因此开发自动的冰下水体探测方法对于极地研究而言是必不可少的。

传统的冰下湖探测方法通常依赖于二维冰雷达信号(B-Scan)数据的统计信息，如回波强度、反射率、粗糙度和水势等，而忽略了一维冰雷达信号(A-Scan) 的固有时频特性，缺乏冰下水体探识别指标导致降低极地冰下水体探测的精度。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于短时傅里叶变换的极地冰下水体探测方法，通过利用基底信号频率响应特征对极地地区的无线电回波测深数据进行处理，构建一个可以描述基底信号形态特征的频率指标，并联合坡度、幅值等指标获得基底存在水的可能性的极地专题信息。

本发明的目的可以通过以下技术方案来实现：

一种基于短时傅里叶变换的极地冰下水体探测方法，包括以下步骤：

1)获取无线电回波测深数据；

2)对原始二维雷达图像进行沿轨道方向的数据模糊处理，获取航向模糊处理后的二维雷达图像；

3)对航向模糊处理后的二维雷达图像进行裁切拼接处理，得到裁切后的二维雷达图像；

4)对裁切后的二维雷达图像中的分界面信号进行波形重整；

5)对波形重整后的二维雷达图像进行短时傅里叶变换处理，获取每一道测线在分界面处的频率响应及幅值；

6)结合归一化后的频率响应、短时傅里叶变换后频率响应幅度以及地形坡度计算冰底物质探测值，进行冰下水体的自动识别。

所述的步骤1)中，无线电回波测深数据具体为冰盖遥感中心提供的不同相干雷达测深仪数据。

所述的步骤2)中，航向模糊处理后的二维雷达图像的表达式为：

其中，BScan表示原始二维雷达图像，w

所述的步骤3)中，进行数据裁切的表达式为：

其中，BScan″表示裁切后的二维雷达图像，

所述的步骤4)具体包括以下步骤：

41)对裁切后的二维雷达图像中的每一道大小为2w

其中，AScan表示裁切后的二维雷达图像中任意一条A-Scan，2w

42)对去直流漂移处理后的二维雷达图像中的分界面波形进行波形重整。

所述的步骤42)中，波形重整具体包括：

421)设定主峰的左右端点，并翻转主峰后分别拼接到主峰左右端点两侧；

422)将其他采样点的值置为0，则有：

其中，

所述的主峰的左右端点的定义为：

将分界面峰值强度的六分之一处设定为阈值，阈值与A-Scan波峰左右两侧最接近波峰位置的两个交点分别作为分界面主峰的左右端点。

所述的步骤5)中，短时傅里叶变换公式为：

其中，x(τ)为待处理的信号，w(*)为窗函数，在t时刻做短时傅里叶变换处理的窗函数为w(τ-t)，ω为角速度。

所述的步骤6)中，冰底物质探测值C的计算式为：

其中，

所述的步骤6)中，进行冰下水体的自动识别具体为：

将冰底物质探测值C作为最终判断冰底物质是否为冰下湖和冰下水体的依据，若

与现有技术相比，本发明具有以下优点：

在现有冰下湖探测技术中，往往依赖连续多道回波强度的绝对大小和相对大小，辅以地形、水势、粗糙度等参数来检测，但是这种检测方法忽略了A-Scan波形的固有特征，缺乏对单道雷达信号的检测能力和解释能力，本发明提出的基于短时傅里叶变换的波形重整方法，提取了额外的特征，以冰底物质探测值作为最终判断冰底物质是否为冰下湖和冰下水体的依据，实现了对单道雷达信号进行冰底物质进行高精度的分辨。

附图说明

图1为基于短时傅里叶变换的极地冰下湖探测流程图。

图2a为二维雷达剖面图像。

图2b为不同分界面的波形图。

图3a为航向模糊后的二维雷达剖面图，曲线表示距冰底分界面两侧各150采样点的位置。

图3b为数据裁切后的二维雷达剖面图。

图4为数据裁切后的A-Scan图像，其中，图(4a)为冰水分界面，图(4b) 冰岩分界面，左侧为波浪线的曲线表示去直流漂移后的波形，左侧为直线的曲线表示波形重整后的波形。

图5a为冰水分界面的短时傅里叶变换结果，图中，(a)表示冰水分界面 A-Scan位置，(b)表示冰水分界面A-Scan波形重整，(c)表示对图(b)做短时傅里叶变换的结果，圆形位置为频率响应。

图5b为冰岩分界面的短时傅里叶变换结果，图中，(d)表示冰岩分界面 A-Scan位置，(e)表示冰岩分界面A-Scan波形重整，(f)表示对图(e)做短时傅里叶变换的结果，圆形位置为频率响应。

图6为测线20081225_02_024所对应的冰雷达图像及探测结果，其中，图(6a) 为二维雷达图像及探测结果(颜色带)，图(6b)为短时傅里叶变换后的频率响应及频率对应幅值，图(6c)为冰底地形坡度及冰底分界面处的回波强度，图(6d) 为冰厚及水势，在图(6a)中，深色半透明矩形为Livingstone等人的冰下湖的边界，浅色半透明矩形为Wolovick等人的冰下湖的边界。

图7为AGAP区域不同雷达剖面图的冰下湖及冰下水体探测结果，其中，图 (7a)为测线20081228_03_006，图(7b)为测线20081225_02_024，图(7c)为测线20090109_01_023，图(7d)为测线20090107_03_015及测线20090107_03_016，图(7e)为测线20081223_01_015，图(7f)为测线20090102_03_001，图(7g) 为测线20090108_01_016及测线20090108_01_017，图(7h)为测线20090106_05_013，图(7i)为测线20090108_03_021，图(7j)为测线20081228_01_020，图(7k) 为测线20081226_02_017，图(7l)为测线20081225_02_018，图(7m)为测线 20090106_05_020及测线20090106_05_021，图(7n)为测线20090107_02_014、测线20090107_02_015、测线20090107_02_016、测线20090107_02_017及测线20090107_02_018，深色半透明矩形为Livingstone等人2022年发表的在CReSIS 测线上的冰下湖的边界，浅色半透明矩形为Wolovick等人2013年发表的在CReSIS 测线上的冰下湖的边界。

图8为2009_Antarctica_TO_Gambit目录的探测结果，其中，AGAP区域位于正中矩形测线区域，灰色圆点为Livingstone等人2022年发表的在CReSIS测线上的冰下湖，白色圆点为Wolovick等人2013年发表的在CReSIS测线上的冰下湖，深色圆点为本发明方法所计算出的冰下湖，圆点大小代表此处为冰下湖或冰下水体的可能性的大小。

图9为CReSIS南极区域数据集探测结果，其中，蓝色实心圆点为本发明方法探测出的冰下湖及冰下水体，红色实心圆点为Livingstone等人总结出的南极地区 675个冰下湖，白色实线为冰雷达测线，图中，(a)表示全南极探测结果，(b)表示AGAP区域探测结果，(c)表示Vostok湖区域探测结果，(d)表示Byrd冰川区域探测结果，(e)表示西南极探测结果，(f)表示罗斯冰架西侧探测结果，(g)表示Coats Land区域探测结果。

图10为本发明探测结果与水势预测结果的对比，其中，(a)表示全南极探测结果对比，(b)表示AGAP区域探测结果对比，(c)表示Vostok湖区域探测结果对比，(d)表示Byrd冰川区域探测结果对比，(e)表示西南极探测结果对比，(f) 表示罗斯冰架西侧探测结果对比，(g)表示Coats Land区域探测结果对比。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

如图1所示，本发明旨在针对冰下水体探测方法中现有冰下水体探识别指标的缺乏，提供一种基于短时傅里叶变换(Short-Time Fourier Transform，STFT)的极地冰下水体探测方法，基于冰架遥感中心(CReSIS)项目整理公开的冰雷达数据，编写了MATLAB算法程序，对全南极和全格陵兰范围内的机载冰雷达信号进行了处理，对冰底物质进行基岩-水的二分类，找出并绘制测区内所有可能存在冰下水体及冰下通道的地理位置。

当电磁波到达冰水界面时，由于水极高的介电常数(∈w≈81)和其光滑的表面，大部分电磁波被反射，只有极少的能量穿过界面进入水中，从而在A-Scan上产生一个急剧变化的特征，相反，岩石或其他非水材料的介电常数并不高(如岩石，∈r≈4～9)，更多的电磁波能进入冰岩界面，并在岩层内由于物质不均匀性产生多次反射，同时粗糙的冰岩界面产生的多次反射会相干相消，这使得A-Scan中分界面处的回波强度逐渐下降。

冰盖下方的电磁波相对功率的变化只取决于反射率的变化，根据反射率计算公式(1)，分界面两侧不同物质的介电常数相差越大，则反射率越高，回波功率变化得越剧烈，A-Scan上的分界面处的特征也就越狭窄且尖锐，图2b展示了不同物质分界面处的A-Scan波形(上方黑线为冰岩界面，下方黑线为冰水界面)。

其中，R

由上可知，雷达信号在冰岩界面和冰水界面的反射系数不同，冰水与冰岩界面的主峰也会显示出不同的特征——冰水界面的主峰狭窄而尖锐，而冰岩界面的主峰较宽，这使得冰下湖泊在雷达图上显示出明显的高频特征，而冰下基岩则显示出低频特征，在现有的冰下湖探测中，往往依赖于连续多道回波强度的绝对大小和相对大小，辅以地形、水势、粗糙度等参数来检测，而这种检测方法忽略了A-Scan波形的固有特征，缺乏对单道雷达信号的检测能力和解释能力，基于此，本发明提出了基于短时傅里叶变换的冰下水体探测方法，用以实现通过单道雷达信号进行冰底物质分辨，具体实施步骤如下：

步骤一、无线电回波测深数据准备；

所用无线电回波测深数据为冰盖遥感中心(CReSIS)提供的不同相干雷达测深仪的数据产品，可免费下载，下载产品数据以测线的覆盖范围为数据位置，提供了包括测线上的经纬度、飞机飞行高度、飞机至冰床的距离、冰厚、冰表面高程、冰床高程、雷达回波图、回波图上冰表面和冰床的位置等。

步骤二、通过对B-Scan数据进行沿轨道方向的数据模糊方法，对原始数据进行降噪和特征增强；

机载雷达测量的每个A-Scan都是相互独立的，但相邻的A-Scan之间又往往存在着一些共同点：

1)地理位置相近的分界面，往往海拔高度也相近，即A-Scan波形中分界面的峰值都位于相近的采样点处；

2)每一道A-Scan中的噪声都是随机噪声。

因此，通过对B-Scan进行航向模糊处理，即对连续的几道A-Scan进行取平均值的操作，突出A-Scan中的相似部分，并在不改变信号最强波峰的固有频率的情况下过滤掉随机噪声，同时，由于冰下湖大部分是水平的，而基岩往往有一定坡度，这会使得重构后的A-Scan冰水分界面的峰更加突出，而冰岩处的峰更加平缓，使得能够进一步通过频率响应来对二者进行区分，航向模糊的公式为：

其中，BScan表示原始二维雷达图像，w

步骤三、对航向模糊处理后的二维雷达图像进行裁剪拼接；

一个完整的A-Scan通常包含直达波、表面反射波、等时层、冰底反射波(基岩或水体)，有些还会包含冰流造成的等时层(internal layer,flow disrupted layer)、冰底融化后重新冻结的冰(frozen-on ice)和雷达产生的类似条形码的干扰条纹等，通常不需要处理从冰盖表面到冰盖底部的整个A-Scan，只需要处理冰底分界面附近的一些采样点即可，本发明采用CReSIS提取的分界面，但选取分界面上下50个采样点范围内的最大值位置作为真实的分界面位置，然后在新的分界面处上下各取个采样点，形成一个新的条带状B-Scan，此时B-Scan内通常仅包含冰岩、冰水分界面及少量等时层，数据裁切的公式为：

其中，BScan″为裁切后的二维雷达图像，

步骤四、对分界面信号进行波形重整；

每个采样点所对应的绝对数值是没有物理意义的，基岩和水的回波强度差距体现在范围内相对数值的变化。此外，对A-Scan进行傅里叶变换可以发现，一段信号的平均值的绝对值较大将会导致其频谱中零频分量非常强，因此，首先对裁切后的B-Scan中的每一道大小为2w

其中，AScan表示裁切后的二维雷达图像中任意一条A-Scan，2w

将直流漂移处理后的分界面峰值强度的六分之一处设定为阈值，阈值与 A-Scan波峰左右两侧最接近波峰位置的两个交点作为分界面主峰的左右端点，如图(4a)和图(4b)中的A、B两点，主峰强度超过阈值的部分为峰高，左右端点的距离为峰宽。

其中，l为主峰与阈值的左交点对应的采样点号，r为主峰与阈值的右交点对应的采样点号，Peak为主峰在阈值之上的部分，

图5中A、B两点间的信号的频率值即可作为分界面信号波形的固有频率值，由于主峰部分本体的均值依然较大，故直接对主峰进行傅里叶变换，其零频分量所占的比重很高，导致其频率响应为0MHz，基于此，本发明对分界面波形进行了波形重整，具体操作如下：

(1)翻转主峰，拼接到主峰左右端点两侧；

(2)将其他采样点的值置为0。

具体公式为：

其中，

波形重整前后对比图如图4所示。

步骤五、通过短时傅里叶变换求出每一道测线在分界面处的频率响应以及幅值，对波形重整后的信号进行短时傅里叶变换处理，公式为：

其中，x(τ)为待处理的信号，w(*)为窗函数，在t时刻做短时傅里叶变换处理的窗函数为w(τ-t)，窗内截取的短时信号x

通过短时傅里叶变换方法的分析，可以得到每个目标A-Scan的二维时域-频域矩阵，并可以总结出冰水界面和冰岩界面的不同特征，不同物质分界面的短时傅里叶变换分析的结果如图5所示，对重构后的波形做短时傅里叶变换，可以看到，重构后的冰水分界面的主峰频率响应在1MHz左右，而冰岩分界面的主峰频率响应依旧在0MHz左右，因此，可以根据频率响应对基岩和冰下湖进行最初步的区分。

步骤六、结合归一化频率值、频率幅值、冰底坡度等多参数，计算出冰底物质的探测值，进行冰下水体的自动识别。

有些冰下含水的沉积物，其频率响应并不为零，为了进一步区分不同物质的分界面，将短时傅里叶变换的结果作为参数值输入到最终判断函数中：

其中，

其中，Δz为相邻两道A-Scan主峰的高程差，Δx为相邻两道A-Scan的水平距离，单位均为m。

本例中，C为平滑前的冰底物质探测值，是一个向量，C′为平滑后的冰底物质探测值，是一个向量；

本发明将冰底物质探测值C作为最终判断冰底物质是否为冰下湖和冰下水体的依据，若

实施例

在CReSIS数据集中，AGAP区域的测线位于2009_Antarctica_TO_Gambit目录下，该目录中还包括一些向Dome Fuji和Vostok方向延伸的测线，以及一条向 Coats Land方向延伸的测线。目录下的有效B-Scan共有889段，有效A-Scan共有 1471608条。同一条测线上相连的B-Scan与B-Scan之间用黑白两色的粗线段进行区分。

AGAP区域的探测结果如图8所示，灰色圆点为Livingstone等人2022年发表的在CReSIS测线上的冰下湖；白色圆点为Wolovick等人2013年发表的在CReSIS 测线上的冰下湖；深色圆点为本发明方法所计算出的冰下湖，圆点大小代表此处为冰下湖或冰下水体的可能性的大小。

图8中白色线段为图7中典型B-Scan例子所在的位置，有些地形横跨多个 B-Scan，本发明将其拼接后进行展示，如图(7m)，每一幅子图中的元素从上至下分别为B-Scan雷达图像、探测值、频率响应、频率幅值、坡度、几何校正后的相对回波强度、冰厚和水势。

本发明提出的检测方法对CReSIS全南极数据集测线的判断结果如图9所示。其中，深色实心圆点为使用本发明方法探测出的冰下水体，深色实心圆点的大小代表是冰下水体的可能性大小，圆点的显示阈值为前文中计算出的阈值，本发明将南极地区的阈值设置为9，探测值超过阈值的会在地图中显示，浅色实心圆点为 Livingstone等人汇总并整理的全南极地区的675个冰下湖，最新版冰下湖清单汇总了往年各个国家、各个项目的测线上的所有冰下湖，比CReSIS数据库中的测线范围更广，因此很多不在CReSIS数据库测线上的湖无法进行交叉对比验证。本发明方法与水势预测的对比结果如图10所示。

以上所述仅为本发明的优选实施方式，本发明的保护范围并不仅限于上述实施方式，凡是属于本发明原理的技术方案均属于本发明的保护范围。对于本领域的技术人员而言，在不脱离本发明的原理的前提下进行的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。