2) 中国天津 300161 中国地震局第一监测中心
2) The First Monitoring and Application Center, China Earthquake Administration, Tianjin 300161, China
磁通门经纬仪(也称DI仪)是一种用于测量地磁偏角D和地磁倾角I绝对值的仪器,它由无磁经纬仪、磁通门传感器、磁通门检测器3部分组成,可用于固定地磁台站测量和野外流动测量,是目前使用最广泛的地磁绝对测量仪器之一(Jankowski et al,1996;宋思璇等,2020)。
根据目前的研究结果,磁通门经纬仪测量误差来源主要有:修正系数(格值)、温度、零点偏移、分辨力,还有与之同步观测的磁通门磁力仪的X、Y、Z三轴正交度、灵敏度(姚远等,2016;王晓美等,2017;张涛等,2018;刘浩等,2022)、偏置量、D分量定向、观测墩漂移、台站供电系统稳定性等(屈文斌,2020;车濛琪等,2020;罗玉芬等,2021;陈贤等,2021;迟铖等,2021)。Marsal等(2007)提出, 日常观测实践中系统效应会引起误差,并分析了磁通门经纬仪由机械结构引入的不确定度,该问题可通过目前已广泛使用的“四位置”测量方式解决。Jankowski等(1996)对单台仪器测量时的影响量进行了分析,但仅在定性层面进行了概述。张策等(2020)对自动磁偏角磁倾角测量仪的误差进行了分析研究,与本文所讨论的人工仪器有一定差别。
所有台站的磁通门经纬仪需要定时进行比对测量(简称比测),其中最重要的指标是“仪器差”,即被测试仪器与标准仪器测量磁偏角和磁倾角结果的差值,属于一种测量误差,但是,针对该指标的定量误差分析研究或者不确定度评定方法,目前尚属空白,所以比测的结果也无相应的误差或者不确定度表示。这导致我国地磁台网对仪器指标进行合格性判定时经常出现错误,尤其当结果在最大允许误差限附近时更容易误判。针对类似问题,解决办法就是对测量结果引入不确定度。传统的磁通门经纬仪比测模型只体现了理想状况,没有综合考虑测量过程中的多种系统因素,并不适用于开展不确定度评定。
综上,本文将建立包含系统因素的磁通门经纬仪比测模型,并通过对测量结果进行不确定度评定、合格判定及数据比较来验证该模型的正确性与实用性。
1 基本原理“近零法”和“指零法”是地磁测量主要采用的2种方法,二者在基本原理上没有区别,均采用了经纬仪的“四位置测量”方法,仅在操作和计算方法上略有不同。文中以近零法的测量为例进行分析。
仪器差测量基本流程如图 1所示。首先得出被测仪器及标准仪器的D的基线值DBX、DBS,二者之差为EDX,其中角标X表示被测仪器,S表示标准仪器、B表示基线值,D、I表示磁偏角和磁倾角,即被测仪器与标准仪器测量磁偏角D的差值,也称D的仪器差。即
$ D_{\mathrm{BX}}-D_{\mathrm{BS}}=E_{D \mathrm{X}} $ | (1) |
同理
$ I_{\mathrm{BX}}-I_{\mathrm{BS}}=E_{I \mathrm{X}} $ | (2) |
EIX即为磁倾角I的仪器差。
需要说明的是,仪器差本质是测量误差,根据计量规范《JJF1001—2011通用计量术语及定义》对“测量误差”的定义,EDX、EIX的正负号应与传统定义相反,本文采用该规范写法。
2 系统因素分析根据前文分析,现有的比测方法并未考虑测量过程中各干扰因素的影响,因此测量结果会有较大不确定度。以下对其中已知的、可分析的不确定度来源,即系统因素进行分析。由未知因素引起的、或不受控的随机因素引起的不确定度,一般可通过重复测量或者剔除异常值的方法消除,不做分析。
2.1 相对记录仪日变影响在理想状态下,基线值为常数,但由于磁通门磁力仪的三轴不能达到完全相互垂直,且在仪器架设时存在定向误差,导致记录的地磁三分量(D、H、Z)变化曲线会产生畸变,因而基线值也发生变化,这种变化称为日变,其在D和I上的变化量即为δRD、δRI,单位为“′”。图 2展示了乾陵地磁台GM4-XL型磁力仪与FGM01型磁力仪日间10 h的基线值,也就是日变记录。
2台仪器都呈现出变化的曲线,因此无论选择哪一台仪器,只要选用不同时间点的基线值进行仪器差计算,就可能引入日变带来的误差,所以,在去除仪器日变引起的误差后,式(1)、式(2)变为
$ D_{\mathrm{BX}}-D_{\mathrm{BS}}-\delta R_D=E_{D \mathrm{X}} $ | (3) |
$ I_{\mathrm{BX}}-I_{\mathrm{BS}}-\delta R_I=E_{I \mathrm{X}} $ | (4) |
以近零法测量磁偏角D为例,测得的偏角D修正值
$ \alpha_{D \mathrm{X}}=\frac{\sum \Delta S_{\mathrm{X}}}{4} \cdot \frac{T \cdot C_{\mathrm{X}}}{H} $ | (5) |
式中,ΣΔSX/4是4个位置检测计示数在赋予正负号后的值的平均,CX是被测仪器的修正系数,T = 3 438′/rad,是弧度转换为角度分的参数,H是测量点位磁场水平分量概值,为常量,单位nT。实际上,由于4个位置测量时有一定时间差,由测量时环境温度变化导致的基线值失稳情况时有发生。
图 3为榆林地磁台观测室2020年的温度变化与仪器零点偏移变化的关系,二者具有明显的相关性。因此在测量时,若温度变化过大,则一定会对零点偏移量产生影响,继而影响到测量结果。根据实际,设测量前后温度为线性单调变化,4个位置温度变化引起的偏移量被算入了修正值,该偏移量可按如下方法计算
$ \left.\begin{array}{l} \Delta \tau_{\mathrm{X} 1}=\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{4} \\ \Delta \tau_{\mathrm{X} 2}=-\frac{2 \Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{4} \\ \Delta \tau_{\mathrm{X} 3}=\frac{3 \Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{4} \\ \Delta \tau_{\mathrm{X} 4}=-\frac{4 \Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{4} \end{array}\right\} $ | (6) |
δθX为仪器传感器及电路部分的综合温度系数,是一个需要预估的变量。ΔtX是测量偏角D前后的环境温度差。则偏移量均值为
$ \frac{\sum \Delta \tau_{\mathrm{X}}}{4}=-\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8} $ | (7) |
此时需要引入一个理论上的真实值ΣΔSXT/4,则
$ \frac{\sum \Delta S_{\mathrm{X}}}{4}=\frac{\sum \Delta S_{\mathrm{XT}}}{4}+\frac{\sum \Delta \tau_{\mathrm{X}}}{4}=\frac{\sum \Delta S_{\mathrm{XT}}}{4}-\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8} $ | (8) |
对于修正系数CX,其通过每年1次的仪器标定得出,每次标定结果的误差以标准偏差来表示。通常情况下,比测与标定相隔数月,修正系数已发生了变化,这种变化对于日常的磁场测量,影响可忽略不计,但对于更高指标要求的仪器比测,则必须加以考虑,因此,若采用上次标定的修正系数进行计算,则
$ C_{\mathrm{X}}=C_{\mathrm{T}}+\delta C_{\mathrm{X}} $ | (9) |
式中,CT是引入的测量时的真实修正系数,δCX为误差。图 4为乾陵地磁台Mingeo型磁通门经纬仪修正系数自2015年以来的变化。
综上,设重复测量N次取平均值,则式(5)替换为
$ \alpha_{D \mathrm{X}}=\left(\frac{\sum \Delta S_{\mathrm{XT}}}{4 N}-\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N}\right) \cdot \frac{\left(C_{\mathrm{T}}+\delta C_{\mathrm{X}}\right) \cdot T}{H} $ | (10) |
式(5)中αDX为实测的修正值,为了分离出真实修正值部分,结合式(8)、式(9),得
$ \alpha_{D \mathrm{X}}=\frac{\sum \Delta S_{\mathrm{XT}}}{4 N} \cdot \frac{C_{\mathrm{T}} \cdot T}{H}-\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N} \cdot \frac{C_{\mathrm{X}} \cdot T}{H}+\left(\frac{\sum \Delta S_{\mathrm{X}}}{4 N}+\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N}\right) \cdot \frac{\delta C_{\mathrm{X}} \cdot T}{H} $ | (11) |
其中
$ \overline{\varepsilon_{D \mathrm{X}}}=\left(\frac{\sum \Delta S_{\mathrm{X}}}{4 N}+\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N}\right) \cdot \frac{\delta C_{\mathrm{X}} \cdot T}{H}-\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N} \cdot \frac{C_{\mathrm{X}} \cdot T}{H} $ | (12) |
同理,计算可得标准仪器及倾角I的修正值误差。
3 建立新的比测模型N次测量所得被测仪基线值均值
$ \begin{aligned} E_{D \mathrm{X}}= & \left(\overline{D_{\mathrm{BX}}}-\overline{\varepsilon_{D \mathrm{X}}}\right)-\left(\overline{D_{\mathrm{BS}}}-\overline{\varepsilon_{D \mathrm{S}}}\right)-\delta R_D \\ = & \overline{D_{\mathrm{BX}}}-\overline{D_{\mathrm{BS}}}-\delta R_D+\left[\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N} C_{\mathrm{X}}-\left(\frac{\sum \Delta S_{\mathrm{X}}}{4 N}+\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N}\right) \delta C_{\mathrm{X}}\right] \cdot \frac{T}{H}- \\ & {\left[\frac{\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N} C_{\mathrm{S}}-\left(\frac{\sum \Delta S_{\mathrm{S}}}{4 N}+\frac{\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N}\right) \delta C_{\mathrm{S}}\right] \cdot \frac{T}{H} } \end{aligned} $ | (13) |
同样得到磁倾角I的仪器差(近零法)为
$ \begin{aligned} E_{\mathrm{IX}}= & \left(\overline{I_{\mathrm{BX}}}-\overline{\varepsilon_{\mathrm{IX}}}\right)-\left(\overline{I_{\mathrm{BS}}}-\overline{\varepsilon_{\mathrm{IS}}}\right)-\delta R_{\mathrm{I}} \\ = & \overline{I_{\mathrm{BX}}}-\overline{I_{\mathrm{BS}}}-\delta R_I+\left[\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N} C_{\mathrm{X}}-\left(\frac{\sum \tilde{\Delta} S_{\mathrm{X}}}{4 N}+\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}}{8 N}\right) \delta C_{\mathrm{X}}\right] \cdot \frac{T}{F}- \\ & {\left[\frac{\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N} C_{\mathrm{S}}-\left(\frac{\sum \tilde{\Delta} S_{\mathrm{S}}}{4 N}+\frac{\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N}\right) \delta C_{\mathrm{S}}\right] \cdot \frac{T}{F} } \end{aligned} $ | (14) |
其中
$ \begin{aligned} E_{D \mathrm{X}} & =\left(\overline{D_{\mathrm{BX}}}-\overline{\varepsilon_{D \mathrm{X}}}\right)-\left(\overline{D_{\mathrm{BS}}}-\overline{\varepsilon_{D \mathrm{S}}}\right)-\delta R_D \\ & =\overline{D_{\mathrm{BX}}}-\overline{D_{\mathrm{BS}}}-\delta R_D+\left(\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}-\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N}\right) \cdot \frac{T}{H} \end{aligned} $ | (15) |
$ \begin{aligned} E_{I \mathrm{X}} & =\left(\overline{I_{\mathrm{BX}}}-\overline{\varepsilon_{I \mathrm{X}}}\right)-\left(\overline{I_{\mathrm{BS}}}-\overline{\varepsilon_{I \mathrm{S}}}\right)-\delta R_I \\ & =\overline{I_{\mathrm{BX}}}-\overline{I_{\mathrm{BS}}}-\delta R_I+\left(\frac{\Delta t_{\mathrm{X}} \cdot \delta \theta_{\mathrm{X}}-\Delta t_{\mathrm{S}} \cdot \delta \theta_{\mathrm{S}}}{8 N}\right) \cdot \frac{T}{F} \end{aligned} $ | (16) |
依托新的比测模型,采用GUM法对仪器差EDX、EIX进行测量不确定度的评定。由于式(13)、式(14)的比测模型为非线性,根据《JJF 1059.1—2012测量不确定度评定与表示》(中华人民共和国国家质量监督检验检疫总局,2013),合成标准不确定度应采用下式
$ u_{\mathrm{c}}^2(y)=\sum\limits_{i=1}^n\left(\frac{\partial f}{\partial x_i}\right)^2 u^2\left(x_i\right)+\sum\limits_{i=1}^n \sum\limits_{j=1}^n\left[\frac{1}{2}\left(\frac{\partial^2 f}{\partial x_i \partial x_j}\right)^2+\frac{\partial f}{\partial x_i} \cdot \frac{\partial^3 f}{\partial x_i \partial x_j^2}\right] u^2\left(x_i\right) u^2\left(x_j\right) $ | (17) |
继而计算出各输入量的标准不确定度u(xi)、灵敏系数ci以及合成不确定度uc。
依照上述方法,以某地磁台的一次磁通门经纬仪比测数据进行计算(近零法,量值溯源至国内基准),各变量估计范围及测量原始参数见表 1。
对各输入量的分布进行估计,继而得到各不确定度分量,见表 2、表 3。
EDX、EIX的标准不确定度uc(EDX) = 0.010 3′,uc(EDX) = 0.006 81′。EDX、EIX接近于正态分布,取k = 2,则EDX、EIX的扩展不确定度
$ U\left(E_{D \mathrm{X}}\right)=k u_{\mathrm{c}}\left(E_{D \mathrm{X}}\right)=0.021^{\prime} $ | (18) |
$ U\left(E_{I \mathrm{X}}\right)=k u_{\mathrm{c}}\left(E_{I \mathrm{X}}\right)=0.014^{\prime} $ | (19) |
最终,该仪器的仪器差
$ E_{D \mathrm{X}}=(0.030 \pm 0.021)^{\prime} \quad E_{I \mathrm{X}}=(-0.020 \pm 0.014)^{\prime} $ | (20) |
依据该结果,对被测仪器的仪器差进行合格性判定,MPE = 0.1′。
如图 5所示,2项结果都合格,显然这种表示和判定方法更可信,也与该仪器多年的比测结果相符合。
参考新的模型,使用2019年全国基准地磁台网比测的部分仪器测试数据(非最终结果,部分仪器进行了复测)进行不确定度评定并重新判定结果,所列结果多数在最大允许误差(MPE = 0.1′)附近,误判的可能性较大,结果如表 4。
按新的方法对12项结果判定后,有8项改判为“待定”,改判率达到67%,说明原判定结果存疑,需要谨慎对待,同时评定过程明确了问题所在,可针对性的查找原因,将系统因素的影响降到最低。例如乌加河的EDX = -0.04′,虽数值较小,但由于系统原因,其测量不确定度达到了0.069′,结果由合格变为待定。在新方法的评定中,
本研究建立了包含多种系统因素的、基于近零法和指零法的磁通门经纬仪比测模型,依托该模型,采用GUM方法,通过实际数据进行了测量不确定度评定实验,相比于传统比测模型和判定方法,结果表明:
(1)新的模型更加符合实际,依据该模型进行测量不确定度评定后,改进了结果的表示形式,优化了判定结果,该模型可以作为传统模型的替代。
(2)根据新模型的公式及不确定度的评定过程,得出了抑制测量过程中系统影响因素的方法:①测量次数N应控制在最佳范围内,N过小导致随机误差增大;若N过大,虽然随机误差减小,但延长了测量时间,导致相对记录仪日变、人员精力下降以及其它环境因素变化引起的不确定度上升。②因日变δR的不确定度对结果影响较大,所以应当科学合理地架设相对仪器,使其达到Z轴最佳的垂直度和最佳定向状态,并根据仪器的日变曲线选择日变化最小的时段进行测量工作。③控制标准仪器修正系数的变化量δC的范围,需要对仪器进行长期的观察,可以通过绘制长期的修正系数变化图来把握其变化规律,如果每年能增加标定次数,则更细化的变化规律有利于缩小变化量范围。④对于标准仪器,温度系数取值范围的确定需要测量员对仪器有所了解,或通过试验测得。对于被测仪器,无法进行试验,温度系数取值范围较难把握,但根据模型计算,温度系数引入的相关不确定度的灵敏系数c会随着ΔtX减小而以2倍数量级的速度减小,因此只要对观测室温差加以控制,温度系数引入的不确定度分量就将明显减小,甚至可忽略。
车濛琪, 应允翔, 汪继林, 等. 地磁观测质量与温度的关系分析[J]. 科技资讯, 2020, 18(33): 50-52. |
陈贤, 成万里, 黄恩贤. 信阳台地磁观测数据干扰分析及观测质量改善[J]. 高原地震, 2021, 33(3): 41-46. |
迟铖, 王丹, 吕俊伟, 等. 基于粒子群遗传算法的三轴磁通门误差校正[J]. 探测与控制学报, 2021, 43(3): 98-102. |
刘浩, 李予国, 丁学振, 等. 三轴磁通门传感器温漂误差分析及其校正方法[J]. 中国海洋大学学报, 2022, 52(5): 107-113. |
罗玉芬, 王建格, 陆镜辉. 肇庆地磁台基线值观测精度分析[J]. 华南地震, 2021, 41(3): 148-154. |
屈文斌. 基于单片机的磁通门传感器非正交性校正研究[J]. 仪表与自动化装置, 2020, 35(4): 62-65. |
宋思璇, 邓明, 陈凯, 等. 数字正交基模磁通门传感器电路[J]. 电子测量与仪器学报, 2020, 34(5): 97-102. |
王晓美, 滕云田, 谭婧, 等. 三轴磁通门传感器水平和定向对地磁日变化观测数据的影响[J]. 地震学报, 2017, 39(3): 429-435. |
姚远, 王晓美, 张平, 等. MINGEO-DI仪磁通门探头更换标定[J]. 地震地磁观测与研究, 2016, 37(5): 126-130. |
张策, 滕云田, 张涛, 等. 自动磁通门经纬仪多参量误差补偿算法[J]. 仪器仪表学报, 2020, 41(6): 85-93. |
张涛, 张策, 滕云田, 等. 地磁偏角倾角绝对测量技术发展现状综述[J]. 仪器仪表学报, 2018, 39(8): 80-90. |
中华人民共和国国家质量监督检验检疫总局. JJF 1059.1—2012测量不确定度评定与表示[S]. 北京: 中国标准出版社, 2013.
|
Jankowski J, Sucksdoroff C. Guide for magnetic measurements and observatory practice[S]. Boulder: International Association of Geomagnetism and Aeronomy, 1996.
|
Marsal S, Torta J M. An evaluation of the uncertainty associated with the measurement of the geomagnetic field with a D/I fluxgate theodolite[J]. Measurement Science and Technology, 2007, 18(7): 2 143-2 156. |