时间域航空电磁数据约束多项式拟合调平方法与流程

本发明涉及一种时间域航空电磁资料处理方法，尤其是带约束的多项式拟合航空电磁数据调平方法。

背景技术：

时间域航空电磁法(ATEM)是一种以飞机为载体的地球物理探测方法，广泛用于矿产资源、地下水及油气探测等方面，具有速度快、效率高等优点。但是航空电磁系统在飞行探测过程中，由于受到诸多的不可避免的干扰因素，如飞行方向、飞行高度、温度等，致使测区的各测线数据基准处于不同的水平上，在平面等值线图上呈现条带状误差，区域与局部异常被扭曲、畸变，产生虚假异常，数据质量变差。为了去除条带状误差，得到准确的异常分布信息，需要对航空电磁数据进行调平(leveling)。

CN105068131A公开了“一种航磁数据调平方法”。该方法沿测线方向和穿越测线方向对航磁数据进行总体平均经验模态分解，得到若干个子数据块，然后在该若干个子数据块中找出干扰子数据块并消除，将其余的子数据块相加得到调平后的航磁数据。该方法通过可视化操作找出调平误差所在，但是并未涉及航空电磁数据各时间道之间的关联约束。

航空数据调平方法在磁数据和电磁数据方面均有应用。Huang(2008)基于航空电磁数据各测线间相关性，以各时间道数据为一个整体，选取数据基准，依次拟合相邻测线调平误差，取得较好的效果；Majid等(2010)提出在一维和二维窗中进行多项式拟合，通过移动窗对数据进行逐点调平，该调平方法能更好的保留地质特征，且通用于规则测线和不规则测线的调平。到目前为止，仍未见航空电磁数据各时间道相互约束的调平方法。所以，本发明提出一种带约束的航空电磁数据调平方法。

技术实现要素：

本文的发明目的在于提供一种准确的时间域航空电磁系统带约束的多项式拟合调平方法，该方法主要针对传统拟合方法忽略各时间道间关系，可能导致调平结果各道之间扭曲的问题。首先由于航空电磁数据的高度敏感性，去除电磁数据的高度影响；其次对于各时间道数据，选择基准测线，依次拟合相邻测线调平误差，得到一次调平结果；最后取一次调平结果，选择基准道，通过各测点的衰减关系建立各时间道间的不等式关系，以此不等式为约束，拟合相邻测线调平误差，得到调平结果。

本方法不仅考虑了各测线间的连续性，而且加入了各时间道间的相关性，使整体调平过程更加的完整。

时间域航空电磁数据约束多项式调平方法，包括以下步骤：

a、录入航空电磁数据；

b、去除数据高度影响；

c、选择基准测线，对各道数据分别进行一次调平；

d、输出一次调平结果；

e、判断一次调平结果各测线不符合衰减规律测点数量是否达到阈值，否，进行下一步；

f、选择基准道，对不符合衰减规律测线进行二次调平；

g、判断调平结果，达到阈值要求，否，返回到上一步；

h、是，把调平结果转换到原数据等级；

i、输出数据调平结果。

步骤a所述的录入时间域航空电磁数据，基于叠加偶极子假设，步骤b所述的为去除原始数据d_t高度影响，将系统接收的瞬态数据变换为数据d

d＝d_t×(h/R)³ (1)

其中，h为飞机飞行高度，R为发射线圈半径。若电磁探测系统发射线圈非圆形，取等效发射线圈半径其中，S为探测系统发射线圈半径。

步骤c所述的基准测线选择，根据数据的平面等值线图，在无明显条带状处，选择基准测线，假设选择测区数据d的第r条测线为基准测线，则基准测线表示为：其中，N表示各测线测点数，M表示各测点的时间道数，各测点表示为x＝(1,2,…,N)，一次调平是对各时间道数据进行独立调平，以第m道为例，基准测线d^r和其相邻测线d^r+1的第m道数据分别可以表示为和

选择测线中包含调平误差f(x)为多项式形式，即，

f(x)＝a₀+a₁x+a₂x²+…+a_kx^k， (2)

其中，a_k为多项式系数，k为多项式拟合阶数，相邻测线的调平结果可获得：

${\mathrm{dl}}_m^{r+1}=d_m^{r+1}-f\left(x\right),---\left(3\right)$

由于相邻测线间数据的相关性，存在

$\left|\right|{\mathrm{dl}}_m^{r+1}-d_m^r|{|}^2=||d_m^{r+1}-f\left(x\right)-d_m^r|{|}^2=\min ---\left(4\right)$

联合式(2)和式(4)，解最小二乘方程，得到调平误差系数a＝(a₀,a₁,…,a_k)^T，代入式(3)，得到调平结果循环取出基准测线d^r和其相邻测线d^r+1的第m(m＝1,2,...,M)道数据，代入公式(2)～(4)得到相邻测线所有道调平结果dl^r+1，以dl^r+1为基准测线，依照此步骤调平其相邻测线d^r+2，从而调平整个测区数据，得到一次调平结果dl。

步骤e判断一次调平结果各测线不符合衰减规律测点数量是否达到阈值，设定阈值为测点数的p倍(N*p,p＜1)，判断一次调平结果dl各测线不符合衰减规律测点数是否小于阈值。若大于设定阈值，对此测线数据进行二次调平，

根据初始数据和一次调平结果衰减曲线，选择较为符合衰减规律的道为基准道，基准道不需要进行二次调平，直接取一次调平结果，以第r+1条测线为例，假设选择各测点的第m-1道作为基准道，则其相邻第m道的二次调平结果数据应该满足：

${\mathrm{dl}}_{m-1}^{r+1}\≥{\mathrm{dls}}_m^{r+1}---\left(5\right)$

同时，该测线和其基准测线(第r条测线)的第m道数据还需要满足：

$\left|\right|{\mathrm{dls}}_m^{r+1}-{\mathrm{dl}}_m^r|{|}^2=||{\mathrm{dl}}_m^{r+1}-f\left(x\right)-{\mathrm{dl}}_m^r|{|}^2=\min ---\left(6\right)$

其中，式中的调平误差f(x)仍取多项式形式，联合式(2)、式(5)和式(6)解得调平误差系数a＝(a₀,a₁,…,a_k)^T，则第r+1条测线第m道数据的二次调平结果为：

${\mathrm{dls}}_m^{r+1}={\mathrm{dl}}_m^{r+1}-f\left(x\right)---\left(7\right)$

此时，第m道可以作为新的基准道去约束其相邻道(第m+1道)，依次循环步骤e得到测区二次调平结果dls。

5、按照权利要求4所述的时间域航空电磁数据约束多项式拟合调平方法，其特征在于，判断二次调平结果dls各测线不符合衰减规律测点数是否小于阈值，如不满足，重步骤e，如满足，取式(1)的逆变换，把此结果转换到原数据等级，结合步骤a取逆变换，得到调平后瞬态数据其中为被调平测点的平均飞行高度。

有益效果：本发明公开的时间域航空电磁数据约束多项式拟合调平方法，(1)可以快速、有效的去除测区数据的条带状误差，如图2(c)所示；(2)本发明由于在调平过程中加入了测点各时间道衰减关系的不等式约束，调平结果的剖面曲线得到很大的改善，避免调平结果各道数据之间的扭曲，如图3(c)和3(f)所示，保证了调平结果的准确性；(3)本发明不需要测线切割线的辅助校正，节省了大笔费用。从图中可以看出，由于约束多项式拟合调平加了各道之间的约束，其调平结果优于多项式拟合一次调平结果。

附图说明

图1时间域航空电磁数据约束多项式拟合调平方法流程图。

图2是测线L310数据平面等值线图。

图2a测线L310初始数据；

图2b测线L310多项式拟合一次调平结果；

图2c测线L310约束多项式拟合调平结果。

图3是测线L310剖面曲线。

图3a测线L310初始数据；

图3b测线L310多项式拟合一次调平结果；

图3c测线L310约束多项式拟合调平结果；

图3d测线L310测道(8-19)测点(4200-4600)初始数据；

图3e测线L310测道(8-19)测点(4200-4600)多项式拟合一次调平结果；

图3f测线L310测道(8-19)测点(4200-4600)约束多项式拟合调平结果。

图4是测线L310测点4400衰减曲线对比。

具体实施方式

下面结合附图和实施例作进一步的详细说明：

以一组Geotech公司在加拿大安大略省西北部的内斯特瀑布区域探测的一组数据为例，该数据于2014年1月26号至3月1号使用Plus系统探测得到。取由西至东测线号为L310-L720的42条测线为基准测线，每条测线对应插值为5748个测点，每个测点抽道为32道数据。

步骤一 读入航空电磁数据B((pV*ms)/(A*m⁴))，记为d_t，读入各测点数据对应高度，把原始数据d_t去除高度影响，转换为数据

d＝d_t×(h/R)³， (1)

其中，发射线圈半径R为26m。

步骤二 数据d的第一道数据(测线L310)等值线图如图2(a)所示(箭头方向为测线方向)，可以看出这组数据中包含有明显的条带状调平误差，其测线L310剖面曲线如图3(a)所示。各测点可以表示为x＝(1,2,…,5748)。根据图2(a)，选择第6条测线(测线L360)为基准测线，选择调平误差方程为线性函数，即：

f(x)＝a₀+a₁x， (2)

其中，a＝(a₀,a₁)^T为待求调平误差系数。以第一道数据为例，基准测线d⁶和其相邻测线d⁷的第一道数据分别可以表示为和

由于相邻测线间数据的相关性，存在：

$\left|\right|{\mathrm{dl}}_1^7-d_1^6|{|}^2=||d_1^7-f\left(x\right)-d_1^6|{|}^2=\min .---\left(3\right)$

相邻测线的调平结果可获得：

${\mathrm{dl}}_1^7=d_1^7-f\left(x\right).---\left(4\right)$

联合式(2)和式(3)，解最小二乘方程，得到调平误差系数a，代入式(4)，得到第7条测线第一道调平结果对于此测线其他各道数据，重复此步骤，得到相邻测线调平结果dl⁷。以dl⁷为基准测线，调平其相邻测线d⁸，以此规律继续，调平整个测区数据，得到一次调平结果dl，一次调平结果如图2(b)，其测线L310剖面曲线如图3(b)所示。

步骤三设定阈值为测点数的1％倍，即1％*5748＝58，判断一次调平结果dl各测线不符合衰减规律测点数是否小于阈值。若大于设定阈值，对此测线数据进行二次调平。

选择各测线的第5道作为基准道，以第7条测线为例，则基准道的相邻道(第6道)的二次调平结果数据应该满足：

${\mathrm{dl}}_5^7\≥{\mathrm{dls}}_6^7.---\left(5\right)$

同时，此测线和其基准测线(第6条测线)的第6道数据还需要满足：

$\left|\right|{\mathrm{dls}}_6^7-{\mathrm{dl}}_6^6|{|}^2=||{\mathrm{dl}}_6^7-f\left(x\right)-{\mathrm{dl}}_6^6|{|}^2=\min ---\left(6\right)$

其中，式中的调平误差f(x)仍取线性多项式形式(如式(2))。联合式(2)、式(5)和式(6)解得调平误差系数a＝(a₀,a₁)^T，则第7条测线第6道数据的二次调平结果为：

${\mathrm{dls}}_6^7={\mathrm{dl}}_6^7-f\left(x\right).---\left(7\right)$

此时，第6道可以作为新的基准道去约束其相邻道(第7道)，依次规律得到测区二次调平结果dls，调平结果如图2(c)所示，其测线L310剖面曲线如图3(c)所示。

步骤四判断二次调平结果dls各测线不符合衰减规律测点数是否小于阈值，本实例数据满足，取式(1)的逆变换，把此结果转换到原数据等级。结合步骤一，得到调平后瞬态数据其中为被调平测点的平均飞行高度。