0 引言

地球重力场及其时间变化是反映地球表面和内部密度分布、质量迁移以及地球圈层之间的质交换和耦合最直接的物理量，其观测室全球和局部地球动力学过程和环境变化研究的重要依据，然而，在重力观测中包含了非常显著的潮汐变化(周江存等，2009)，其中包括重力固体潮和海洋潮汐的负荷效应.重力潮汐变化是固体地球在日、月及近地行星等天体引潮力作用下整体形变导致的地球重力场周期性变化(徐建桥等，2014).

重力固体潮主要来自于日、月天体(非常微小的部分来自于离地球较近的行星)的引潮力的作用(周江存等，2009).固体潮是唯一能够预先计算出的地球物理现象，一直以来有Kelvin引进的对比研究法则一直为国内外固体潮研究所遵循(周挚等，2008)，并且已有学者严格给出固体潮理论值得基本公式(郗钦文，1982；郗钦文和侯天敏，1986)，近年来很多学者不断改进，提出了更加精确的重力固体潮的理论值计算(董良等，2015).目前趋向于根用观测固体潮对比正常固体潮得出固体潮异常，进而研究固体潮异常的时空分布规律(周挚等，2008).固体潮主要包含地倾斜固体潮及重力固体潮两部分.近年来许多学者提出固体潮与地震具有一定的触发关系，固体潮应变力一定程度会触发地震，(周江存等，2013；孙长青等，2014)因此深入研究固体潮是及其重要的.重力固体潮信号中含有丰富的潮汐谐波分量，潮汐波按其周期可分成长周期波、日波和半日波，根据杜森公式展开可求得相应谐波分量的频率可以进行类别划分(北京大学地球物理系等，1982).目前广泛使用的固体潮分析方法Venedikov调和分析方法和ETERNA调和分析方法，近年来较为流行的方法有HHT分析法.此类方法都能分解出重力固体潮信号的一些谐波分量，但不能将这些谐波分量与重力固体潮产生机制一一对应.

本文提通过分析重力固体潮信号产生机制，提出将EDA和ICA算法结合，结合重力固体潮正交分解模型(李云飞，2015)进行对重力固体潮信号分析.ICA (Independent Component Analysis, 独立成分分析)是从多元数据中寻找其内在具有统计独立和非高斯的因子或成分的一种盲源分离方法(Hyvarinen等，2007).重力固体潮是一种复杂的多变的混合信号，因此可以使用此种方法进行处理得到解混信号，再通过Fourier变换得到观测数据的频率特征，分析其细节性的时频特征.EDA (Estimate Distribution Algorithm，分布估计算法)和传统的遗传算法不同，在对问题求解过程中并没有进行交叉和变异等操作，反而采用的是概率模型的学习和采样，通过群体的自我更新来达到全局最优(周树德和孙增圻，2007).

通过仿真实验及验证可知，本文提出的改进算法可以有效地将重力固体潮信号中含有的长周期波、日波和半日波谐波分量成功分离，同时能将分离的三种谐波对应到重力固体潮信号正交分解模型的各个方向上，从产生机制上划分潮汐谐波.同时得到的谐波时序特征清晰，可用于后序与实际观测值对比研究.

1 重力固体潮正交分解模型

在图 1中，A为地球上的某一观测点，其受到的太阳引潮力力和月球引力主要分解为地倾斜固体潮信号F_h和重力固体潮信号F_g.本文只分析重力固体潮信号F_g，对其进行正交分解，分解为一个平行于赤道平面的信号分量F₁，一个平行于地球自转轴的信号分量F₂(李云飞，2015).F₁与赤道平面平行，该信号分量仅受地球自转的影响，因此其主要含地球自转引起的独立成分谐波分量，即F₁₁日波、F₁₂半日波等.F₂和地球自转轴平行，F₂方向上的信号分量不受地球自转的影响，因此该方向主要含受月球和太阳引潮力作用产生的独立谐波成分，即长周期波.

2 改进算法对重力固体潮信号的处理 2.1 EDA算法

分布估计算法通过一个概率模型描述候选解在空间的分布，采用统计学习手段从群体宏观的角度建立一个描述解分布的概率模型，然后对概率模型随机采样产生新的种群，如此反复，实现种群的进化.根据改路模型的复杂程度以及不同的采样方法，分布估计算法主要以下两个步骤：

1)构建描述解空间的概率模型，通过对种群的评估，选择优秀的个体集合，采用统计学习等手段构造一个描述当前街机的概率模型.

2)由概率O型随机采样产生新的种群.一般采用蒙特卡洛方法，对概率模型采样得到新的种群(周树德和孙增圻，2007).

本文在此采用的EDA算法中的一类，变量无关的分布估计算法，PBIL (Population-Based Incremental Learning，基于群体增量学习)，PBIL由美国卡耐基梅隆大学的Baluja在1994年提出，用来解决二进制编码的优化问题.

该算法采用的优化变量是0-1编码，它们之间相互独立.设优化变量x=(x₁, x₂, …, x_n)，x_i∈{0,  1}，在没带选择群体中，对个体分别统计每一位上1出现的概率n_i/N，其中N是优势群体的规模，是优势群体中第i位出现1的个体数.PBIL在第k-1代概率密度函数的基础上进行线性调整，得到第k代估计的概率密度函数为

(1)

其中初始概率p₀(x_i)=0.5, 0 < α < 1为学习速率.在产生新个体是，第i位上按概率p_k(x_i)产生1，由此构成一个新的0-1编码的个体(王丽芳，2012).PBIL算法的流程图(Gao et al, 2012)，如下图 1所示：

2.2 ICA算法

ICA (Independent Component Analysis, 独立成分分析)是伴随BSS (Blind Source Separation，盲源分离)问题发展起来的信号分离算法，该算法在于寻找一个线性坐标系统，使产生的信号尽可能的彼此统计独立( 余先川和胡丹，2011 ).

ICA算法对接收到的观测值信号进行预处理、球化、再经过解混矩阵处理，得到源信号各个分量的近似估计值.进行ICA处理的目的就是找到混合信号x的一个线性变换矩阵W，使得输出尽可能的独立，即

(2)

ICA的计算模型估计通常是通过选择一个合适的目标函数，然后对其进行最小化或最大化.本文选择四阶累积量(Hyvarinen等，2007)作为目标函数，四阶累计量尽管小，但是为了经精度的需求，也是需要考虑的(周江存等，2013).四阶累积量kurt (y)定义为

(3)

式中所有的随机变量都假定是零均值的.峭度的值可正可负，非高斯性通常由峭度的绝对值来度量.峭度值越大，其信号的非高斯性就越强，独立性就越高.

2.3 基于分布估计算法对重力固体潮信号独立成分的提取

ICA算法中的关键环节是目标函数和最优化算法两个方面.本文的采用四阶累积量作为目标函数，EDA作为优化算法来实现ICA计算模型估计获得最优解混矩阵W.

其应用在重力固体潮信号分析中步骤如下：

1)获取观测点处的重力固体潮信号观测值.通过模型分析可知，重力固体潮信号由地球受到自身自转、太阳和月球的引潮力作用引起，所以分析重力固体潮信号选取3路信号X₁, X₂, X₃.

2)重力固体潮观测值预处理.对3路信号采用公式(4)去均值、公式(5)进行白化等操作.公式为

(4)

(5)

在(5)中，V是白化矩阵，D=diag[λ₁λ₂·λ_n]; E=[c₁c₂·c_n], λ_i协方差矩阵C_x=E{x(t)x^T(t)}的是第i个最大特征值(Hyvarinen等, 2007).

3)利用分布估计算法优化解混矩阵W

(1)初始化粒子的数量M及粒子的维度D，v=[v_1:D, v_2:D, …, v_m:D]^T, i=1, 2, …m.并设置迭代次数最大值N.

(2)开始计算迭代循环，令i=1，将M个二进制的粒子转变为十进制，W=(w_{M, 1}, w_{M, 2}, …, w_{M, V})

(3)利用公式(1)，，对Y使用公式(2)进行峭度计算，并按从大到小排序，得到F=(f₁, f₂, …, f_M)

(4)从F中按照峭度选择原则选择K个最优粒子，x_1:M¹, …x_i:M¹, …, x_K:M¹，利用公式(1)得到概率向量p₂(x)，利用p₂(x)进行更新粒子的分布，用于下一次的计算.

(5)进行循环计算，当i < N，返回第(2)步.否则，结束循环.得到最优分离矩阵W.

4)根据公式(2)，得到输出信号y(t).

5)得到重力固体潮信号的各谐波信号分量，将其做频谱分析.将频谱分析结果与重力固体潮正交分解模型以及相应的谐波理论值进行对比.

3 实验结果分析

根据重力固体潮正交分解模型的分析，对于观测点的选择应具有相同的经度和不同的纬度.在实验过程中，我们选择了地球的多个观测点，其分析结果是一样的.所以，本文选择2010年1月到2010年11月，同一经度上不同纬度的三个观测点1.[110°E, 20°N]，2.[110°E, 30°N]，3.[110°E, 40°N]进行重力固体潮信号分析，理论值数据由固体潮理论值封闭公式(刘序俨和李平，1986；董良等，2015)求得，其观测点分布如图 3所示：

实验中，一些参数设置如表 1所示：

观测点处重力固体潮信号波形图，如图 4所示.通过PBIL-ICA算法处理后的信号波形如图 5所示.

为了检验结果的正确性，并获得各输出信号中的谐波分量，本文对输出信号进行了频谱分析，频谱分析结果如图 6，分析过程中只研究正频部分，为了更好地观察和分析，对图 6进行局部放大.根据杜森公式展开求得的潮汐波的角频率作为理论值对比，其对比结果在表 2中给出，并划分了测量点得到的谐波分量类别.其中e代表 10的幂次方.

从表 2和图 3可知：(1)3路输出解混信号所含谐波分量主要频率与重力固体潮信号理论谐波分量频率基本一致；(2)改进的算法分离效果较好，在重力固体潮正交分解模型的基础上，算法能够将重力固体潮信号按照谐波分量系自动分离；(3)输出的3路重力固体潮信号能够与正交分解模型的相应方向对应，互不干扰；(4)日波和半日波、月波与半月波之间有明显的调制关系.

4 结果与讨论 4.1

本文提出了一种EDA和ICA的混合优化算法，在提高算法分离性能的结果基础上，用来分析重力固体潮信号.得到的解混输出信号中所含谐波分量频率信息与理论值进行对比，可知结合算法可以将重力固体潮信号分解为三个部分，相对应的为长周期波，日波和半日波.因此，本文提出的算法，在重力固体潮正交分解模型的基础上，能够将重力固体潮信号自动的按照相应的谐波系进行分离.省去了求得所有谐波系再进行划分类别的步骤，简化了对了重力固体潮信号的分析.

4.2

本文得到的谐波效果清晰，谐波包含丰富的谐波频率分量，因此可以用于与实际观测值(经过一定处理)进行对比，得到一些与地震信息有关的物理量，如强地震前会观测到震颤异常波(蒋骏等，2012).因此针对本文提出的一些观点和问题，下一步的研究方向是利用理论值和实际观测值进行对比，发现其中差异，求取地震的前兆信息或积累的异常信息.

致谢 感谢审稿老师对本文给予的指导和宝贵意见.