CN110244367B - 一种基于地面多基站的ztem系统姿态补偿方法 - Google Patents

Abstract

本发明涉及电磁探测技术领域，涉及一种基于地面多基站的ZTEM系统姿态补偿方法，包括以下步骤：在地面多测点处设置基站，测量地面磁场水平分量；将测得的地面磁场水平分量应用电磁场向上延拓定理，延拓至线圈高度，得到稀疏测点x,y方向的空中磁场水平分量；将稀疏测点的空中磁场水平分量进行高精度插值，得到密集测点x,y方向的空中磁场水平分量；用插值后的稠密空中磁场水平分量数据与实测的空中磁场垂直分量建立矩阵，计算姿态补偿后的空中磁场垂直分量；根据补偿后的空中磁场垂直分量与地面磁场水平分量得到补偿后的倾子响应和补偿后的倾子散度的计算公式；将实测姿态角度数据代入到补偿后的倾子响应和补偿后的倾子散度的计算公式，进行姿态补偿。

Description

一种基于地面多基站的ZTEM系统姿态补偿方法

技术领域

本发明涉及电磁探测技术领域，尤其涉及到一种基于地面多基站的ZTEM系统姿态补偿方，属于频率域大地电磁资料处理方法。

背景技术

ZTEM是一种新型的天然源航空电磁法，该方法通过在地面设置一个基站来测量地面磁场水平分量，直升机所吊线圈沿测线飞行来测量空中磁场垂直分量，并通过倾子T将磁场的垂直分量与水平分量联系起来，从而获得地下结构信息。ZTEM航空电磁探测系统在实际飞行测量中，由于飞行速度的不稳定和风速变化，机体吊挂的接收线圈在作业过程中会发生姿态变化，使线圈所测得的磁场垂直分量中耦合部分空中磁场水平分量，影响数据反演解释。

CN103675927A公开了一种“固定翼航空电磁系统接收吊舱摆动角度的校正方法”该方法中接收线圈通过未校正的感应电动势与测量感应电动势数据间拟合误差判断校正的必要性，求出感应系数，并对测量数据除以感应系数来完成校正，该方法可以在一定程度上实现姿态校正。但是该发明未涉及到针对ZTEM系统，在地面设置多基站来测量地面磁场水平分量，在对空中线圈进行姿态补偿的方法。

嵇艳菊等，《直升机TEM系统中心回线线圈姿态校正的理论研究》地球物理学报，2010年，第1期，分析了直升机航空系统在不同姿态变化时的电磁响应计算公式，给出了不同姿态变化时的电磁响应系数，但未给出姿态校正公式。

王卫平等(2015)《频率域航空电磁系统线圈姿态变化影响及校正方法》(地球科学-中国地质大学学报，2015年7期)采用三维频率域有限差分模拟计算航空电磁系统的电磁响应，分析了不同收发线圈姿态变化对水平面和垂直轴线圈装置的电磁响应影响，根据姿态误差几何校正方法进行了电磁数据校正，有效去除了误差影响。到目前为止，仍未见用地面设置多基站的方法来对ZTEM系统线圈进行姿态补偿方面的研究。

发明内容

本发明所要解决的技术问题在于提供一种基于地面多基站的ZTEM系统姿态补偿方法，解决

本发明是这样实现的，

一种基于地面多基站的ZTEM系统姿态补偿方法，该方法以下步骤：

步骤a、在地面多测点处设置基站，测量地面磁场水平分量，其中水平第一分量表示为H_x-ground，水平第二分量表示为H_y-ground；

步骤b、将测得的地面磁场水平分量应用电磁场向上延拓定理，延拓至线圈高度，得到x,y方向的稀疏测点空中磁场水平分量分别为：H_x-air和H_y-air；

步骤c、将稀疏测点的空中磁场水平分量进行高精度插值，得到x,y方向的密集测点空中磁场水平分量：H_x-air-inter和H_y-air-inter；

步骤d、用插值后的稠密空中磁场水平分量数据与实测的空中磁场垂直分量建立矩阵，计算姿态补偿后的空中磁场垂直分量；根据补偿后的空中磁场垂直分量与地面磁场水平分量得到补偿后的倾子响应和补偿后的倾子散度的计算公式；

步骤e、将实测姿态角度数据代入到补偿后的倾子响应和补偿后的倾子散度的计算公式，进行姿态补偿。

进一步地，所述步骤a具体包括在整个测区内，设置i个基站，并覆盖整个测区，地面基站处测得的磁场水平分量分别为水平第一分量：H_x-ground(1)、H_x-ground(2)、...H_x-ground(i)，水平第二分量：H_y-ground(1)、H_y-ground(2)、H_y-ground(3)...H_y-ground(i)。

进一步地，所述步骤b具体包括：将地面磁场水平分量延拓至空中h米高度，先将地表的磁场水平分量做傅里叶变换，将它从空间域变换到波数域，应用向上延拓理论，可利用z＝0处的磁场水平分量延拓出高度为h米的空中磁场数据，h为线圈距离地面的高度。

进一步地，所述步骤c将空中稀疏测点磁场水平分量进行高精度插值包括：假设研究区域A上研究变量Z(x)，在点x_i∈A，i＝1，2，…，n处属性值为Z(x_i)，则待插点x₀∈A处的属性值Z(x₀)的克里金插值结果是已知采样点属性值Z(x_i)，i＝1，2，…，n的加权和：

其中：

Z(x_i)为第i个位置处的测量值；

λ_i为第i个位置处的测量值的未知权重；

x₀为预测位置；

n为测量值数；

根据

将稀疏测点空中磁场水平分量替代Z(x_i)求取得到插入值，与所述稀疏测点空中磁场水平分量形成密集测点空中磁场水平分量。

进一步地，用插值后得到的空中磁场水平分量，与实测空中磁场垂直分量H′_z-air建立矩阵，将建立得到的矩阵左乘姿态变化逆矩阵的方式将线圈姿态调整到水平状态，得到姿态补偿后的空中总磁场值H_compensate，姿态补偿后的空中磁场垂直分量为H_z-compensate，具体计算公式如下：

其中，H_x-air、H_y-air分别为x，y方向的空中磁场水平分量，H_x-air-inter和H_y-air-inter为密集测点磁场数据，

为线圈绕x轴旋转即发生横滚时产生的姿态角；θ为线圈绕y轴旋转即发生俯仰时产生的姿态角；ψ为线圈绕z轴旋转发生偏航时产生的姿态角。

进一步地，用补偿后的空中磁场垂直分量分别与地面磁场水分分量做比，得到补偿后的第一倾子响应T_{zx-compensate}，和第二倾子响应T_{zy-compensαte}，再分别将第一倾子响应T_{zx-compensate}，第二倾子响T_{zy-compensate}对x、y求偏导后进行加和，得到补偿后的倾子散度，具体公式如下所示。

本发明与现有技术相比，有益效果在于：利用地面多测点的磁场水平分量，应用电磁场向上延拓定理和克里金插值法，得到空中密集测点磁场数据，并对空中线圈姿态进行补偿。该方法能够在系统仅测量空中磁场垂直分量的情况下，利用地面磁场水平分量计算得到空中磁场水平分量数据，更为准确的对空中线圈姿态进行补偿，本发明采用了不同与现有技术中其他航空系统的工作原理。补偿的效果更好精确。应用仿真计算，验证该姿态补偿方法的具有很好有效性。

附图说明

图1：基于地面多基站的ZTEM系统姿态补偿方法流程图；

图2：为低阻模型示意图；

图2(a)为低阻模型x-z方向示意图；

图2(b)为低阻模型x-y方向示意图；

图3为基站位置示意图；

图4为实测姿态角度；

图5为ZTEM倾子散度等值线图；图5(a)ZTEM倾子散度实部等值线图；图5(b)ZTEM倾子散度虚部等值线图；

图6：发生姿态变化后的ZTEM倾子散度等值线图，图6(a)ZTEM倾子散度实部等值线图；图6(b)ZTEM倾子散度虚部等值线图；

图7：姿态补偿后的ZTEM倾子散度等值线图，图7(a)ZTEM倾子散度实部等值线图；图7(b)ZTEM倾子散度虚部等值线图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参见图1所示，一种基于多基站的ZTEM系统姿态补偿方法，包括以下步骤：

步骤a、在地面多测点处设置基站，测量地面磁场水平分量，其中水平第一分量表示为H_x-ground，水平第二分量表示为H_y-ground；

步骤b、将测得的地面磁场水平分量应用电磁场向上延拓定理，延拓至线圈高度，得到稀疏测点x，y方向的空中磁场水平分量分别为：H_x-air和H_y-air；

步骤c、将稀疏测点的空中磁场水平分量进行高精度插值，得到密集测点x，y方向的空中磁场水平分量：H_{x-air-int er}和H_{y-airint er}；

步骤d、用插值后的稠密空中磁场水平分量数据与实测的空中磁场垂直分量建立矩阵，推导得出变量为姿态角度的姿态补偿后的空中磁场垂直分量，根据补偿后的空中磁场垂直分量与地面磁场水分分量计算得到补偿后的倾子响应和补偿后的倾子散度计算公式；

步骤e、录入实测姿态角度数据，进行姿态补偿。

其中，步骤a，首先在整个测区内，设置i个基站，并覆盖整个测区，地面基站处测得的磁场水平分量分别为H_x-ground(1)、H_x-ground(2)、...H_x-ground(i)，H_y-ground(1)、H_y-ground(2)、H_y-ground(3)...H_y-ground(i)。

步骤b，应用电磁场向上延拓定理，将地面磁场水平分量延拓至空中h米高度，先将地表的磁场水平分量H_x-ground(x，y，z＝0)，H_y-ground(x，y，z＝0)分别做傅里叶变换，将它从空间域变换到波数域，如公式(1)(2)所示。

其中H_x(k_x，k_y，z＝0)，H_y(k_x，k_y，z＝0)为H_x-ground(x，y，z＝0)，H_y-ground(x，y，z＝0)

在波数域的表达式，k_x，k_y对应x，y的波数域变量。应用向上延拓理论，可利用z＝0处的H_x(k_ｘ，k_y，z＝0)和z＝0处的H_y(k_ｘ，k_y，z＝0)延拓出高度为h米的空中磁场数据H_x(k_x，k_y，z＝-h′)，H_y(k_ｘ，k_y，z＝-h′)。h为线圈距离地面的高度，延拓关系式如下：

将式(3)(4)做傅里叶反变换可得到延拓h米后的空中磁场水平分量h_x(k_x，k_y，z＝-h′)，h_y(k_x，k_y，z＝-h′)如公式(5)(6)所示。

地面多基站同时测量地面磁场水平分量H_x-ground(1)、H_x-ground(2)、...H_x-ground(i)，H_y-ground(1)、H_y-ground(2)、H_y-ground(3)...H_y-ground(i)。经过电磁场向上延拓后，得到的(5)和(6)式子的结果对应测点的空中磁场水平分量，表示为：H_x-air(1)、H_x-air(2)、H_x-air(3)...H_x-air(i)，H_y-air(1)、H_y-air(2)、H_y-air(3)...H_y-air(i)。

步骤c，将空中稀疏测点磁场水平分量进行高精度插值，选用克里金插值法。假设研究区域A上研究变量Z(x)，在点x_i∈A(i＝1，2，…，n)处属性值为Z(x_i)，则待插点x₀∈A处的属性值Z(x₀)的克里金插值结果是已知采样点属性值Z(x_i)(i＝1，2，…，n)的加权和，即：

Z(x_i)＝第i个位置处的测量值

λ_i＝第i个位置处的测量值的未知权重

x₀＝预测位置

n＝测量值数；

根据

将稀疏测点空中磁场水平分量替代Z(x_i)求取得到插入值，与所述稀疏测点空中磁场水平分量形成密集测点空中磁场水平分量H_x-air(3)...H_x-air(i)，H_y-air(1)、H_y-air(2)、H_y-air(3)...H_y-air(i)。

应用克里金插值法将地面稀疏测点数据插值后得到的密集测点磁场数据表示为H_x-air-inter，H_y-air-inter用于姿态补偿。

步骤d，姿态补偿公式推导，用插值后得到的空中磁场水平分量H_x-air-inter，H_y-air-inter与实测空中磁场垂直分量H′_z-air建立矩阵，为将耦合进的空中磁场水平分量进行去除，将建立得到的矩阵左乘姿态变化逆矩阵的方式将线圈姿态调整到水平状态，从而得到姿态补偿后的空中总磁场值，H_compensate。姿态补偿后的空中磁场垂直分量为H_{z-compensatｅ}，具体计算公式如下：

用补偿后的空中磁场垂直分量分别与地面磁场水分分量做比，得到补偿后的倾子响应T_{zx-compensate}，T_{zy-compensate}.再分别将T_{zx-compensate}，T_{zy-compensate}对x、y求偏导后进行加和，即得到补偿后的倾子散度，具体公式如下所示。

步骤e，在测得姿态角度的同时，代入到公式(8)和(9)得到姿态补偿后的空中总磁场值，H_compensate。姿态补偿后的空中磁场垂直分量为H_z-compensate，进行姿态补偿，判断时候计算完成，完成后输出ZTEM倾子散度数据。

采用仿真模型对本发明的方法进行仿真：

首先，设计一个仿真模型，如图2所示，模型为一个低阻异常体，大小为2200*2300*600(m³)，测区大小为16000*16000(m²)，在测区上覆盖多个测点，本实例中选择安置25个基站来测量地面磁场水平分量，基站位置如图3所示。基站处测得的磁场水平分量数据为：H_x-ground(1)、H_x-ground(2)、H_x-ground(3)...H_x-grouna(25)，H_y-ground(1)、H_y-ground(2)、...H_y-ground(25)。

其次，应用电磁场向上延拓定理，将地面磁场水平分量延拓至空中80米高度，先将地表的磁场水平分量(H_x-ground(x，y，z＝0)，H_y-ground(x，y，z＝0))做傅里叶变换，将它从空间域变换到波数域，如公式(1)(2)所示。

其中H_x(k_x，k_y，z＝0)H_y(k_x，k_y，z＝0)为(H_x-ground(x，y，z＝0)，H_y-ground(x，y，z＝0))在波数域的表达式，k_x，k_y对应x，y的波数域变量。应用向上延拓理论，可利用z＝0处的H_x(k_x，k_y，z＝0)和z＝0处的H_y(k_x，k_y，z＝0)延拓出高度为80米的空中磁场数据，H_x(k_x，k_y，z＝-h′)，H_y(k_x，k_y，z＝-h′)。延拓关系式如下公式(3)(4)所示：

将式(3)(4)做傅里叶反变换即可得到延拓80米后的空中磁场水平分量h_x(k_x，k_y，z＝-h′)，hy(k_x，k_y，z＝-h′)，如公式(5)(6)所示。

地面多基站同时测量地面磁场水平分量H_x-ground(1)、H_x-ground(2)、H_x-ground(3)...H_x-ground(25)，H_y-ground(1)、H_y-ground(2)、H_y-ground(3)...H_y-ground(25)。应用电磁场向上延拓后定理得到对应测点的的空中磁场值为H_x-air(1)、H_x-air(2)、H_x-air(3)...H_x-air(25)，H_y-air(1)、H_y-air(2)、H_y-air(3)...H_y-air(25)。

第三步，将延拓后的空中磁场数据进行高精度插值，插值方法选用克里金插值法。假设研究区域A上研究变量Z(x)，在点x_i∈A(i＝1，2，…，n)处属性值为Z(x_i)，则待插点x₀∈A处的属性值Z(x₀)的克里金插值结果是已知采样点属性值Z(x_i)(i＝1，2，…，n)的加权和，即：

Z(x_i)＝第i个位置处的测量值；

λ_i＝第i个位置处的测量值的未知权重；

x₀＝预测位置；

n＝测量值数；

应用克里金插值法将地面稀疏测点数据插值后得到的密集测点磁场数据、H_x-air-inter，H_y-air-inter用于姿态补偿。

第四步，当飞机所吊线圈由于机体的不稳定及飞行速度的变化或天气等因素的影响而产生姿态变化时，线圈接收到的空中磁场垂直分量中会耦合部分空中磁场水平分量，引入测量误差，发生姿态变化时空中磁场垂直分量大小如公式(8)所示：

其中

为线圈绕x轴旋转即发生横滚(roll)时产生的姿态角；θ为线圈绕y轴旋转即发生俯仰(pitch)时产生的姿态角；ψ为线圈绕z轴旋转即发生偏航(yaw)时产生的姿态角；用插值后得到的空中磁场水平分量与发生姿态变化后线圈测得的空中磁场垂直分量建立矩阵，并左乘旋转矩阵的逆矩阵的方式将对线圈姿态进行补偿，得到补偿后的空中总磁场为H_compensate，补偿后的空中磁场垂直分量为H_z-compensate，计算公式如(9)(10)所示。用补偿后的空中磁场垂直分量带入公式即可求得补偿后的ZTEM倾子与倾子散度，如公式(10)(11)所示：

第五步，如图4所示，为一组实测从线圈姿态角度，角度变化范围如图所示，其中横滚和俯仰的姿态变化角度

θ的范围大概-10～10度，偏航的姿态变化角度ψ的范围大概在-15～15度。当无姿态变化时，正演计算得到ZTEM倾子散度的等值线图如图5所示，(a)为散度实部，(b)为虚部。从图中看出，当无姿态变化时，倾子散度能够很直观的反映地下异常体的大小及位置。图6为发生姿态变化后，正演计算得到的倾子散度等值线图，(a)ZTEM倾子散度实部等值线图；(b)ZTEM倾子散度虚部等值线图，从图中看出，姿态变化对ZTEM系统影响很大，异常体完全淹没在背景中，不能识别出其大小及位置。图7为用该补偿方法对线圈进行姿态补偿后，(a)ZTEM倾子散度实部等值线图；图(b)ZTEM倾子散度虚部等值线图，计算得到的ZTEM倾子散度等值线图，从图中可以看出，无论是散度实部还是虚部等值线图都能够清晰的识别出异常体的大小及位置，本发明的有效性得以验证。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于地面多基站的ZTEM系统姿态补偿方法，其特征在于，包括以下步骤：

步骤a、在地面多测点处设置基站，测量地面磁场水平分量，其中水平第一分量表示为H_x-ground，水平第二分量表示为H_y-ground；

步骤b、将测得的地面磁场水平分量应用电磁场向上延拓定理，延拓至线圈高度，得到x,y方向的稀疏测点的空中磁场水平分量：H_x-air和H_y-air；

步骤c、将稀疏测点的空中磁场水平分量进行高精度插值，得到x,y方向的密集测点空中磁场水平分量：H_x-air-inter和H_y-air-inter；

所述步骤c将稀疏测点的空中磁场水平分量进行高精度插值包括：假设研究区域A上研究变量Z(x)在点x_i处属性值为Z(x_i)，则待插点x₀处的属性值Z(x₀)的克里金插值结果为已知采样点属性值Z(x_i)的加权和：

其中：

Z(x_i)为研究变量Z(x)在点x_i处属性值，x_i∈A,i＝1,2,...,n；

λ_i为第i个位置处的测量值的未知权重；

x₀∈A；

n为测量值数；

根据公式

用稀疏测点的空中磁场水平分量替代Z(x_i)处取得的插值结果，形成密集测点空中磁场水平分量；

步骤d、用插值后的密集测点空中磁场水平分量与实测的空中磁场垂直分量建立矩阵，计算姿态补偿后的空中磁场垂直分量；根据补偿后的空中磁场垂直分量与地面磁场水平分量得到补偿后的倾子响应和补偿后的倾子散度的计算公式；

步骤e、将实测姿态角度数据代入到补偿后的倾子响应和补偿后的倾子散度的计算公式，进行姿态补偿。

2.按照权利要求1所述的方法，其特征在于，所述步骤a具体包括：在整个测区内设置i个基站，并覆盖整个测区，基站处测得的磁场水平分量分别为水平第一分量：H_x-ground(1)、H_x-ground(2)、...H_x-ground(i),水平第二分量：H_y-ground(1)、H_y-ground(2)、H_y-ground(3)…H_y-ground(i)。

3.按照权利要求1所述的方法，其特征在于，所述步骤b具体包括：将地面磁场水平分量做傅里叶变换，将它从空间域变换到波数域，然后，根据电磁场向上延拓理论，将z＝0处的磁场水平分量延拓出线圈距离地面的高度为h米的空中磁场数据。

4.按照权利要求1所述的方法，其特征在于，用插值后得到的密集测点空中磁场水平分量，与实测的空中磁场垂直分量H′_z-air建立矩阵，通过将建立得到的矩阵左乘姿态变化逆矩阵，把线圈姿态调整到水平状态，得到姿态补偿后的空中总磁场值H_compensate；补偿后的空中磁场垂直分量为H_z-compensate，具体计算公式如下：

其中，

为线圈绕x轴旋转即发生横滚时产生的姿态角；θ为线圈绕y轴旋转即发生俯仰时产生的姿态角；ψ为线圈绕z轴旋转发生偏航时产生的姿态角。

5.按照权利要求4所述的方法，其特征在于，用补偿后的空中磁场垂直分量分别与地面磁场水平分量做比，得到补偿后的第一倾子响应T_{zx-compensate}和第二倾子响应T_{zy-compensate}，再分别将第一倾子响应T_{zx-compensate}和第二倾子响应T_{zy-compensate}对x、y求偏导后进行加权求和，得到补偿后的倾子散度，具体公式如下所示：

Images