CN109283526A - 一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法 - Google Patents

Abstract

本发明公开一种基于全极化雷达数据的冰盖内部组构和冰流场分布特征识别方法，因冰盖内部冰晶体的二阶组构方位张量和介电常数张量存在相关关系，本发明通过逐层递推法从全极化雷达回波信号求取冰盖内部各反射层的介电常数张量及其旋转角，进而计算得到冰组构二阶方位张量特征值并识别组构类型，再根据组构、应力变形和冰流的相关关系提取冰流方向及其空间分布特征信息。该方法具有快捷、高效、无损且费用投入低等特点，且雷达探测的测量方式灵活、覆盖范围大；通过极化雷达测量可快速获得大范围冰盖内组构与冰流场空间分布信息。

Description

一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别 方法

技术领域

本发明属于雷达勘探技术领域，涉及一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法。

背景技术

冰晶组构是一种用于描述冰晶体光轴(C轴)取向几何分布统计特征的专业术语。冰晶组构是冰盖流动和演化历史的产物。冰盖流动、演化过程中冰晶体的位置和深度都会发生变化，相应地，所受应力大小和方向也随之变化，促使冰晶体发生形变，C轴聚集方向也随之改变，从而形成多种组构类型(随机型、单极大型、拉长单极大型和垂直带状结构等)。冰晶组构保存着冰盖流动、演化过程中应力、变形的历史信息，对于恢复和重建冰盖演化过程具有重要意义。另一方面，冰晶组构表现出明显的各向异性，不同组构类型受应力作用后变形程度具有明显的差异性，宏观表现为不同组构冰的黏度差异明显，对冰盖流动过程又产生反馈性影响。因此冰晶组构是冰盖流动、演化及其动力过程的关键指示物，对于理解和预测冰盖的稳定性及其演化趋势具有重要意义，成为冰盖模式研究的关键约束参数。

然而目前冰晶组构信息获取的途径只能来源于对大量深冰芯样品的实验室测试分析：通过对冰样品薄片的偏振光分析，统计冰晶体C轴方向的分布特征得到组构类型。由于极地钻探费时费力、耗资巨大，南北极目前只有为数不多的几个深冰芯钻探孔；对于广袤冰盖而言，仅从寥寥几个钻孔难以把握冰晶组构在整个冰盖内的空间分布特征和演化规律。使得组构在冰盖内的分布和演化机制成为冰盖研究中最不确定的因素。

冰下地形的起伏变化和冰晶组构类型的多样性造就了冰盖内部冰流场的复杂多样性。冰流速度和流动方向是冰盖演化过程研究的重要特征指标。地表的冰流方向和速度可以通过测量定位仪器直接观测得到，但冰盖内部冰流状况则始终是一个谜，目前没有特别有效的测试手段，只能通过冰盖模式(数值模拟)方法加以预测。但冰流信息本身就是冰盖模式研究中需要用到的关键约束参数。目前只能先作假定，将模式计算结果中冰表面流速与现场测量结果比较，验证和评价冰表部分的吻合程度；但冰盖内部应力和冰流状况分布更加复杂，特别是冰床底部附近受冰下地形影响突出。冰盖内部冰流信息的“空白”成为冰盖稳定性评价和预测过程又一个不确定的因素。

最为关键的是，由于冰下地形的起伏变化、气候变迁导致的冰晶体尺寸差异使得冰盖内部冰流场和组构分布存在明显的垂向不均一性(软/硬冰层)。而这种垂向不均一性又直接影响冰盖内部层的横向扩展和流动速度，形成更加复杂的组构类型，进而影响冰盖演化过程。而现有的技术又难以有效把握这种不均一性的分布和演化规律，成为当前冰盖研究亟待解决的关键点和技术难点。虽然目前常用的单极化冰雷达在冰盖研究中取得巨大成功，可以得到冰内等时层、冰下地形等几何参数，但对于冰盖内部组构分布、应力/应变状态以及冰流特征研究却无能为力。多极化雷达具有研究冰盖内部组构各向异性特征的潜力，但由于电磁波在多层各向异性冰层中传播理论过于复杂，严重制约着多极化雷达在冰盖研究中的应用。

发明内容

本发明针对现有技术的不足，提供一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法，具体技术方案如下：

一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法，其特征在于，该方法包括如下步骤：

S1：通过正交方式布置的全极化雷达收发天线获得四个通道雷达数据

S2：对步骤一获得的雷达数据进行归一化处理，得到四个通道的雷达反射回波记录R_HH、R_HV、R_VH、R_VV；

其中，表示H方向的入射波电场强度，表示V方向的入射波电场强度；

S3：对四个通道雷达反射回波记录R_HH、R_HV、R_VH、R_VV进行反褶积处理，得到四个通道的广义反射系数序列合在一起记为

S4：采用逐层递推法从上往下依次求取每一层的介电常数张量的特征值和偏转角，从第1层到第n-1层的推导如公式(2)所示，而界面反射系数矩阵公式和透射系数矩阵公式、逆向透射系数矩阵公式分别如式(3)-(5)所示：

其中，表示第i层与第i+1层之间的界面反射系数矩阵，表示第i层与第i+1层之间的界面透射系数矩阵，表示第i层与第i+1层之间的界面逆向透射系数矩阵，表示雷达回波记录中第i层与第i+1层之间界面的广义反射系数矩阵，表示第i层的相位因子矩阵，d_i表示第i层的厚度；因第一层地表以上为各向同性空气层，介电常数张量的特征值均为1，旋转角为0，即ε₁₁＝1、ε₁₂＝1、ε₁₃＝1和将其代入公式(1)-(3)，求解得到第二层的介电常数张量的特征值ε₂₁、ε₂₂和旋转角再通过公式(4)-(5)计算界面透射系数矩阵和逆向透射系数矩阵然后依次递推，计算出n层的介电常数张量的特征值和旋转角；

S5：根据S4得到的每层的介电常数张量的特征值ε_i1、ε_i2以及旋转角以及根据冰晶组构方位张量和介电常数张量的相关关系的公式(6)和(7)，计算每一层的二阶组构方位张量特征值a₁ ⁽²⁾、a₂ ⁽²⁾和a₃ ⁽²⁾

其中，ε_⊥和ε_//分别为垂直于C轴方向和平行于的C轴方向的介电常数，其由实验室测得；根据公式(6)、(7)计算得到的二阶组构方位张量即可用于判断冰晶组构类型；

S6：因冰晶体C轴方向往往向正压力轴方向聚集、垂直于剪切力轴、远离拉伸力轴，因此S5中得到的冰晶组构类型和C轴方向特征也反映了晶体的受力状况和形变程度。根据组构类型和C轴分布特征可以判断冰晶体的受力状况，其中拉伸力的方向即是冰流方向；逐层提取不同深度的冰层的流动方向，进一步获得冰流场的空间分布特征。

进一步地，所述的S3中的反褶积处理采用稀疏脉冲反褶积方法。

进一步地，所述的层界面反射系数矩阵公式(3)可变形为用于求解下层的介电常数张量的特征值ε_i+1,1、ε_i+1,2的公式(8)或(9)以及用于求解下层的旋转角的公式(10)或(11)

其中，公式(8)-(11)中，和为公式(3)求得的第i层与第i+1层之间的界面反射系数矩阵的四个分量。

与现有技术相比，本发明的有益效果如下：

1.快捷、高效、无损且费用投入低：本发明只需要利用全极化雷达探测数据，通过极化雷达数据处理技术直接提取组构参数和冰流信息，避免了传统深冰芯样品测试所需的高投入、长耗时；

2.雷达探测的测量方式灵活、覆盖范围大，可以得到组构三维空间分布信息。传统通过深冰芯分析只能获取钻孔所在位置一个点在深度方向的一维组构分布信息，而冰流信息即使深冰芯钻探也无法准确获取。由于深冰芯钻探点都是孤立、稀疏的，难以总结出组构在整个冰盖中的空间分布特征和演化规律。而雷达测量可以灵活布置测线(包括航空测量)，易于实现大范围、三维测量；

3.充分利用全极化雷达的技术优势，获取信息多。全极化雷达野外采集方式与传统雷达类似，施工效率高，易于实施大范围快速测量。全极化雷达测量不仅可以得到常规雷达的冰内等时层、冰厚和冰下地形等信息，还可以进一步得到组构、冰流场在冰盖内的空间分布等重要信息。这些关键信息的有机融合，有助于进一步理解冰盖的演化机制和动力过程。更重要的是，雷达测线可以将孤立的深冰芯钻探点有机连接起来，充分利用已有钻孔资料实现对比验证、校核和相互补充。

附图说明

图1为本发明基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法的流程示意图；

图2为本发明用到的冰盖内部应力分布与组构演化示意图；

图3为本发明用到的冰盖内部组构演化与冰流方向关系示意图；

图4为本发明用到的全极化雷达正交天线及其天线组合方式；

图5为本发明用到的全极化雷达分量剖面对比图。

具体实施方式

下面根据附图和优选实例详细描述本发明，本发明的目的和效果将变得更加明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，为基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法的流程图。识别技术1的关键是利用微观组构方位张量和介电常数张量的相关关系，建立宏观全极化雷达得到的反射系数和微观组构物性的关系表达式。这样通过全极化雷达探测数据直接计算并识别出组构类型，识别技术2是基于微观冰晶组构的形成原因是应力作用和变形的结果。利用应力、组构和冰流的相关关系，建立冰流方向的识别方法，继而得到冰流的空间分布特征。

1.冰晶组构的两种物性表达方式及其相关关系

宏观巨厚的冰盖微观上是由数量巨大的冰晶体集聚而成。最初的降雪中冰晶体C轴往往随机分布。由于应力作用弱，在深度小于700m的冰盖浅部区域，组构类型往往仍是随机型。随着重力和侧向剪切力的逐渐增强，使得冰晶体C轴逐渐趋于优势方向(与应力方向相关)，出现多种组构类型：单极大、强单极大、拉长单极大以及垂直带状结构等，如图2所示。不同的组构类型往往出现在冰盖特定的区域和深度：如单极大往往出现在冰穹区域，拉长单极大是由于垂向正压力和微小侧向剪切力共同作用的结果。垂直带状结构多出现在冰盖中下游，是强烈侧向剪切力作用的结果。冰晶组构是冰盖内部冰晶体受长期应力作用的产物，反映了冰盖内部动力过程和变形历史。单极大和垂直带状结构是组构演化过程中侧向剪切作用的两种极端情况，而拉长单极大则介于两者之间，其拉长程度则反映所受侧向剪切力作用程度和作用历史，如图2-3所示。实验结果表明，冰晶体C轴方向和垂直于C轴方向抗剪切力的作用具有数十倍量级的差异，垂直于C轴方向更容易滑动变形。因此冰晶组构、应力/形变和冰流场分布之间相互密切关联。

就光学性质而言，冰晶组构可以用晶体C轴方向的统计参数——二阶方位张量来表征：9个参数中只有6个参数独立的二阶对称矩阵。

该二阶对称矩阵的特征向量指示主轴方向，特征值a₁ ⁽²⁾,a₂ ⁽²⁾和a₃ ⁽²⁾就是主轴坐标系中三个主轴分量，且满足以下关系

就电磁性质而言，冰晶组构又可以用二阶介电常数张量来表征

类似地，在主轴坐标系中，介电常数张量只保留主对角线方向的三个介电常数张量特征值ε₁,ε₂和ε₃

表1组构类型与组构方位张量及介电常数张量的关系

对于组构识别而言，方位张量和介电常数张量的三个特征值都可以用于判别组构类型。两者具有相似的参数表达形式，且存在以下相关关系

其中为折射率张量，为组构方位张量，分别为冰晶体C轴方向以及垂直于C轴方向的折射率(对应温度-21℃时的测量值)，Δn＝n_//-n_⊥为两者之差。折射率和介电常数的关系是

在主轴方向，公式(5)简化为组构方位张量特征值a₁ ⁽²⁾,a₂ ⁽²⁾和a₃ ⁽²⁾和介电常数张量特征值ε₁,ε₂和ε₃的相关关系表达式

其中，ε_⊥和ε_//分别为垂直于C轴方向和平行于的C轴方向的介电常数，其由实验室测得。ε_ij表示第i层介电常数张量j方向特征值。

因此只要得到介电常数张量ε_i1和ε_i2，即可通过公式(7)换算出组构方位张量特征值a₁ ⁽²⁾,a₂ ⁽²⁾。再根据(8)式求取a₃ ⁽²⁾。对于组构类型识别而言，换算成组构方位因子或者直接用介电常数(或折射率)值都可以用于判断组构类型，冰晶组构光学和电磁物性的相关关系，即公式(5)-(8)构成了极化雷达组构识别技术的理论基础。

2.全极化雷达组构和冰流特征参数提取和识别技术

全极化雷达收发天线采用正交方式布置，如图4所示，通过收发天线组合得到四个通道的雷达数据冰晶体具有明显的各向异性特征：不同极化方向的雷达回波呈现明显的幅度差异，利用这种幅度差异就可以逆向反演出冰盖内部的介电常数张量，进而得到组构类型及其空间分布。图5a和b显示VV和VH正交方向雷达回波由于冰组构各向异性引起的不同极化方向雷达回波的能量差异，图5c为正交极化方向的平均值剖面，模糊了正交极化方向的能量差异。图5d为Stokes向量S1分量(与组构相关的参数)剖面图，显示组构在冰盖内的分布及起伏变化，特别是由于冰期间冰期旋回导致的软硬冰交替层的存在。但Stokes向量只能定性说明组构的空间分布特征，目前还不能给出用于直接判别组构类型的量化相关关系。

2.1组构类型识别技术

全极化雷达的雷达回波可以用反射系数矩阵和入射波乘积来表示。

其中，表示H方向的入射波电场强度，表示V方向的入射波电场强度；公式(9)中反射回波的电场强度除了与入射波电场强度有关外，还与反射系数相关，而反射系数是由地下介质的物性以及天线与介质的几何位置关系决定。反射波形成的原因是上下层波阻抗差异(冰盖内主要是介电常数差异)。因此，通过反射系数序列可以逆向求取各层的介电常数参数。

极化雷达接收到的反射系数序列一般表达式为

其中表示第i和i+1层之间界面的回波反射系数矩阵，表示t_i时刻的冲击脉冲响应矩阵，t_i表示该界面回波的延迟时间。广义反射系数可以用多层介质之间多次反射、透射的几何光学级数来表示，由于冰盖内等时冰层介电常数差异较小，可以忽略多次波，只考虑一次反射波，则各层的广义反射系数频域表达式具有如下形式：

层界面反射系数、透射系数以及逆向透射系数计算公式分别如式(13)、(14)-(15)所示：

其中，表示第i层与第i+1层之间的界面反射系数矩阵，表示第i层与第i+1层之间的界面透射系数矩阵，表示第i层与第i+1层之间的界面逆向透射系数矩阵，表示雷达回波记录中第i层与第i+1层之间界面的广义反射系数矩阵，表示第i层的相位因子矩阵，d_i表示第i层的厚度；这里采用逐层递推法，从上往下递推，首先从第一层和第二层之间的反射界面开始，因第一层地表以上为各向同性空气层，介电常数张量的特征值为1，旋转角为0，即ε₁₁＝1、ε₁₂＝1、ε₁₃＝1和将其代入层界面反射系数公式，联立公式(10)、(12)-(13)，可得到第二层介电常数张量的特征值ε₂₁、ε₂₂和旋转角再通过公式(14)、(15)求得界面透射系数矩阵和逆向透射系数矩阵将求得的以及第二层参数介电常数张量ε₂₁、ε₂₂和旋转角和已知的带入公式(12)求得再带入层界面反射系数公式(13)求解方程组，得到第三层参数介电常数张量ε₃₁、ε₃₂和旋转角再根据公式(13)-(15)顺序求得R₃₂和再和已知的带入公式(12)求得依次递推，计算出n层的介电常数张量的特征值和旋转角；

根据上述得到的每层的介电常数张量的特征值ε_i1、ε_i2以及旋转角根据冰晶组构方位张量和介电常数张量的相关关系的公式(7)和(8)计算每一层的二阶组构方位张量特征值a₁ ⁽²⁾、a₂ ⁽²⁾和a₃ ⁽²⁾，并根据表1识别出组构类型。各层的旋转角用于构建和指示组构和冰流的方向和空间分布特征。

3.冰流方向与物性参数的相关关系

冰的流动是应力作用下微观冰晶体产生形变的宏观表象。应力作用于冰晶体，一方面促使冰晶体产生形变，形成流动现象；另一方面促使冰晶体组构演化，形成新的组构类型。因此冰的流动过程和流动方向与组构演化和晶体变形密切相关，都是应力作用于冰晶体的结果。二阶组构方位张量(光学物性参数)和二阶介电常数张量(电磁物性参数)都将随着冰盖内应力场的变化而变化，也隐含着可以用这两种参数判断应力作用方向和冰流方向。

冰盖内主要存在两种力的共同作用：由于重力产生的垂向正压力和由于坡度而产生的侧向剪切力(或侧向拉伸力)。前人实验研究表明：冰晶体C轴方向往往向正压力轴方向聚集、垂直于剪切力轴、远离拉伸力轴。因此组构的形态也反映了晶体的受力程度和形变程度，受力方向则指示着冰的流动方向。

图2和图3显示冰盖内部冰晶体在正应力和侧向剪切力作用下组构的变化趋势。在冰穹区域，侧向剪切力很小，重力作用的正应力使得晶体C轴逐渐聚集，趋于优势方向，形成一个极端：与表1中单极大情形吻合。此时，介电常数张量和组构张量的XY轴与水平面重合、C轴逐渐趋于垂向聚集(Z轴方向)，水平面内介电常数值相等。随着冰往下游流动，侧向剪切力逐步增强(主要为侧向拉伸)，剪切力拉长了C轴面(由圆形拉长为椭圆)，椭圆短轴指示侧向拉伸力和冰流方向(X轴方向)。当侧向拉伸变形足够大(或组构拉伸到一定程度)，由拉长单极大走向另一个极端，形成垂直带状结构。此时，C轴在垂直于冰流方向的YZ面内任意分布(图3中X方向为侧向拉伸力方向和冰流方向)。因此可以根据组构的形态，即组构Schmidt图中椭圆的短轴方向以及带状结构中垂直于条带的方向来判断冰流方向。

4.冰流方向识别技术

逐层递推法得到各层组构方位张量特征值a₁ ⁽²⁾,a₂ ⁽²⁾和a₃ ⁽²⁾以及旋转角不仅给出了组构类型，也给出了组构的偏转方向。因组构的拉伸方向与应力变形以及冰流方向密切相关，也指示了侧向拉伸力的方向和冰流方向，根据计算得到的a₁ ⁽²⁾,a₂ ⁽²⁾大小对比以及旋转角可进一步给出每一层的冰流方向。全极化雷达往往是剖面或三维测量，将所有测点深度方向的冰流方向合在一起就可以构建和分析冰盖内部冰流方向的空间分布特征。

Claims (3)

1.一种基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法，其特征在于，该方法包括如下步骤：

S1：通过正交方式布置的全极化雷达收发天线获得四个通道雷达数据

S2：对步骤一获得的雷达数据进行归一化处理，得到四个通道的雷达反射回波记录R_HH、R_HV、R_VH、R_VV；

其中，表示H方向的入射波电场强度，表示V方向的入射波电场强度；

S3：对四个通道雷达反射回波记录R_HH、R_HV、R_VH、R_VV进行反褶积处理，得到四个通道的广义反射系数序列合在一起记为

S4：采用逐层递推法从上往下依次求取每一层的介电常数张量的特征值和偏转角，从第1层到第n-1层的推导如公式(2)所示，而界面反射系数矩阵公式和透射系数矩阵公式、逆向透射系数矩阵公式分别如式(3)-(5)所示：

其中，表示第i层与第i+1层之间的界面反射系数矩阵，表示第i层与第i+1层之间的界面透射系数矩阵，表示第i层与第i+1层之间的界面逆向透射系数矩阵，表示雷达回波记录中第i层与第i+1层之间界面的广义反射系数矩阵，表示第i层的相位因子矩阵，d_i表示第i层的厚度；因第一层地表以上为各向同性空气层，介电常数张量的特征值均为1，旋转角为0，即ε₁₁＝1、ε₁₂＝1、ε₁₃＝1和将其代入公式(1)-(3)，求解得到第二层的介电常数张量的特征值ε₂₁、ε₂₂和旋转角再通过公式(4)-(5)计算界面透射系数矩阵和逆向透射系数矩阵然后依次递推，计算出n层的介电常数张量的特征值和旋转角；

S5：根据S4得到的每层的介电常数张量的特征值ε_i1、ε_i2以及旋转角以及根据冰晶组构方位张量和介电常数张量的相关关系的公式(6)和(7)，计算每一层的二阶组构方位张量特征值a₁ ⁽²⁾、a₂ ⁽²⁾和a₃ ⁽²⁾

其中，ε_⊥和ε_//分别为垂直于C轴方向和平行于的C轴方向的介电常数，其由实验室测得；根据公式(6)、(7)计算得到的二阶组构方位张量即可用于判断冰晶组构类型；

S6：因冰晶体C轴方向往往向正压力轴方向聚集、垂直于剪切力轴、远离拉伸力轴，因此S5中得到的冰晶组构类型和C轴方向特征也反映了晶体的受力状况和形变程度。根据组构类型和C轴分布特征可以判断冰晶体的受力状况，其中拉伸力的方向即是冰流方向；逐层提取不同深度的冰层的流动方向，进一步获得冰流场的空间分布特征。

2.根据权利要求1所述的基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法，其特征在于，所述的S3中的反褶积处理采用稀疏脉冲反褶积方法。

3.根据权利要求1所述的基于全极化雷达的冰盖内部组构和冰流场分布特征识别方法，其特征在于，所述的层界面反射系数矩阵公式(3)可变形为用于求解下层的介电常数张量的特征值ε_i+1,1、ε_i+1,2的公式(8)或(9)以及用于求解下层的旋转角的公式(10)或(11)

其中，公式(8)-(11)中，和为公式(3)求得的第i层与第i+1层之间的界面反射系数矩阵的四个分量。