航空重力梯度仪实时重力梯度解调方法

图1 旋转加速度计重力梯度仪测量原理

Fig.1 Schematic of the rotating accelerometer gravity gradiometer

2 重力梯度解调相位角分析

重力梯度仪样机是由复杂机械结构、电气线路和控制算法等构成的高精密测量仪器。由于重力梯度仪所处的环境非常恶劣,如果直接对输出信号进行重力梯度解调,必然会引入额外的干扰,导致重力梯度解调过程中造成数据的不准确,因此,在进行重力梯度信号解调前,需要实时确定出精确的解调相位角。重力梯度解调相位角结构如图2所示。

图2

图2 重力梯度解调相位角结构

Fig.2 Schematic of the gravity gradient demodulation phase angle

4只加速度计输出信号经过信号调理处理,得到模拟信号 $E_{o u t} (t)$ ,角度信息转换器输出同步脉冲信号给A/D转换器,A/D转换器对模拟信号进行A/D转换得到第(i-1)Δt时刻数字信号 $E_{o u t} (i)$ ,其表达式为^[9]:

(1)

\{\begin{array}{l} E_{o u t} (i) = A_{k} s i n [(2 k + 1) (ω Δ t (i - 1) + α - θ)] + \\ B_{k} c o s [(2 k + 1) (ω Δ t (i - 1) + α - θ)] + \\ C_{k} s i n [(2 k + 2) (ω Δ t (i - 1) + α - θ)] + D \overset{\cdot}{ω} \\ α = φ - ψ \end{array}

式中:A_k、B_k分别为2k+1次谐波正弦、余弦信号幅值;C_k为2k+2次谐波信号幅值;D为圆盘旋转角加速度幅值; $α$ 为重力梯度解调相位角; $φ$ 为调理通路信号延迟导致的延时角; $φ$ 为加速度计A₁与重力梯度测量坐标系x轴之间的初始夹角; $θ$ 为环境物体引起的方位角;Δt为采样时间; $ω$ 为重力梯度仪圆盘旋转角频率; $\overset{\cdot}{ω}$ 为重力梯度仪圆盘旋转角加速度。

参考信号源数字正余弦信号,它由DSP发出^[9],其表达式为:

(2)

\{\begin{array}{l} S_{s} (i, j) = s i n 2 [ω Δ t (i - 1) + ϑ_{j}] \\ S_{c} (i, j) = c o s 2 [ω Δ t (i - 1) + ζ_{j}] \end{array}

式中: $ϑ_{j}$ 、 $ζ_{j}$ 分别为正弦、余弦信号参考源的初始相角在第j次的修正角。

利用参考信号源数字正余弦信号对A/D转换器的输出信号进行2倍频幅值解调,通过改变其初始相角直到得到最大的幅值为止,最终的初始相角即为重力梯度解调相位角^[9],其计算表达式为:

(3)

\{\begin{array}{l} ϑ_{q} = ϑ_{i} |m a x \overset{N q}{\underset{k = N q - N + 1}{Σ}} E_{o u t} (k) S_{s} (k, i) \\ ζ_{q} = ζ_{i} |m a x \overset{N q}{\underset{k = N q - N + 1}{Σ}} E_{o u t} (k) S_{c} (k, i) \end{array}

式中: N为圆盘每旋转一圈角度信息转换器输出的脉冲数, $ϑ_{q}$ 、 $ζ_{q}$ 分别为第q次得到的重力梯度解调相位角;对得到的重力梯度解调相位角 $ϑ_{q}$ 、 $ζ_{q}$ 进行组合,然后把组合的数据通过零相位低通滤波器进行滤波处理,最后对处理后的相位角进行均值处理,最终得到重力梯度解调相位角 $α$ 。

3 重力梯度实时解调分析

为了实现航空重力梯度实时测量,在求解出解调相位角的基础上,添加重力梯度解调部分^[10],结构图如图3所示。

图3

图3 航空重力梯度仪实时重力梯度解调结构

Fig.3 Structure of the real-time gravity gradient demodulation for the airborne gravity gradiometer

将4只加速度计输出信号经过信号调理和滤波电路进行处理,得到模拟信号 $E_{o u t} (t)$ ,角度信息转换器输出同步脉冲信号给A/D转换器,A/D转换器对模拟信号 $E_{o u t} (t)$ 进行A/D转换得到第(k-1)Δt时刻数字信号 $E_{o u t} (t)$ ^[10],其表达式为

(4)

E_{o u t} (k) = 2 K_{I} K R \{(Γ_{y y} - Γ_{x x}) s i n 2 [ω Δ t (k - 1) - α] + 2 α_{x y} c o s 2 [ω Δ t (k - 1) - α]\}

式中:k为数据序列索引; K为系统信号放大增益; Δt为采样时间; K_I为加速度计标度因数; R为圆盘半径; $ω$ 为圆盘旋转角频率; $α$ 为重力梯度解调相位角; Г_xx为x轴上的重力加速度分量在x轴方向上的空间导数; Г_yy为y轴上的重力加速度分量在y轴方向上的空间导数; (Г_yy - Г_xx)为重力梯度Г_yy与Г_xx之差; Г_xy为x(y)轴上的重力加速度分量在y(x)轴方向上的空间导数。

圆盘光栅编码器上的角度信息转换器实时输出圆盘角位置信息为 $θ_{p, k}$ ,利用角位置信息 $θ_{p, k}$ 对重力梯度数字信号进行实时重力梯度解调,得到具有时空信息的重力梯度数据,重力梯度解调表达式为:

(5)

\{\begin{array}{l} (Г_{y y} - Г_{x x}) (g p s, Δ t p k, p, k) = \frac{\overset{M}{\underset{p = 1}{Σ}} [\overset{N + p - 1}{\underset{k = p}{Σ}} E (k) s i n [2 ω Δ t (k - p) - α + 2 θ_{p, k}]]}{M N R K K_{I}} \\ Г_{x y} (g p s, Δ t p k, p, k) = \frac{\overset{M}{\underset{p = 1}{Σ}} [\overset{N + p - 1}{\underset{k = p}{Σ}} E (k) c o s [2 ω Δ t (k - p) - α + 2 θ_{p, k}]]}{2 M N R K K_{I}} \end{array}

式中: $q$ 为当前时刻的圆盘旋转圈数序列索引; M为圆盘旋转总圈数; gps为GPS提供的位置信息。

4 仿真验证

根据旋转加速度计重力梯度仪工作原理,采用高性能计算机、可编程高精度电流源、低噪声电流放大器、低噪声电压放大器、多路切换开关和高精度数字电压表等仪器组建了重力梯度半实物仿真分析系统。计算机把仿真数据发送给仿真分析系统,转换为接近真实的重力梯度仪工作环境下的重力梯度信号,通过改变各种实验参数和实验数据,可以全方位模拟重力梯度信号,有助于进一步对重力梯度仪的信号特征、误差分析、加速度计标度因子反馈调整算法以及信号处理方法等进行详细研究 ${^{[}}^{11]}$ ,重力梯度半实物仿真分析系统实物图如图4所示。

图4

图4 重力梯度仪半物理仿真系统实物

Fig.4 Hardware-in-the-loop simulation platform of the gravity gradiometer

设定圆盘基线距离为1 m,加速度计标度因数为32 mA/g,重力加速度为9.81 m/s²,圆盘旋转周期为20 s,采样率为2 Hz,初始时刻,加速度计A₁与x轴之间的夹角为0°,环境物体引起的方位角为0°,对重力梯度信号进行数据采集并处理,角度变动步长为0.01°,解调相位角结果如图5所示。从图5中可以看到,重力梯度解调相位角在-50.3°左右,零相位滤波处理前的解调相位角波动较大,标准差为0.14°,通过零相位滤波处理后降低了波动范围,标准差变为0.046°,最后做均值处理,得到重力梯度解调相位角均值为-50.29°。如果采用0°初始相位角作为重力梯度解调相位角,则会导致较大的重力梯度测量误差,从而影响重力梯度测量精度^[10]。

图5

图5 重力梯度解调相位角在零相位滤波前后比较波形

Fig.5 Comparison of gravity gradient demodulation phase angle before and after zero-phase filtering

通过计算得到该环境物体引起的重力梯度分量(Г_yy -Г_xx)和Г_xy的理论值分别为-335 Eu和0 Eu,如果不做重力梯度解调相位角处理,采用0°初始相位角作为重力梯度解调相位角,解调得到的重力梯度分量(Г_yy -Г_xx)和Г_xy分别为61.89 Eu和164.62 Eu,与理论梯度值之间的误差达到了几百Eu,因此直接对原始重力梯度数据解调会产生非常大的测量误差,无法进行航空重力梯度测量;若采用计算得到的解调相位角作为重力梯度解调初始相位角对重力梯度信号进行解调,解调得到重力梯度分量(Г_yy -Г_xx)和Г_xy分别为-335.008 5 Eu和0.000 55 Eu,与理论梯度值之间的误差小于0.01 Eu,能够满足航空重力梯度测量需求,重力梯度解调相位补偿前后的梯度信息比较如图6、图7所示。

图6

图6 相位补偿前后的梯度信息(Г_yy-Г_xx)比较的波形

Fig.6 Comparison of the gradient component (Г_yy-Г_xx) before and after phase compensation

图7

图7 相位补偿前后的梯度信息Г_xy比较的波形

Fig.7 Comparison of the gradient component Г_xy before and after phase compensation

5 结论

本文针对航空重力梯度仪在测量过程中解调相位角对重力梯度测量精度的影响,提出了实时确定解调相位角和重力梯度解调求解的方法。首先提出了解调相位角的确定方法和求解算式;然后给出了实时重力梯度解调算式;接着利用解调相位角对重力梯度仪输出数据进行实时重力梯度解调,得到更高精度的重力梯度数据;最后为了验证所提方法的有效性,在重力梯度半物理仿真系统上进行半物理仿真实验测试,仿真结果表明,所给的初始相位角解算式和重力梯度实时解调算式能够明显提高重力梯度测量精度,可以为重力梯度仪的研制提供技术参考。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

钱学武, 蔡体菁, 尹航.

一种旋转加速度计重力梯度仪重力梯度解调方法

[J]. 物探与化探, 2015, 39(S1):80-83.

Qian

X W

, Cai

T J

, Yin

A gravity gradient demodulation method for rotating accelerometer gravity gradiometer

[J]. Geophysical and Geochemical Exploration, 2015, 39(S1):80-83.

[2]

蔡体菁, 钱学武, 丁昊.

旋转加速度计重力梯度仪重力梯度信号仿真

[J]. 物探与化探, 2015, 39(S1):80-83.

Cai

T J

, Qian

X W

, Ding

Signal simulation of gravity gradient for gravity gradiometer of rotating accelerometer

[J]. Geophysical and Geochemical Exploration, 2015, 39(S1): 80-83.

[3]

Difrancesco

Advances and challenges in the development and deployment of gravity gradiometer systems

[C]// EGM 2007 International Workshop, Innovation in EM, Grav. and Mag.Methods:A new Perspective for Exploration, Italy: Capri, 2007.

[4]

Roberts

, Chowdhury

P R

, Lowe

S J

, et al.

Airborne gravity gradiometer surveying of petroleum systems under Lake Tanganyika, Tanzania

[J]. Exploration Geophysics, 2016, 47(3):228-236.

DOI:10.1071/EG15075 URL [本文引用: 1]

[5]

DiFrancesco

, Meyer

, Christensen

, et al.

Gravity gradiometry-today and tomorrow

[C]// 11^th SAGA Biennial Technical Meeting and Exhibition, 2009.

[6]

Lococo

J J

The Falcon airborne gravity gradiometer for engineering applications

[C]// Environment and Engineering Geophysical Society, 2010.

[7]

Yamanaka

HeliFalcon© airborne gravity gradiometer data acquisition in rugged terrain

[C]// Global Meeting Abstracts Society of Exploration Geophysicists, 2013.

[8]

钱学武, 蔡体菁.

旋转加速度计重力梯度仪数据处理方法

[J]. 东南大学学报:自然科学版, 2016, 46(4):708-712.

Qian

X W

, Cai

T J

Data processing methods for rotating accelerometer gravity gradiometer

[J]. Journal of Southeast University Natural Science Edition, 2016, 46(4): 708-712.

[9]

钱学武, 赵立业.

旋转加速度计重力梯度仪重力梯度解调相位角确定方法及装置

[P]. 中国专利, CN 111624671A. 2020.09.04

Qian

X W

, Zhao

L Y

Method and device for determining phase angle of gravity gradient demodulation for rotating accelerometer gravity gradiometer

[P]. Patent, CN 111624671A. 2020.09.04.

[10]

钱学武, 赵立业.

一种航空重力梯度仪实时重力梯度解调方法及装置

[P]. 中国专利, CN 111650664A. 2020.09.11.

Qian

X W

, Zhao

L Y

Method and device of real-time gravity gradient demodulation for an airborne gravity gradiometer

[P]. Patent, CN 111650664A. 2020.09.11.