CN103577607A - 一种基于地磁异常数据形态特征的边界补偿方法 - Google Patents

Abstract

本发明属于地磁导航领域，特别涉及一种基于地磁异常数据形态特征的边界补偿方法。基于地磁异常数据形态特征的边界补偿方法，包括：对实测地磁异常数据进行规则网格化处理；平移地磁异常网格数据集，使均值为零；对平移后的零均值网格数据集进行二维经验模态分解；对二维经验模态分解结果进行形态特征提取；分割数据块单元，实时消除相邻单元的特征相似性冗余；消除所有数据块单元间的特征相似性冗余；构建多层级地磁异常数据边界补偿数据库；对待分析地磁异常数据进行边界补偿。本发明能够有效的改善区域地磁异常数据分析中的边界效应问题，相比与其他方法有更好的适用性和使用上的便捷性。

Description

一种基于地磁异常数据形态特征的边界补偿方法

技术领域

本发明属于地磁导航领域，特别涉及一种基于地磁异常数据形态特征的边界补偿方法。

背景技术

地磁导航是一种新兴的无源、自主导航方式，相比传统导航方法具有全天候、全地域、隐蔽性好、无积累误差等优点。地磁异常场是地磁导航中的主要参考量，由地壳表面分布和局部地理特征如矿产、岩石、人造磁场等产生，空间结构复杂，磁异常强度随空间位置变化明显，且在时间上变化非常稳定，由许多大小尺度不等、正负相间的磁场分布区组成，精确的地磁异常数据分析是实现高精度地磁导航的关键。由于进行精细的磁异常测量难度较大、测量过程非常复杂，因此目前所掌握的地磁异常场数据多呈区域性分布，从而导致利用实测数据插值重构导航用基准图或进行多尺度分解操作时面临着边界效应的重大影响。边界效应是指对区域性地磁异常数据进行插值或分解等操作时，由于边界区域缺少参考信息，从而造成边界区域分析结果误差过大甚至发散的现象，而这种误差甚至会深入到数据内部，从而大大降低分析结果的利用价值。因此，本专利发明一种边界补偿方法，旨在增加地磁异常数据边界处的参考信息，减小边界区域的分析误差，提高地磁导航的可靠性。

解决边界效应问题的本质方法就是在所研究数据周围添加适量可靠的参考点，很多学者对该问题进行了研究。乔玉坤分别使用了支持向量机和BP神经网络方法预测导航区域边缘位置的基准数据，从而减弱边界效应的影响；谭斌利用径向基函数进行边界插值约束，减小了所构建区域地磁场模型的预测偏差；杨云涛采用边界插值约束的Tayler多项式拟合法进行区域地磁场建模，有效改善了建模过程中存在的边界效应问题；在空间数据多尺度分解方面，经验模态分解法是非常优秀且常用的方法，在二维空间数据分解过程中，受边界效应的影响，所构建的包络面产生畸变，从而导致多尺度分解准确性的降低，镜像法、神经网络预测、多项式外延法、灰度预测等方法被用来解决该问题。但是这些方法都存在很强的局限性，尤其是对于随机分布特性较强的地磁异常场，无论是支持向量机、神经网络、灰度预测等方法的模型预测能力，还是镜像法、多项式外延等补偿方法的准确性都将大大降低。

本专利发明了一种全新的基于地磁异常数据形态相似性的边界补偿方法。地磁异常场作为一种广泛分布的空间变量，分布不规则形强，同时又有着很强的自相似性，Spector等进行了航磁数据谱的统计分析，表明地磁异常具有明显的分形特征。本专利首先在对所研究地磁异常数据进行零均值平移的基础上，进行了形态特征的提取，避免了不同数据集之间由于磁异常强度不同从而导致的特征提取结果的差异；将数据集分割成众多数据块单元，在消除单元之间相似性冗余的基础上，建立了具有普遍适用性的多层级地磁异常数据补偿单元数据库。对于任何待分析的地磁异常数据，均可以从该数据库中提取补偿单元进行边界补偿。同时该数据库具有良好的可扩充性，随着实测数据量的增多，可不断对其进行丰富。本方法极大的提高了边界补偿的便捷性；同时相比于已有方法，能够更加有效的改善地磁异常数据分析时的边界效应问题，增强分析结果的可靠性。

发明内容

本发明的目的在于提供一种进一步提高边界补偿的准确性和可靠性的基于地磁异常数据形态特征的边界补偿方法。

本发明的目的是这样实现的：

一种基于地磁异常数据形态特征的边界补偿方法，

（1）对实测地磁异常数据进行规则网格化处理：

选用克里金法作为网格化插值方法，网格化后所得规则网格数据集D₀(x,y)在测线方向上间隔为0.001°，在测线垂直方向上的间隔为0.002°；

（2）平移地磁异常网格数据集，使均值为零：

计算规则网格数据集D₀(x,y)的均值M₀，将数据集在xy平面垂直方向上平移|M₀|距离，使得平移后所得零均值网格数据集D(x,y)的均值M为零；

（3）对平移后的零均值网格数据集进行二维经验模态分解：

对平移后的零均值网格数据集进行二维经验模态分解，得到本征模态函数分量和一个残余项；

（4）对二维经验模态分解结果进行形态特征提取：

利用直接微分法对分解所得前两个本征模态函数分量进行双向瞬时频率和幅值特征的提取，利用其他分量进行趋势项幅值特征的计算；

（5）分割数据块单元，实时消除相邻单元的特征相似性冗余：

去除D(x,y)的边界区域，确定构建边界补偿数据库数据块单元的大小，在剩余区域中逐个进行单元的分割；实时将分割出的每一个单元与周围所有已分割出的单元块进行形态特征对比；

（6）消除所有数据块单元间的特征相似性冗余：

对所有样本数据集完成特征提取和单元块分割之后，将所有单元构成一个序列，逐个进行单元之间的形态特征差异性对比分析，逐步消除单元间的所有冗余信息；

（7）构建多层级地磁异常数据边界补偿数据库；

对最终所得的所有数据块单元计算地磁异常强度分布范围、平均振幅、极值点数三个整体特征，分别对这三个特征划分区间范围，构建多层级地磁异常数据边界补偿数据库；

（8）对待分析地磁异常数据进行边界补偿；

按照步骤（1）至步骤（4）提取待分析数据的形态特征，去除零均值网格数据集的边界部分，在剩余数据中划分边界补偿区域，对每个补偿区域，在地磁异常边界补偿数据库中提取数据块单元与其进行特征匹配，直至完成边界补偿。

对平移后的零均值网格数据集进行二维经验模态分解包括，

（3.1）初始化：r_i-1(x,y)＝D(x,y)，i＝1；

（3.2）抽取imf_i(x,y)，即第i个本征模态函数分量；

（3.2.1）初始化：h₀(x,y)＝r_j-1(x,y)，j＝1；

（3.2.2）计算h_j-1(x,y)的平均包络：

a)确定h_j-1(x,y)的极大值集合和极小值集合，并分别用径向基函数插值法构建极大值包络面E_upper(x,y)和极小值包络面E_lower(x,y)；

b)计算平均包络m_j-1(x,y)＝(E_upper(x,y)+E_lower(x,y))/2；

（3.2.3）h_j(x,y)＝h_j-1(x,y)-m_j-1(x,y)，令j＝j+1；

（3.2.4）当满足终止条件 $\mathrm{SD}:=\underset{x=0}{\overset{X}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{Y}{\mathrm{\&Sigma;}}}\left[\frac{{|h_{j-1}(x,y)-h_j(x,y)|}^2}{h_{j-1}^2(x,y)}\right]<0.05$ 或SD值增大，令imf_i(x,y)＝h_j(x,y)，否则令j＝j+1，重新执行步骤2）；

（3.3）r_i(x,y)＝r_i-1(x,y)-imf_i(x,y)；

（3.4）判断r_i(x,y)在θ或θ的垂直方向上是否存在单调一维采样：

v_1,c(x)＝u(x,[tanθ]x+c)，   0≤θ＜π/2

$v_{2,c}\left(x\right)=\left\{\begin{array}{cc}u(\&CenterDot;,x)& \mathrm{\&theta;}=0\\ u(x,[\tan (\mathrm{\&theta;}+\mathrm{\&pi;}/2)]x+c),& 0<\mathrm{\&theta;}<\mathrm{\&pi;}/2\end{array}\right.$

其中

v_1,c(x)和v_2,c(x)为r_i(x,y)在两个方向上的一维采样，如果存在，则分解过程终止，否则令i＝i+1，重新执行步骤2，直到得到n个本征模态函数分量IMF₁,IMF₂,…,IMF_n和一个残余项RES。

对二维经验模态分解结果进行形态特征提取包括：

获取在IMF₁和IMF₂中所有数据点在横、纵双向上的瞬时幅值和瞬时频率特征

s(t)为IMF₁和IMF₂中的网格线数据集，确定网格线数据集的极大值点集合和极小值点集合，并分别构建上包络线L_upper和下包络线L_lower，A(t)表示此网格线方向上的瞬时振幅，t表示x或y方向；

$A\left(t\right)=\left\{\begin{array}{cc}{L}_{\mathrm{upper}}\left(t\right)& s\left(t\right)>0\\ L_{\mathrm{lower}}\left(t\right)& s\left(t\right)<0\end{array}\right.,S\left(t\right)=s\left(t\right)/A\left(t\right)=\cos \mathrm{\&phi;}\left(t\right),$

S'(t)＝-φ'(t)sinφ(t)＝-ωsinφ(t)，

当sinφ(t)≠0，由ω的非负性可得，

$\mathrm{\&omega;}=-S^{\&prime;}\left(t\right)/\sin \mathrm{\&phi;}\left(t\right)=\left|S^{\&prime;}\left(t\right)\right|/\sqrt{1-S^2\left(t\right)}$

当sinφ(t)＝0，即|S(t)|＝1时，

S''(t)＝-φ''(t)sinφ(t)-[φ'(t)]²cosφ(t)＝-ω²S(t)

$\mathrm{\&omega;}=\sqrt{|S^{\&prime;\&prime;}\left(t\right)/S\left(t\right)|}=\sqrt{\left|S^{\&prime;\&prime;}\left(t\right)\right|},$ ω为瞬时频率，

将其余本征模态函数分量IMF₃,IMF₄,…,IMF_n和残余项RES相加作为趋势项，

$\mathrm{Trend}=\underset{i=3}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{IMF}}_i+\mathrm{RES}$

将Trend值作为对应数据点的趋势项幅值特征。

分割数据块单元，实时消除相邻单元的特征相似性冗余包括：

在零均值网格数据集中划定边界区域，分割成数据块单元，并保证IMF₁和IMF₂中每条网格线数据集的两端在边界区域内至少有一个极小值点和极大指点，将边界区域去除，保留中心区域作为有效区域，确定数据块单元在两个方向上均包含10个数据点，采用弓字形分割方式进行数据块单元的分割，且相邻数据块单元的对应边界之间为1倍的网格间距，计算其与相邻已分割出数据块之间形态特征的差异性：

$\mathrm{diff}=\frac{1}{3}\&CenterDot;(\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(F_{0\mathrm{ij}}^k-F_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(F_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(A_{0\mathrm{ij}}^k-A_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(A_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{0\mathrm{ij}}^k-T_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_{0\mathrm{ij}}^k\right)}^2}})$

M和N分别表示数据块单元在横向和纵向上数据点数，

和

分别表示待比较单元块与相邻单元块中对应数据点在IMF1和IMF2中瞬时频率特征，

和表示幅值特征，

和

表示趋势项幅值特征，若diff＜0.05，则将新分割的数据块忽略。

消除所有数据块单元间的特征相似性冗余包括：

（5.1）将数据块单元构建成一个序列U₁,U₂,…,U_n；

（5.2）将U1作为基准节点，依次与后面所有节点进行特征差异性比较，如果差异性小于0.05，则将被比较的节点舍弃，否则保留，比较完成后得到一新序列U₁,U₂,…,U_m，m≤n；

（5.3）依次在新序列中将U₂,U₃,…作为基准节点，与后面所有节点进行特征差异性比较，当所有节点都作为基准节点被比较过之后，得到最终序列U₁,U₂,…,U_k。

构建多层级地磁异常数据边界补偿数据库包括：

（6.1）若对测线垂直方向的边界进行补偿，则将数据块单元划分为上下两个5×10的匹配区域，即区域1、区域3；

（6.2）若对测线方向的边界进行补偿，则将数据块单元划分为左右两个5×10的匹配区域，即区域2、区域4；

（6.3）若对数据集四个拐角处进行补偿，则将数据块单元划分为4个5×5的匹配区域，即区域5、区域6、区域7、区域8，对每个数据块单元进行强度分布范围AR、平均振幅AM，、极值点数EN三个总体特征计算，

AR＝max(f(x_i,y_j))-min(f(x_i,y_j))，1≤i≤m，1≤j≤n；

$\mathrm{AM}=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|f(x_i,y_i)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(f(x_i,y_i)\right)/(m\&CenterDot;n)|$

按照强度分布范围、平均振幅、极值点数的顺序建立多层级数据库：

1）强度分布范围按100nT作为间隔进行区间划分，即[0 100]，[100 200]，……，分割得到的所有数据块单元的强度分布范围为m个区间范围；

2）在步骤1）得到的每个区间范围的基础上，按平均振幅进行区间划分，以50nT作为间隔进行区间划分，即[0 50]，[50 100]，……，划分n个区间范围；

3）在步骤2）得到的每个区间范围的基础上，按每个数据块单元中所包含的极值点数进行区间划分，以2作为间隔进行区间划分，即[0 1]，[2 3]，……，划分为k个区间范围，最终分割所得所有单元划分为m·n·k个数据块单元集合。

本发明的有益效果在于：本发明能够有效的改善区域地磁异常数据分析中的边界效应问题，相比与其他方法有更好的适用性和使用上的便捷性。在数据库构建过程中，充分利用掌握的所有样本数据集，通过提取不依赖于地磁异常强度值的形态特征，极大地丰富了边界补偿资源，为所有待分析地磁异常数据提供了统一的边界补偿来源。该数据库具有良好的扩充性，随着实测地磁异常数据的增加，可不断丰富此数据库所包含的数据资源，进一步提高边界补偿的准确性和可靠性。

附图说明

图1是本发明的方法流程图；

图2是地磁异常数据单元块分割方式示意图；

图3是单元块间的特征相似性冗余消除过程示意图；

图4是数据块单元特征匹配区域划分示意图，(a)是对边界处5×10矩形区域补偿时的匹配区域分割方式；(b)对拐角处5×5正方形区域补偿时的匹配区域分割方式；

图5是地磁异常数据边界补偿数据库的多层级结构图；

图6是地磁异常数据边界补偿示意图。

具体实施方式

一种基于地磁异常数据形态特征的边界补偿方法，包括以下几个步骤:

步骤一：对实测地磁异常数据进行规则网格化处理；

确定横向、纵向上的网格间距，采用高精度插值方法对实测数据进行规则网格化处理；

步骤二：平移地磁异常网格数据集，使均值为零；

在垂直于数据平面方向上，平移数据集，使得平移后所得数据集D(x,y)的均值为零；

步骤三：对平移后的数据集进行二维经验模态分解；

对平移后的数据集进行二维经验模态分解（BEMD），得到若干本征模态函数分量和一个残余项；

步骤四：采用直接微分法对分解结果进行形态特征提取；

利用直接微分法对分解所得前两个本征模态函数分量进行双向瞬时频率和幅值特征的提取，利用其他分量进行趋势项幅值特征的计算；

步骤五：分割数据块单元，实时消除相邻单元的特征相似性冗余；

去除D(x,y)的边界区域，确定构建边界补偿数据库数据块单元的大小，在剩余区域中逐个进行单元的分割；实时将分割出的每一个单元与周围所有已分割出的单元块进行形态特征对比，当差异性小于一定值时，则忽略该单元；

步骤六：消除所有数据块单元间的特征相似性冗余；

对所有样本数据集完成特征提取和单元块分割之后，将所有单元构成一个序列，逐个进行单元之间的形态特征差异性对比分析，逐步消除单元间的所有冗余信息；

步骤七：构建多层级地磁异常数据边界补偿数据库；

对最终所得的所有数据块单元计算地磁异常强度分布范围、平均振幅、极值点数三个整体特征，分别对这三个特征划分区间范围，构建多层级地磁异常数据边界补偿数据库；

步骤八：对待分析地磁异常数据进行边界补偿；

按照步骤一至步骤四提取待分析数据的形态特征，去除零均值网格数据集的边界部分，在剩余数据中划分边界补偿区域，对每个补偿区域，在地磁异常边界补偿数据库中提取数据块单元与其进行特征匹配，直至完成边界补偿。

步骤二通过将地磁异常数据集进行零均值平移，从而避免了具有相同形态的数据集由于磁异常强度不同导致的形态特征提取的差异。

步骤三利用二维经验模态分解方法对地磁异常数据集进行多尺度分解，从而得到若干本征模态函数分量和一个残余项。

步骤四根据个本征模态函数分量中所含极点数量，确定将分解所得的前两个分量作为特征提取分量，利用直接微分法对分量中每个点在两个方向上进行瞬时幅值和瞬时频率的提取，以第一个分量IMF₁为例，设s(t)为其中一网格线数据集，确定该数据集的极大值点集合和极小值点集合，并分别利用样条函数法构建上、下包络线L_upper和下L_lower，令A(t)（t表示x或y方向）表示此网格线方向上的瞬时振幅，如式(3)所示；

$A\left(t\right)=\left\{\begin{array}{cc}{L}_{\mathrm{upper}}\left(t\right)& s\left(t\right)>0\\ L_{\mathrm{lower}}\left(t\right)& s\left(t\right)<0\end{array}\right.---\left(3\right)$

令S(t)＝f(t)/A(t)＝cosφ(t)，方程两边对t求导，有

S'(t)＝-φ'(t)sinφ(t)＝-ωsinφ(t)   (4)

当sinφ(t)≠0，由ω的非负性可得，

$\mathrm{\&omega;}=-S^{\&prime;}\left(t\right)/\sin \mathrm{\&phi;}\left(t\right)=\left|S^{\&prime;}\left(t\right)\right|/\sqrt{1-S^2\left(t\right)}---\left(5\right)$

当sinφ(t)＝0，即|S(t)|＝1时，对(4)两边再次求导得，

S''(t)＝-φ''(t)sinφ(t)-[φ'(t)]²cosφ(t)＝-ω²S(t)   (6)

从而可得，

$\mathrm{\&omega;}=\sqrt{|S^{\&prime;\&prime;}\left(t\right)/S\left(t\right)|}=\sqrt{\left|S^{\&prime;\&prime;}\left(t\right)\right|}---\left(7\right)$

根据式(3)可得分量中所有数据点在其所在的两条网格线方向上的瞬时振幅特征，根据(5)或(7)可得瞬时频率特征。利用此方法，可完成对IMF₁和IMF₂中所有网格线数据集的特征提取，即得到IMF₁和IMF₂中所有数据点在横、纵双向上的瞬时幅值和瞬时频率特征。

将其余几个本征模态函数分量IMF₃,IMF₄,…,IMF_n和残余项RES直接相加作为趋势项，如式(8)，

$\mathrm{Trend}=\underset{i=3}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{IMF}}_i+\mathrm{RES}---\left(8\right)$

将Trend值作为对应数据点的趋势项幅值特征，从而所研究数据集中每个数据点具有9个特征。

步骤五对零均值网格数据集中划定边界区域，并保证IMF₁和IMF₂中每条网格线数据集的两端在该边界区域内至少有一个极小值点和极大指点，将该边界区域去除，保留中心区域作为有效区域。确定数据块单元在两个方向上均包含10个数据点，采用“弓”字形分割方法进行数据块的分割，且相邻数据块的对应边界之间为1倍的网格间距。

实时消除相邻数据块单元之间的特征相似性冗余，当分割出一个数据块后，根据式(9)计算其与相邻已分割出数据块之间形态特征的差异性，

$\mathrm{diff}=\frac{1}{3}\&CenterDot;(\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(F_{0\mathrm{ij}}^k-F_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(F_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(A_{0\mathrm{ij}}^k-A_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(A_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{0\mathrm{ij}}^k-T_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_{0\mathrm{ij}}^k\right)}^2}})---\left(9\right)$

式中，M和N分别表示数据块单元在横向和纵向上数据点数，

和

分别表示待比较单元块与某个相邻单元块中对应数据点在IMF1和IMF2中瞬时频率特征，相应的，和表示幅值特征，

和

表示趋势项幅值特征，若diff＜0.05，则将新分割的数据块忽略。

步骤六消除所有数据块单元之间的特征相似性冗余，包括以下步骤，

设将所有地磁异常数据集进行分割并消除相邻数据块之间冗余后共得n个数据块单元，

1.将这些单元构建成一个序列U₁,U₂,…,U_n；

2.将U₁作为基准节点，根据式(9)，依次计算其与后面所有节点进行特征差异性比较，如果差异性小于0.05，则将被比较的节点舍弃，否则保留，比较完成后得到一新序列U₁,U₂,…,U_m，m≤n；

3.根据2的比较方法，依次在新序列中将U₂,U₃,…作为基准节点，不断进行筛选，当所有节点都作为基准节点被比较过之后，得到最终序列U₁,U₂,…,U_k。

步骤七构建了多层级地磁异常数据边界补偿数据库。

确定边界补偿时采用5行网格数据进行特征匹配，根据不同边界位置的补偿需求，按三种方式对数据块单元中进行匹配区域划分：（1）若对测线垂直方向的边界进行补偿，则将数据块单元划分为上下两个5×10的匹配区域，即区域1、3；（2）若对测线方向的边界进行补偿，则将数据块单元划分为左右两个5×10的匹配区域，即区域2、4；（3）若对数据集四个拐角处进行补偿，则将数据块单元划分为4个5×5的匹配区域，即区域5、6、7、8。为提高数据库的层次性和边界补偿时对数据库的搜索效率，对每个数据块单元进行强度分布范围（AR）、平均振幅（AM，在零均值前提下计算）、极值点数（EN）三个总体特征计算，AR和AM计算方式分别如式(10)、(11)所示，

AR＝max(f(x_i,y_j))-min(f(x_i,y_j))，1≤i≤m，1≤j≤n   (10)

$\mathrm{AM}=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|f(x_i,y_i)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(f(x_i,y_i)\right)/(m\&CenterDot;n)|---\left(11\right)$

按照强度分布范围、平均振幅、极值点数的顺序建立多层级数据库，规则如下：

1.强度分布范围按100nT作为间隔进行区间划分，即[0 100]，[100 200]，……，设分割得到的所有单元的强度分布范围可划分为m个区间范围；

2.在1所得每个区间范围的基础上，进一步按平均振幅进行区间划分，以50nT作为间隔进行区间划分，即[0 50]，[50 100]，……，设共划分n个区间范围；

3.在2所得每个区间范围的基础上，进一步按每个单元中所包含的极值点数进行区间划分，以2作为间隔进行区间划分，即[0 1]，[2 3]，……，设共划分k个区间范围。

综上，可将分割所得所有单元划分为m·n·k个单元集合。

步骤八对待分析地磁异常数据的边界区域作为特征匹配区域进行补偿，具体方法如下：

1.按照步骤一至步骤四方法对待分析数据进行规则网格化、零均值平移、二维经验模态分解和形态特征的提取；

2.在数据集的边界处划定包含5行网格数据的边界区域，由于在1的分析过程中，边界效应会对该区域带来一定的误差，因此在补偿过程中将该区域暂时去除；

3.在剩余区域中，将边界区域划分为包含5×10网格数据的补偿区域，在四个拐角处划分出四个5×5的补偿区域，用于和数据库中的数据块单元进行特征匹配；

4.按一定方向，依次对所划分出的补偿区域进行匹配补偿，具体如下：

1)选择一补偿区域，按步骤七中方法计算其强度分布范围ar、平均振幅am、极值点数en三个总体特征，并定位ar，am，2·en（对于5×10匹配区域）或4·en（对于5×5匹配区域）这三个特征在数据库中所在的单元集合；

2)在所确定的数据块单元集合中选择数据块单元，若补偿区域处于测线垂直方向，则通过旋转数据块单元，分别将其中的匹配区域1和3与该补偿区域进行形态特征对比，结果为diff₁、diff₃；若矩形区域处于测线方向，则通过旋转数据块单元，分别将其中的匹配区域2和4与该补偿区域进行形态特征对比，结果为diff₂、diff₄；若补偿区域处于四个拐角处，则通过旋转数据块单元，分别将其中的匹配区域5、6、7、8与该补偿区域进行形态对比，结果为diff₅、diff₆、diff₇、diff₈，匹配公式如式(9)，其中，M＝5，对于5×5的补偿区域，N＝5，对于5×10的补偿区域，N＝10。

3)设定匹配性指标δ＝0.05，对于2）中的每种补偿情况，若其中有一项匹配差异性小于δ（以第一种情况为例，若diff₁＜δ或diff₃＜δ），则对于该边界补偿区域，匹配过程结束，否则，选择下一数据块单元继续进行匹配。若相应单元集合中的单元均不符合要求，则选择差异性最小的进行补偿。补偿时，补偿区域内原始数值不变，数据块单元中匹配区域以外的数据补偿到边界外；

4)按以上方法，完成1轮补偿后，得到新的矩形数据集，进一步可继续进行若干轮补偿，最终得到的数据集D_c相比原始数据集D₀在最外层多出5·(n-1)层网格数据，n为总共补偿的轮数；

5.将2中暂时去除掉边界区域恢复到D_c中，针对1中的零均值平移操作，将所得数据集反向平移d₀，则补偿过程完毕。

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于地磁异常数据形态特征的边界补偿方法，流程如图1所示，包括以下几个步骤：

步骤一：对实测地磁异常数据进行规则网格化处理；

在对地磁数据进行实际测量时，受测量环境、运载器自身运动姿态等影响，测量结果无法呈现十分规则的分布，在进行数据分析和边界补偿操作之前，首先对测量数据进行规则网格化处理。

开展航空磁场测量时，测线之间间隔在0.002—0.011°（纬度）之间，测线上测点之间的间隔在0.0004—0.0008°（经度）之间。选用克里金法作为网格化插值方法，网格化后所得规则数据集D₀(x,y)在测线方向上间隔为0.001°，在测线垂直方向上的间隔为0.002°。

步骤二：平移地磁异常网格数据集，使均值为零；

计算规则网格数据集D₀(x,y)的均值M₀，将数据集在xy平面垂直方向上平移|M₀|距离，使得平移后所得数据集D(x,y)的均值M为零，从而可以消除不同数据集之间在幅值上的差异性。

步骤三：对平移后的数据集进行二维经验模态分解；

对D(x,y)进行二维经验模态分解（BEMD），具体如下，

1.初始化：r_i-1(x,y)＝D(x,y)，i＝1；

2.抽取imf_i(x,y)，即第i个本征模态函数分量（IMF）；

1)初始化：h₀(x,y)＝r_i-1(x,y)，j＝1；

2)计算h_j-1(x,y)的平均包络：

c)确定h_j-1(x,y)的极大值集合和极小值集合，并分别用径向基函数插值法构建极大值包络面E_upper(x,y)和极小值包络面E_lower(x,y)；

d)计算平均包络m_j-1(x,y)＝(E_upper(x,y)+E_lower(x,y))/2；

3)h_j(x,y)＝h_j-1(x,y)-m_j-1(x,y)，令j＝j+1；

4)一旦终止条件 $(\mathrm{SD}:=\underset{x=0}{\overset{X}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{Y}{\mathrm{\&Sigma;}}}[\frac{{|h_{j-1}(x,y)-h_j(x,y)|}^2}{h_{j-1}^2(x,y)}]<0.05),$ 满足或SD开始增大，令imf_i(x,y)＝h_j(x,y)，否则令j＝j+1，转至2）。

3.r_i(x,y)＝r_i-1(x,y)-imf_i(x,y)；

4.根据式(1)、(2)判断若r_i(x,y)在θ或其垂直方向上是否存在单调一维采样，

v_1,c(x)＝u(x,[tanθ]x+c)，   0≤θ＜π/2   (1)

$v_{2,c}\left(x\right)=\left\{\begin{array}{cc}u(\&CenterDot;,x)& \mathrm{\&theta;}=0\\ u(x,[\tan (\mathrm{\&theta;}+\mathrm{\&pi;}/2)]x+c)& 0<\mathrm{\&theta;}<\mathrm{\&pi;}/2\end{array}\right.---\left(2\right)$

其中v_1,c(x)和v_2,c(x)为r_i(x,y)在两个方向上的一维采样。如果存在，则分解过程终止，否则令i＝i+1，转至2。

经过以上分解过程，可得n个本征模态函数分量IMF₁,IMF₂,…,IMF_n和一个残余项RES。

步骤四：采用直接微分法对分解结果进行形态特征提取；

分解所得本征模态函数分量中极值点数量体现了其所包含数据形态特征的多少，且分解所得的每个本征模态函数分量中极值点数逐层递减，取IMF₁,IMF₂为特征提取分量。

采用直接微分法对本征模态函数分量进行形态特征提取的方法如下：

以IMF₁为例，设s(t)为其中一网格线数据集，确定该数据集的极大值点集合和极小值点集合，并分别利用样条函数法构建上包络线L_upper和下包络线L_lower，令A(t)表示此网格线方向上的瞬时振幅t表示x或y方向；

$A\left(t\right)=\left\{\begin{array}{cc}{L}_{\mathrm{upper}}\left(t\right)& s\left(t\right)>0\\ L_{\mathrm{lower}}\left(t\right)& s\left(t\right)<0\end{array}\right.---\left(3\right)$

令S(t)＝s(t)/A(t)＝cosφ(t)，方程两边对t求导，有

S'(t)＝-φ'(t)sinφ(t)＝-ωsinφ(t)   (4)

当sinφ(t)≠0，由ω的非负性可得，

$\mathrm{\&omega;}=-S^{\&prime;}\left(t\right)/\sin \mathrm{\&phi;}\left(t\right)=\left|S^{\&prime;}\left(t\right)\right|/\sqrt{1-S^2\left(t\right)}---\left(5\right)$

当sinφ(t)＝0，即|S(t)|＝1时，对(4)两边再次求导得，

S''(t)＝-φ''(t)sinφ(t)-[φ'(t)]²cosφ(t)＝-ω²S(t)   (6)

从而可得，

$\mathrm{\&omega;}=\sqrt{|S^{\&prime;\&prime;}\left(t\right)/S\left(t\right)|}=\sqrt{\left|S^{\&prime;\&prime;}\left(t\right)\right|}---\left(7\right)$

根据式(3)可得分量中所有数据点在其所在的两条网格线方向上的瞬时振幅特征，根据(5)或(7)可得瞬时频率特征。利用此方法，可完成对IMF₁和IMF₂中所有网格线数据集的特征提取，即得到IMF₁和IMF₂中所有数据点在横、纵双向上的瞬时幅值和瞬时频率特征。

将其余几个本征模态函数分量IMF₃,IMF₄,…,IMF_n和残余项RES直接相加作为趋势项，如式(8)，

$\mathrm{Trend}=\underset{i=3}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{IMF}}_i+\mathrm{RES}---\left(8\right)$

将Trend值作为对应数据点的趋势项幅值特征，从而所研究数据集中每个数据点具有9个特征。

步骤五：分割数据块单元，实时消除相邻单元的特征相似性冗余；

为避免步骤一至步骤四中边界效应对分析结果的影响，在零均值网格数据集中划定边界区域，并保证IMF₁和IMF₂中每条网格线数据集的两端在该边界区域内至少有一个极小值点和极大指点，将该边界区域去除，保留中心区域作为有效区域。确定数据块单元在两个方向上均包含10个数据点，采用“弓”字形分割方法进行数据块的分割，且相邻数据块的对应边界之间为1倍的网格间距，如图2所示。

采用实时方法消除相邻数据块单元之间的特征相似性冗余，当分割出一个数据块后，根据式(9)计算其与相邻已分割出数据块之间形态特征的差异性，

$\mathrm{diff}=\frac{1}{3}\&CenterDot;(\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(F_{0\mathrm{ij}}^k-F_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(F_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(A_{0\mathrm{ij}}^k-A_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(A_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{0\mathrm{ij}}^k-T_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_{0\mathrm{ij}}^k\right)}^2}})---\left(9\right)$

式中，M和N分别表示数据块单元在横向和纵向上数据点数，

和

分别表示待比较单元块与某个相邻单元块中对应数据点在IMF1和IMF2中瞬时频率特征，相应的，

和

表示幅值特征，

和

表示趋势项幅值特征，若diff＜0.05，则将新分割的数据块忽略。

步骤六：消除所有数据块单元间的特征相似性冗余；

不仅相邻数据块之间存在相似性，即使没有位置关系的数据块单元之间也会存在形态相似性，需要对这类相似性冗余进行消除。设将所有地磁异常数据集进行分割并消除相邻数据块之间冗余后共得n个数据块单元，根据图3，冗余消除过程如下，

4.将这些单元构建成一个序列U₁,U₂,…,U_n；

5.将U₁作为基准节点，根据式(9)，依次计算其与后面所有节点进行特征差异性比较，如果差异性小于0.05，则将被比较的节点舍弃，否则保留，比较完成后得到一新序列U₁,U₂,…,U_m，m≤n；

6.根据2的比较方法，依次在新序列中将U₂,U₃,…作为基准节点，不断进行筛选，当所有节点都作为基准节点被比较过之后，得到最终序列U₁,U₂,…,U_k。

步骤七：构建多层级地磁异常数据边界补偿数据库；

首先确定边界补偿时采用5行网格数据进行特征匹配，根据不同边界位置的补偿需求，按三种方式对数据块单元中进行匹配区域划分：（1）若对测线垂直方向的边界进行补偿，则将数据块单元划分为上下两个5×10的匹配区域，即区域1、3；（2）若对测线方向的边界进行补偿，则将数据块单元划分为左右两个5×10的匹配区域，即区域2、4；（3）若对数据集四个拐角处进行补偿，则将数据块单元划分为4个5×5的匹配区域，即区域5、6、7、8，如图4所示。为提高数据库的层次性和边界补偿时对数据库的搜索效率，对每个数据块单元进行强度分布范围（AR）、平均振幅（AM，在零均值前提下计算）、极值点数（EN）三个总体特征计算，AR和AM计算方式分别如式(10)、(11)所示，

AR＝max(f(x_i,y_j))-min(f(x_i,y_j))，1≤i≤m，1≤j≤n   (10)

$\mathrm{AM}=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|f(x_i,y_i)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(f(x_i,y_i)\right)/(m\&CenterDot;n)|---\left(11\right)$

按照强度分布范围、平均振幅、极值点数的顺序建立多层级数据库，如图5所示，规则如下：

4.强度分布范围按100nT作为间隔进行区间划分，即[0 100]，[100 200]，……，设分割得到的所有单元的强度分布范围可划分为m个区间范围；

5.在1所得每个区间范围的基础上，进一步按平均振幅进行区间划分，以50nT作为间隔进行区间划分，即[0 50]，[50 100]，……，设共划分n个区间范围；

6.在2所得每个区间范围的基础上，进一步按每个单元中所包含的极值点数进行区间划分，以2作为间隔进行区间划分，即[0 1]，[2 3]，……，设共划分k个区间范围。

综上，可将分割所得所有单元划分为m·n·k个单元集合。

步骤八：对待分析地磁异常数据进行边界补偿

对待分析数据进行补偿时，将边界区域作为特征匹配区域进行补偿，根据图6，具体方法如下：

6.按照步骤一至步骤四方法对待分析数据进行规则网格化、零均值平移、二维经验模态分解和形态特征的提取，并记零均值平移距离为d₀；

7.在数据集的边界处划定包含5行网格数据的边界区域，由于在1的分析过程中，边界效应会对该区域带来一定的误差，因此在补偿过程中将该区域暂时去除；

8.在剩余区域中，将边界区域划分为包含5×10网格数据的补偿区域，在四个拐角处划分出四个5×5的补偿区域，如图6所示，用于和数据库中的数据块单元进行特征匹配；

9.按一定方向，依次对所划分出的补偿区域进行匹配补偿，具体如下：

5)选择补偿区域，按步骤七中方法计算其强度分布范围ar、平均振幅am、极值点数en三个总体特征，并定位ar，am，2·en（对于5×10匹配区域）或4·en（对于5×5匹配区域）这三个特征在数据库中所在的单元集合；

6)在所确定的数据块单元集合中选择数据块单元，若补偿区域处于测线垂直方向，则通过旋转数据块单元，分别将其中的匹配区域1和3与该补偿区域进行形态特征对比，结果为diff₁、diff₃；若矩形区域处于测线方向，则通过旋转数据块单元，分别将其中的匹配区域2和4与该补偿区域进行形态特征对比，结果为diff₂、diff₄；若补偿区域处于四个拐角处，则通过旋转数据块单元，分别将其中的匹配区域5、6、7、8与该补偿区域进行形态对比，结果为diff₅、diff₆、diff₇、diff₈，匹配公式如式(9)，其中，M＝5，对于5×5的补偿区域，N＝5，对于5×10的补偿区域，N＝10。

7)设定匹配性指标δ＝0.05，对于2）中的每种补偿情况，若其中有一项匹配差异性小于δ（以第一种情况为例，若diff₁＜δ或diff₃＜δ），则对于该边界补偿区域，匹配过程结束；如果不满足指标要求，则选择下一数据块单元继续进行匹配。若相应单元集合中的单元均不符合要求，则选择差异性最小的进行补偿。补偿时，补偿区域内原始数值不变，数据块单元中匹配区域以外的数据补偿到边界外；

8)按以上方法，完成1轮补偿后，得到新的矩形数据集，进一步，可继续进行若干轮补偿，最终得到的数据集D_c相比原始数据集D₀在最外层多出5·(n-1)层网格数据，n为总共补偿的轮数；

10.将2中暂时去除掉边界区域恢复到D_c中，针对1中的零均值平移操作，将所得数据集反向平移d₀，补偿过程完毕。

Claims (6)

1.一种基于地磁异常数据形态特征的边界补偿方法，其特征在于：

（1）对实测地磁异常数据进行规则网格化处理：

选用克里金法作为网格化插值方法，网格化后所得规则网格数据集D₀(x,y)在测线方向上间隔为0.001°，在测线垂直方向上的间隔为0.002°；

（2）平移地磁异常网格数据集，使均值为零：

计算规则网格数据集D₀(x,y)的均值M₀，将数据集在xy平面垂直方向上平移|M₀|距离，使得平移后所得零均值网格数据集D(x,y)的均值M为零；

（3）对平移后的零均值网格数据集进行二维经验模态分解：

对平移后的零均值网格数据集进行二维经验模态分解，得到本征模态函数分量和一个残余项；

（4）对二维经验模态分解结果进行形态特征提取：

利用直接微分法对分解所得前两个本征模态函数分量进行双向瞬时频率和幅值特征的提取，利用其他分量进行趋势项幅值特征的计算；

（5）分割数据块单元，实时消除相邻单元的特征相似性冗余：

去除D(x,y)的边界区域，确定构建边界补偿数据库数据块单元的大小，在剩余区域中逐个进行单元的分割；实时将分割出的每一个单元与周围所有已分割出的单元块进行形态特征对比；

（6）消除所有数据块单元间的特征相似性冗余：

对所有样本数据集完成特征提取和单元块分割之后，将所有单元构成一个序列，逐个进行单元之间的形态特征差异性对比分析，逐步消除单元间的所有冗余信息；

（7）构建多层级地磁异常数据边界补偿数据库：

对最终所得的所有数据块单元计算地磁异常强度分布范围、平均振幅、极值点数三个整体特征，分别对这三个特征划分区间范围，构建多层级地磁异常数据边界补偿数据库；

（8）对待分析地磁异常数据进行边界补偿：

按照步骤（1）至步骤（4）提取待分析数据的形态特征，去除零均值网格数据集的边界部分，在剩余数据中划分边界补偿区域，对每个补偿区域，在地磁异常边界补偿数据库中提取数据块单元与其进行特征匹配，直至完成边界补偿。

2.根据权利要求1所述的一种基于地磁异常数据形态特征的边界补偿方法，其特征在于：所述的对平移后的零均值网格数据集进行二维经验模态分解包括，

（3.1）初始化：r_i-1(x,y)＝D(x,y)，i＝1；

（3.2）抽取imf_i(x,y)，即第i个本征模态函数分量；

（3.2.1）初始化：h₀(x,y)＝r_j-1(x,y)，j＝1；

（3.2.2）计算h_j-1(x,y)的平均包络：

a)确定h_j-1(x,y)的极大值集合和极小值集合，并分别用径向基函数插值法构建极大值包络面E_upper(x,y)和极小值包络面E_lower(x,y)；

b)计算平均包络m_j-1(x,y)＝(E_upper(x,y)+E_lower(x,y))/2；

（3.2.3）h_j(x,y)＝h_j-1(x,y)-m_j-1(x,y)，令j＝j+1；

（3.2.4）当满足终止条件 $\mathrm{SD}:=\underset{x=0}{\overset{X}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{Y}{\mathrm{\&Sigma;}}}\left[\frac{{|h_{j-1}(x,y)-h_j(x,y)|}^2}{h_{j-1}^2(x,y)}\right]<0.05$ 或SD值增大，令imf_i(x,y)＝h_j(x,y)，否则令j＝j+1，重新执行步骤2）；

（3.3）r_i(x,y)＝r_i-1(x,y)-imf_i(x,y)；

（3.4）判断r_i(x,y)在θ或θ的垂直方向上是否存在单调一维采样：

v_1,c(x)＝u(x,[tanθ]x+c)，   0≤θ＜π/2

$v_{2,c}\left(x\right)=\left\{\begin{array}{cc}u(\&CenterDot;,x)& \mathrm{\&theta;}=0\\ u(x,[\tan (\mathrm{\&theta;}+\mathrm{\&pi;}/2)]x+c),& 0<\mathrm{\&theta;}<\mathrm{\&pi;}/2\end{array}\right.$

其中v_1,c(x)和v_2,c(x)为r_i(x,y)在两个方向上的一维采样，如果存在，则分解过程终止，否则令i＝i+1，重新执行步骤2，直到得到n个本征模态函数分量IMF₁,IMF₂,…,IMF_n和一个残余项RES。

3.根据权利要求1或2所述的一种基于地磁异常数据形态特征的边界补偿方法，其特征在于：所述的对二维经验模态分解结果进行形态特征提取包括：

获取在IMF₁和IMF₂中所有数据点在横、纵双向上的瞬时幅值和瞬时频率特征

s(t)为IMF₁和IMF₂中的网格线数据集，确定网格线数据集的极大值点集合和极小值点集合，并分别构建上包络线L_upper和下包络线L_lower，A(t)表示此网格线方向上的瞬时振幅，t表示x或y方向；

$A\left(t\right)=\left\{\begin{array}{cc}{L}_{\mathrm{upper}}\left(t\right)& s\left(t\right)>0\\ L_{\mathrm{lower}}\left(t\right)& s\left(t\right)<0\end{array}\right.,S\left(t\right)=s\left(t\right)/A\left(t\right)=\cos \mathrm{\&phi;}\left(t\right),$

S'(t)＝-φ'(t)sinφ(t)＝-ωsinφ(t)，

当sinφ(t)≠0，由ω的非负性可得，

$\mathrm{\&omega;}=-S^{\&prime;}\left(t\right)/\sin \mathrm{\&phi;}\left(t\right)=\left|S^{\&prime;}\left(t\right)\right|/\sqrt{1-S^2\left(t\right)}$

当sinφ(t)＝0，即|S(t)|＝1时，

S''(t)＝-φ''(t)sinφ(t)-[φ'(t)]²cosφ(t)＝-ω²S(t)

$\mathrm{\&omega;}=\sqrt{|S^{\&prime;\&prime;}\left(t\right)/S\left(t\right)|}=\sqrt{\left|S^{\&prime;\&prime;}\left(t\right)\right|},$ ω为瞬时频率，

将其余本征模态函数分量IMF₃,IMF₄,…,IMF_n和残余项RES相加作为趋势项，

$\mathrm{Trend}=\underset{i=3}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{IMF}}_i+\mathrm{RES}$

将Trend值作为对应数据点的趋势项幅值特征。

4.根据权利要求3所述的一种基于地磁异常数据形态特征的边界补偿方法，其特征在于，所述的分割数据块单元，实时消除相邻单元的特征相似性冗余包括：

在零均值网格数据集中划定边界区域，分割成数据块单元，并保证IMF₁和IMF₂中每条网格线数据集的两端在边界区域内至少有一个极小值点和极大指点，将边界区域去除，保留中心区域作为有效区域，确定数据块单元在两个方向上均包含10个数据点，采用弓字形分割方式进行数据块单元的分割，且相邻数据块单元的对应边界之间为1倍的网格间距，计算其与相邻已分割出数据块之间形态特征的差异性：

$\mathrm{diff}=\frac{1}{3}\&CenterDot;(\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(F_{0\mathrm{ij}}^k-F_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(F_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{(A_{0\mathrm{ij}}^k-A_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(A_{0\mathrm{ij}}^k\right)}^2}}+\sqrt{\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(T_{0\mathrm{ij}}^k-T_{1\mathrm{ij}}^k)}^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(T_{0\mathrm{ij}}^k\right)}^2}})$

M和N分别表示数据块单元在横向和纵向上数据点数，

和分别表示待比较单元块与相邻单元块中对应数据点在IMF1和IMF2中瞬时频率特征，

和

表示幅值特征，

和

表示趋势项幅值特征，若diff＜0.05，则将新分割的数据块忽略。

5.根据权利要求1或2所述的一种基于地磁异常数据形态特征的边界补偿方法，其特征在于，所述消除所有数据块单元间的特征相似性冗余包括：

（5.1）将数据块单元构建成一个序列U₁,U₂,…,U_n；

（5.2）将U₁作为基准节点，依次与后面所有节点进行特征差异性比较，如果差异性小于0.05，则将被比较的节点舍弃，否则保留，比较完成后得到一新序列U₁,U₂,…,U_m，m≤n；

（5.3）依次在新序列中将U₂,U₃,…作为基准节点，与后面所有节点进行特征差异性比较，当所有节点都作为基准节点被比较过之后，得到最终序列U₁,U₂,…,U_k。

6.根据权利要求1或2所述的一种基于地磁异常数据形态特征的边界补偿方法，其特征在于，所述构建多层级地磁异常数据边界补偿数据库包括：

（6.1）若对测线垂直方向的边界进行补偿，则将数据块单元划分为上下两个5×10的匹配区域，即区域1、区域3；

（6.2）若对测线方向的边界进行补偿，则将数据块单元划分为左右两个5×10的匹配区域，即区域2、区域4；

（6.3）若对数据集四个拐角处进行补偿，则将数据块单元划分为4个5×5的匹配区域，即区域5、区域6、区域7、区域8，对每个数据块单元进行强度分布范围AR、平均振幅AM，、极值点数EN三个总体特征计算，

AR＝max(f(x_i,y_j))-min(f(x_i,y_j))，1≤i≤m，1≤j≤n；

$\mathrm{AM}=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|f(x_i,y_i)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(f(x_i,y_i)\right)/(m\&CenterDot;n)|$

按照强度分布范围、平均振幅、极值点数的顺序建立多层级数据库：

1）强度分布范围按100nT作为间隔进行区间划分，即[0 100]，[100 200]，……，分割得到的所有数据块单元的强度分布范围为m个区间范围；

2）在步骤1）得到的每个区间范围的基础上，按平均振幅进行区间划分，以50nT作为间隔进行区间划分，即[0 50]，[50 100]，……，划分n个区间范围；

3）在步骤2）得到的每个区间范围的基础上，按每个数据块单元中所包含的极值点数进行区间划分，以2作为间隔进行区间划分，即[0 1]，[2 3]，……，划分为k个区间范围，最终分割所得所有单元划分为m·n·k个数据块单元集合。

Images