CN108802805A - 六自由度宽频带一体化地震监测系统 - Google Patents
六自由度宽频带一体化地震监测系统 Download PDFInfo
- Publication number
- CN108802805A CN108802805A CN201810580041.0A CN201810580041A CN108802805A CN 108802805 A CN108802805 A CN 108802805A CN 201810580041 A CN201810580041 A CN 201810580041A CN 108802805 A CN108802805 A CN 108802805A
- Authority
- CN
- China
- Prior art keywords
- data
- parameter
- gnss
- acceleration
- angular speed
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000012544 monitoring process Methods 0.000 title claims abstract description 23
- 230000010354 integration Effects 0.000 title claims abstract description 19
- 230000001133 acceleration Effects 0.000 claims abstract description 95
- 230000004927 fusion Effects 0.000 claims abstract description 43
- 238000000034 method Methods 0.000 claims abstract description 31
- 230000004807 localization Effects 0.000 claims abstract description 15
- 238000005259 measurement Methods 0.000 claims abstract description 14
- 238000004364 calculation method Methods 0.000 claims abstract description 13
- 230000015572 biosynthetic process Effects 0.000 claims abstract description 10
- 238000003786 synthesis reaction Methods 0.000 claims abstract description 10
- 230000001360 synchronised effect Effects 0.000 claims abstract description 8
- 238000006073 displacement reaction Methods 0.000 claims description 47
- 238000012545 processing Methods 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 10
- 230000003068 static effect Effects 0.000 claims description 8
- 241001269238 Data Species 0.000 claims description 5
- 238000005267 amalgamation Methods 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 3
- 241000208340 Araliaceae Species 0.000 claims description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 235000008434 ginseng Nutrition 0.000 claims description 2
- 239000007787 solid Substances 0.000 claims 1
- 230000033001 locomotion Effects 0.000 description 7
- 230000009897 systematic effect Effects 0.000 description 7
- 238000001914 filtration Methods 0.000 description 6
- 230000007704 transition Effects 0.000 description 5
- 238000012937 correction Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 238000013519 translation Methods 0.000 description 3
- 230000005540 biological transmission Effects 0.000 description 2
- 230000000295 complement effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 235000008331 Pinus X rigitaeda Nutrition 0.000 description 1
- 235000011613 Pinus brutia Nutrition 0.000 description 1
- 241000018646 Pinus brutia Species 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000000632 dystonic effect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 230000001568 sexual effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/01—Measuring or predicting earthquakes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2200/00—Details of seismic or acoustic prospecting or detecting in general
- G01V2200/10—Miscellaneous details
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Gyroscopes (AREA)
- Navigation (AREA)
Abstract
本发明公开了六自由度宽频带一体化地震监测系统，至少包括传感器集成单元和数据融合单元；传感器集成单元包括同址安置的GNSS接收机、三轴加速度计和三轴陀螺仪；GNSS接收机、三轴加速度计、三轴陀螺仪用来同步采集GNSS观测数据、加速度观测数据、角速度观测数据并存储；数据融合单元包括GNSS数据定位处理模块、集成融合测量数据模块和综合滤波器；GNSS数据定位处理模块用来根据GNSS观测数据进行定位解算；集成融合测量数据模块用来将GNSS位置信息、加速度观测数据和角速度观测数据进行融合；综合滤波器用来根据融合后的数据进行解算。本发明能实现单站六自由度地震监测，可拓宽地震监测带宽，并提升地震监测能力。
Description
技术领域
本发明涉及一种六自由度宽频带一体化地震监测系统。
背景技术
完整准确地获取和恢复近场强震事件在各时空频段上的地震波信号是深入研究和理解大地震的重要基础。然而，当前近场强震研究的一个重要缺憾是极低频观测信息(如20秒以上)可能淹没于各类噪声或系统误差之中，不易提取可靠的地震波信号，成为理解大地震震源物理机制的关键限制性因素之一。
近场强震动监测中，理论上加速度计和高频GNSS都能够捕捉到大地震在低频段上所释放的能量，但实践中它们在观测质量和硬件性能上存在不足，无法完整获取低频段地震波信号。从仪器的观测质量上看，加速度计的噪声可低至0.6μm/s2/√Hz，也具备测量“零频率”信号的能力，充分满足地震观测的精度和带宽需求，但强烈的地震波会引发加速度计的显著倾斜，使得重力加速度投影到平面方向的观测值中，积分时便产生分米级甚至米级的基线漂移误差，阻碍了极低频率上的地震平移观测量的准确恢复。因此，在加速度计数据处理中通常被迫实施高通滤波以消除低频段的系统误差，同时也丧失了包括永久形变量在内的重要地震信号。另一方面，以中国北斗卫星系统为代表所支撑的高频GNSS对地观测技术直接测量动态位移，没有基线漂移误差，其平面精度一般可达2cm～4cm，垂直精度约5cm～8cm(即2倍中误差)，是一种监测近场强震动平移量的重要手段。
并置GNSS接收机和加速度计这两种仪器，并融合其输出位移和加速度获得可靠的宽频带(包含直流DC到20秒频段)地震平移量观测，既保留GNSS获取的较为可信的低频位移信号，也保持加速度计所捕获的相对高精度的高频信号。
这种通过现成仪器的非联通过网络传输或就地解算。接性并置，并在数据处理层面松组合高频GNSS和加速度计平移观测量的技术存在以下不足：
(1)加速度计观测所受到的地震旋转量污染并没有得到补偿，以致与GNSS融合时必须降权，使得输出位移仍然受到GNSS低频系统误差的干扰；
(2)GNSS和加速度计的数据融合策略及研究较多，但都不能推广到陀螺仪数据。现有的各种GNSS、加速度计融合策略普遍将地震旋转量造成的加速度计基线漂移作为一种无法补偿的误差，因此期冀从多源数据权重比的角度来解决这个问题，导致这些方法不能兼容地震旋转量观测数据，限制了与陀螺仪观测数据的进一步耦合。惯性传感器陀螺仪观测量的低频部分也会受到各种系统误差的影响。因此，积分陀螺仪观测量时我们也需要首先实施高通滤波，但同时失去了永久旋转形变量等重要的地震信号；
理论上，加速度计和陀螺仪是互补的传感器。如校正陀螺仪的线加速度敏感性离不开加速度计，而后者的旋转量修正又需要前者的角度观测量，因此将它们集成为一种六自由度监测仪可以校正其各自的系统误差。然而，这二者的系统误差相互依赖于对方准确的平移和旋转观测量，原理实现上其实是自我矛盾的，该问题的根源在于这种经典的六自由度仪器缺失了一个作为参考基准的非惯性检核传感器。在近场强震动情况下，仅仅依靠传统惯性传感器所形成的六自由度监测平台难以克服仪器自身的系统误差，在极低频段上不具有可靠性。GNSS可以直接获取DC的形变信息，与加速度计、陀螺仪协同，能够互相补充。
发明内容
本发明的目的是提供一种组合了GNSS接收机、加速度计和陀螺仪的六自由度宽频带一体化地震监测系统。
本发明能够在同址直接获取GNSS接收机、加速度计和陀螺仪的观测数据，并在数据处理层面将三种观测数据融合，最终输出地震宽频带、六自由度完备观测量。本发明通过多源观测数据融合，来弥补各传感器在观测噪声和系统误差上的相对劣势，以获取单一传感器所不易或不能获得的地震波信息，尤其是极低频段上的地震波信息,从而获得精确可靠的六自由度、宽频带地震位移和旋转量。
本发明提供的六自由度宽频带一体化地震监测系统，至少包括传感器集成单元和数据融合单元；
所述传感器集成单元进一步包括同址安置的GNSS接收机、三轴加速度计和三轴陀螺仪；GNSS接收机、三轴加速度计、三轴陀螺仪用来同步采集GNSS观测数据、加速度观测数据、角速度观测数据并存储；
所述数据融合单元进一步包括GNSS数据定位处理模块、集成融合测量数据模块和综合滤波器；其中，GNSS数据定位处理模块用来根据GNSS观测数据进行定位解算，获得地心地固坐标系下的GNSS位置信息；集成融合测量数据模块用来将GNSS位置信息、加速度观测数据和角速度观测数据进行融合；综合滤波器用来根据融合后的数据进行解算，获得融合解算后的待求解参数，所述待求解参数包括位置信息、速度信息、加速度信息、姿态角信息和角速度信息。
进一步的，本发明系统还包括信息输出单元，用来输出传感器集成单元同步采集的原始观测数据、以及根据融合后的位置信息、速度信息、加速度信息、姿态角信息和角速度信息获得的含零频位移时间序列、速度时间序列和旋转性形变时间序列。
进一步的，GNSS数据定位处理模块利用相对定位法或精密单点定位法进行定位解算。
进一步的，集成融合测量数据模块用来将GNSS位置信息、加速度观测数据和角速度观测数据进行融合，进一步包括：
用来将地心地固坐标系下的GNSS位置信息转换至测站“东北高”坐标系下；
用来对加速度观测数据进行零偏和线性趋势估计；
用来对角速度观测数据进行零偏估计；
用来根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，计算GNSS接收机、加速度计和陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在测站“东北高”坐标系三轴方向的权重。
进一步的，综合滤波器用来根据融合后的数据进行解算，进一步包括：
用来根据融合后的数据，结合位置、加速度、角速度观测方程、待求解参数间的递推关系以及GNSS位置信息、加速度观测数据、角速度观测数据的权重，解算出待求解参数；
所述位置观测方程的构建为：利用测站“东北高”坐标系下的GNSS位置信息和GNSS位置误差，来表示待求解位置参数；
所述加速度观测方程的构建为：利用加速度观测数据和加速度计误差，来表示待求解加速度参数和待求解姿态角参数；
所述角速度观测方程的构建为：利用角速度观测数据和陀螺仪误差，来表示待求解角速度参数；
所述待求解参数间的递推关系的构建为：利用第i历元的待求解速度参数和待求解加速度参数表示第i+1历元的待求解速度参数，构建待求解速度参数的递推关系；利用第i历元的待求解位置参数、待求解速度参数和待求解加速度参数表示第i+1历元的待求解位置参数，构建待求解位置参数的递推关系；利用第i历元的待求解姿态角参数和待求解角速度参数表示第i+1历元的待求解姿态角参数，构建待求解姿态角参数的递推关系；
所述权重的获得为：根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，计算GNSS接收机、加速度计和陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在测站“东北高”坐标系三轴方向的权重。
本发明提供的GNSS接收机、加速度计和陀螺仪观测数据的融合处理方法，包括步骤：
S100：对原始观测数据进行预处理，所述原始观测数据包括GNSS接收机、加速度计和陀螺仪同步采集的GNSS观测数据、加速度观测数据和角速度观测数据；
S200：根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，对GNSS位置信息、加速度观测数据、角速度观测数据进行定权；
S300：构建位置、加速度、角速度观测方程；本步骤进一步包括：
S310：利用预处理后的GNSS位置信息和GNSS位置误差，构建位置观测方程其中，[xe xn xu]T表示预处理后的GNSS位置信息，[ξe ξn ξu]T表示 GNSS位置误差，[pe pn pu]T表示待求解的位置参数；
S320：利用预处理后的加速度观测数据和加速度计误差，构建加速度观测方程其中，[Ae'An'Au']T表示预处理后的加速度观测数据，[εe εn εu]T表示加速度计误差，[ae an au]T表示待求解的加速度参数，[φ θ ψ]T表示待求解的姿态角参数，g表示重力加速度；
S330：利用预处理后的角速度观测数据和陀螺仪误差，构建角速度观测方程其中，[re rn ru]T表示预处理后的角速度观测数据，[ζe ζn ζu]T表示陀螺仪误差，[ge gn gu]T表示待求解的角速度参数；
S400：构建相邻历元的待求解参数间的递推关系，所述待求解参数包括待求解速度参数、待求解加速度参数、待求解位置参数、待求解姿态角参数和待求解角速度参数；
所构建的待求解参数间的递推关系包括：
其中，pi和pi+1分别表示第i历元和第i+1历元下的待求解位置参数；vi和vi+1分别表示第i历元和第i+1历元下的待求解速度参数；τ表示采样间隔；ai表示第i历元下的待求解加速度参数；[φi θi ψi]T和[φi+1 θi+1 ψi+1]T分别表示第i历元和第i+1历元下的待求解姿态角参数；[ge gn gu]T表示待求解角速度参数；式(3)下标中i和i+1分别表示第 i历元和第i+1历元对应的待求解参数；表示3行3列的单位矩阵；表示3行3列的零矩阵；v＝[ve vnvu]T，表示待求解速度参数；
S500：结合位置、加速度、角速度观测方程、待求解参数间的递推关系以及GNSS 位置信息、加速度观测数据、角速度观测数据的权重，解算出待求解参数。
进一步的，所述预处理进一步包括：
S110：根据GNSS观测数据解算地心地固坐标系下的GNSS位置信息，并将GNSS 位置信息转换至测站“东北高”坐标系下；
S120：对加速度观测数据进行零偏和线性趋势估计；
S130：对角速度观测数据进行零偏估计。
所述子步骤S110，具体为：
S111：地心地固坐标系下的GNSS位置信息转换到大地坐标系下；
S112：利用事后精密单点定位方法计算的单天解坐标计算大地坐标系到测站“东北高”坐标系的转换矩阵；
S113：利用转换矩阵将大地坐标系下的GNSS位置信息转换到至测站“东北高”坐标系下。
进一步的，步骤S200，具体为：
计算GNSS接收机、加速度计和陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在测站“东北高”坐标系三轴方向的权重。
本发明的有益效果如下：
(1)汲取了GNSS和强震仪各自的优点，避免了强震仪容易受到倾斜和旋转导致的基线漂移的影响，能够在静态、低频震动和强震动多种条件、低频到高频多种频带上，获取到高精度的三方向位移和速度信息。
(2)组合了陀螺仪，能够直接获取到地震动的角速度信息，可以直接对强震仪的加速度数据进行倾斜、旋转改正，避免了像传统数据融合方式那样对加速度数据盲目降权。同时，三轴旋转观测量，结合三向加速度观测量，以及GNSS观测量，能够获取六自由度的地震监测结果，极大地扩充了地震的可监测范围，使得地震单站六自由度监测成为可能，为地震震源物理机制等地球科学研究提供了数据观测基础。
(3)实现了三种测量仪器数据融合，能够实现单站六自由度地震监测，获取准确的宽频带的地震平移量和地震旋转量，极大拓宽地震监测带宽，提升了地震的监测能力。
附图说明
图1是本发明系统具体的构架图；
图2是本发明系统具体的数据解算流程图；
图3是实施例中本发明系统与加速度计阵列处理得到的旋转角对比图；
图4是实施例中GPS加速度计两种仪器融合与GPS单独获得的位移对比图；
图5是实施例中本发明系统融合位移与GPS单独获得的位移对比图
图6是实施例中振动过程中振动台记录的振动方向的真实位移图；
图7是实施例中加速度计高通滤波并积分得到的位移；
图8是实施例中GPS加速度计两种仪器融合位移与真值差异以及本发明系统位移与真值差异的对比图；
图9是实施例中加速度计积分位移和GPS加速度计两种仪器融合与振动台真值的差异对比图。
具体实施方式
为了更清楚地说明本发明和/或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。
应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
如图1所示，本发明一种具体的系统架构包括传感器集成单元、数据融合单元和信息输出单元。传感器集成单元包括GNSS接收机、加速度计和陀螺仪；数据融合单元包括GNSS数据定位处理模块、集成融合测量数据模块和综合滤波器；信息输出单元用来输出数据，所输出数据包括但不限于传感器集成单元输出的原始观测数据、以及综合滤波器输出的含零频位移时间序列、速度时间序列和旋转性形变时间序列。
需要注意的是，本发明中，加速度计指三轴加速度计，陀螺仪指三轴陀螺仪。
本发明系统的具体实施过程如下：
(1)系统安装与初始化。
对GNSS接收机、加速度计和陀螺仪进行同址安置，使得GNSS接收机、加速度计和陀螺仪位于同一坐标基准下，且时间同步。
首先，正确确定地理北方向，即地心地固坐标系的北方向。对于GNSS接收机及其天线，应将天线通过刚性连接杆固定在观测墩上，保持天线水平并将天线按照规范指向特定方向，并正确量取天线高。对于加速度计和陀螺仪两种惯性仪器，应将二者安置于水平台面上，并将仪器坐标系北方向对准地理北方向，保证加速度计和陀螺仪自身的惯性坐标系和测站“东北高”坐标系对齐，从而保证所有的数据处理在同一坐标基准下进行。
本发明中，仪器坐标系指仪器自身所带的测量坐标系，调整加速度计和陀螺仪指北轴指向地理北方向，即实现仪器坐标系北方向对准地理北方向。大地测量学中，测站“东北高”坐标系即站心坐标系，站心坐标系中X轴指东，Y轴指北，Z轴指天顶方向。本发明中，站心坐标系记为ENU坐标系。
三种传感器仪器的设置与初始化：
设置GNSS接收机采样率；可选1Hz、5Hz、10Hz、20Hz或更高；设置加速度计和陀螺仪的采样率，可选100Hz或更高。其中，加速度计和陀螺仪在安置前应当调试妥当。将GNSS授时天线连接好，确定系统各部分工作正常。本实施例中，GNSS接收机的采样率设为1Hz，加速度计的采样率设为100Hz，陀螺仪的采样率设为100Hz。
(2)数据采集、传输与存储。
GNSS接收机用来接收GNSS观测数据，所述GNSS观测数据包括GNSS伪距观测值和/或GNSS载波相位观测值。根据GNSS观测数据可进行相对定位或精密单点定位数据解算，可获得测站位置。当进行精密单点定位数据解算时，需利用通过外部网络发送过来的精密卫星星历和钟差文件。加速度计用来获取加速度观测数据，陀螺仪用来获取角速度观测数据。
GNSS接收机接收的GNSS观测数据、加速度计获得的加速度观测数据、以及陀螺仪获得的角速度观测数据统记为原始数据。对原始数据进行同步采集和存储，并传输到数据融合单元。
(3)数据融合处理方法：
数据融合进一步包括GNSS数据定位处理和集成融合测量数据。GNSS数据定位处理可采用相对定位法或精密单点定位法，前者需要基准站或改正数作为参考，后者需要精密卫星星历和钟差。解算后，获得测站在地心地固坐标系下精确的绝对坐标。由于同震位移监测一般采用测站的站心坐标系，即需要利用坐标系旋转矩阵将绝对坐标转换到测站“东北高”坐标系，此处需要精确的测站已知坐标作为参考坐标，该参考坐标可以利用事后精密单点定位方法计算的单天解坐标。
下面结合实例和附图对本发明的集成融合测量处理作进一步详细的说明。
如图2所示的数据融合流程，具体步骤如下：
步骤A：获取原始数据。
通过系统获取包括GNSS位置信息、加速度观测数据以及角速度观测数据。GNSS 位置信息为利用相对定位法或精密单点定位法解算出的地心地固坐标系下测站的绝对坐标。加速度计数据是标有GNSS时标的仪器坐标系下三个轴方向的加速度观测数据，单位为gal，即cm/s2。未经任何处理的原始加速度观测数据在三个轴向上分别包含零偏和漂移，需要在预处理的阶段估计。角速度观测数据是标有GNSS时标的仪器坐标系下绕三个轴方向的角速度数据，单位为rad/s。未经任何处理的原始角速度观测数据在三个轴向上包含零偏，也需要在预处理阶段估计。
步骤B：原始数据的预处理。
由于加速度计获取的是仪器自身惯性坐标系下的加速度，陀螺仪获取的是自身惯性坐标系下的三轴角速度，陀螺仪和加速度计的观测噪声，是后续数据融合的必要参数，因此需要进行记录和计算。由于自身和环境噪声影响，加速度计原始的加速度观测数据存在零漂等多种误差，因此，使用前需要对其进行预处理，即对零漂、线性趋势等进行去除或者估计。去除零漂可以将数据统一扣除静止时一段时间的观测平均值，去除线性趋势道理也类似，即先拟合再扣除的方式。陀螺仪原始的角速度观测数据本身也是具有多种观测误差的，因而，使用前也应进行去除零点漂移的处理，并且要根据一段时间的观测值，估计陀螺仪本身的噪声。
下面将逐一描述原始数据的预处理过程，包括绝对坐标的转换、以及零偏和线性趋势的估计。
第1步：将GNSS位置信息转换至测站“东北高”坐标系的ENU方向。
首先，将地心地固坐标系下绝对坐标转换成大地坐标系下的“大地纬度、大地经度、大地高”坐标，即BLH坐标；之后，计算BLH坐标到ENU坐标的转换矩阵rot，转换矩阵rot见式(1)：
式(1)中，B、L表示大地坐标系下的大地纬度、大地经度。
利用转换矩阵rot将BLH坐标转换至ENU坐标。
第2步：对加速度观测数据的零偏和线性趋势进行估计。
由于加速度计输出的原始数据(即加速度观测数据)存在零偏和线性趋势，所以需要预先估计加速度计原始数据中的零偏和线性趋势，具体可采用式(2)估计：
式(2)中，Ae、An、Au表示加速度计输出的三个轴向的原始数据；t表示时间；kten1、kten2、kten3表示加速度计三个轴向上的线性趋势；bias1、bias2、bias3表示加速度计三个轴向上的零偏。
选取震动开始前10分钟静止状态下的加速度观测数据，对该静止状态下的加速度观测数据Ae、An、Au取平均和进行线性拟合，即可估计kten1、kten2、kten3、bias1、bias2、bias3。
第3步：对角速度观测数据的零偏进行估计。
与加速度计一样，陀螺仪输出的原始数据(即角速度观测数据)中也存在零偏，同样需要预先估计角速度观测数据中的零偏。估计零偏的方法为：选取震动开始前10分钟静止状态下的角速度观测数据，对该静止状态下的角速度观测数据取平均，即可估计陀螺仪的零偏。
步骤C：对参数进行定权。
本发明监测系统包含15个参数，分别为3个位置参数、3个速度参数、3个加速度参数、3个姿态角参数，3个角速度参数。定权之前，用震动开始前10分钟静止状态下的GNSS位置、加速度和角速度数据，分别计算GNSS接收机、加速度计、陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在ENU三轴方向的权重。将GNSS位置信息在ENU三轴方向的权重记为位置权重，将加速度观测数据在ENU三轴方向的权重记为加速度权重，将角速度观测数据在ENU三轴方向的权重记为角速度权重。
步骤D：将预处理后的数据加入法方程，构建位置、加速度、角速度的观测方程。
第1步：建立位置观测方程。
待求解位置参数与位置观测值间的关系见图3，也即所构建的位置观测方程：
式(3)中，[xe xn xu]T表示测站“东北高”坐标系下的GNSS位置信息，[ξe ξn ξu]T表示GNSS位置误差，[pe pn pu]T表示待求解的位置参数。
第2步：建立加速度观测方程。
待求解加速度参数、待求解姿态角参数与加速度观测值间的关系见式(4)，也即所构建的加速度观测方程：
式(4)中，[Ae' An' Au']T表示预处理后的加速度观测数据，[εe εn εu]T表示加速度计误差，[ae an au]T表示待求解的加速度参数，[φ θ ψ]T表示待求解的姿态角参数，g表示重力加速度。
第3步：建立陀螺仪观测值角速度观测方程。
待求解角速度参数与陀螺仪角速度观测数据之间的关系见式(5)，也即所构建的角速度观测方程：
式(5)中，[re rn ru]T表示预处理后的角速度观测数据，[ζe ζn ζu]T表示陀螺仪误差， [ge gn gu]T表示待求解的角速度参数。
步骤E：建立待求解参数间的递推关系。
建立位置、速度、加速度、姿态角与角速度参数之间的关系，以将参数向前递推一个历元。
构建待求解位置参数、待求解速度参数、待求解加速度参数间的数学关系，如下：
式(6)中，pi和pi+1分别表示第i历元和第i+1历元下的待求解位置参数；vi和vi+1分别表示第i历元和第i+1历元下的待求解速度参数；τ表示采样间隔，即采样率的倒数，单位：秒；ai表示第i历元下的待求解加速度参数。
构建姿态角参数与角速度参数间的数学关系，如下：
式(7)中，[φi θi ψi]T和[φi+1 θi+1 ψi+1]T分别表示第i历元和第i+1历元下的待求解姿态角参数；[ge gn gu]T表示待求解角速度参数。
由式(6)～(7)构建待求解参数的递推关系，如下：
式(8)中，下标中i和i+1分别表示第i历元和第i+1历元对应的待求解参数。
步骤F：解算各参数。
联合方程(3)～(8)，采用最小二乘方法解算各待求解参数的最优值，并计算最优值和观测值之差，即残差。
在系统获得GNSS位置信息、加速度计数据以及陀螺仪测量得到的角速率数据之后，按照以上数据融合处理算法步骤对数据进行处理，可以输出得到融合之后的位置信息、速度信息、加速度信息以及仪器本身的旋转信息，以及原始数据经过解算过后的残差。
实施例
图3至图9是在2018年4月10日于武汉地区(114.2927°E，30.3056°N)，用振动台及GPS接收机、加速度计、陀螺仪模拟地震条件观测得到的试验数据。
如图3所示，为本发明系统处理得到的在测站“东北高(ENU)”坐标系下角度旋转量与4台加速度计阵列计算得到的3个轴向上的角度旋转量的对比，图中，X、Y、Z 轴分别对应ENU坐标系的三轴。图中，实线表示本发明系统计算得到的绕三轴的角度旋转量，虚线表示4台加速度计阵列计算得到的绕三轴的角度旋转量，可以认为是旋转角大小的真值。由图中对比可以看出，本发明系统计算得到的角度旋转量与加速度计阵列计算得到的角度很接近。
图4所示为GPS与加速度计两种仪器融合计算得到的位移与GPS单独得到的位移的对比，图中，实线表示GPS与加速度计两种仪器融合计算的E、N、U三个方向的位移曲线，虚线表示GPS单独得到的E、N、U三个方向的位移曲线。图5所示为本发明系统计算得到的位移与GPS单独得到的位移的对比，图中，实线表示本发明系统计算的E、N、U三个方向的位移曲线，虚线表示GPS单独得到的E、N、U三个方向的位移曲线。从直观上看，本发明系统得到的位移与GPS加速度计两者融合得到的位移没有明显差别。可以发现，利用本发明系统能够获得振动方向E方向的永久位移，约为 2cm，这是图7加速度计高通滤波积分无法实现的。但是振动台记录了振动方向较为精确的位移，所以，需要在振动方向上将GPS加速度计两者融合的位移与本发明系统计算得到的位移与振动台记录的位移作比较。
图6所示为振动台记录的较精确的振动方向的位移。图8所示为发明系统融合位移与振动台精确位移求差值后的对比结果，图8中，虚线为GPS和加速度计两者融合获得的位移相对振动台输出位移的差异，实线为本发明系统计算的位移相对于振动台输出位移的差异。从图中可以看出，虽然只有毫米级的差异，但是本发明系统的计算结果在整体上跟接近振动台的真实位移。
图7所示为加速度计结合高通滤波并积分获得的E、N、U三个方向的位移时间序列。同样将加速度计积分得到的位移与振动台输出位移在振动方向上作对比，图9为对比结果。图9中，虚线为加速度计结合高通滤波并积分后的位移与振动台输出位移的差异，实线为本发明系统计算的位移相对于振动台输出位移的差异。从图中可以明显的对比出，本发明系统的解算结果要更精确。
经过上述的数据融合处理，对于一次完整的地震事件，本发明系统可得到精确的测站位移时间序列、速度时间序列以及测站旋转角度时间序列。与此同时，还可以将高采样率的GNSS、加速度计、陀螺仪原始数据作为备份，可以用作后续的精细处理和研究。
以上给出了本发明涉及的具体实施方式，但本发明不局限于所描述的实施方式。在本发明给出的思路下，采用对本领域技术人员而言容易想到的方式对上述实施例中的技术手段进行变换、替换、修改，并且起到的作用与本发明中的相应技术手段基本相同、实现的发明目的也基本相同，这样形成的技术方案是对上述实施例进行微调形成的，这种技术方案仍落入本发明的保护范围内。
Claims (7)
1.六自由度宽频带一体化地震监测系统，其特征是：
至少包括传感器集成单元和数据融合单元；
所述传感器集成单元进一步包括同址安置的GNSS接收机、三轴加速度计和三轴陀螺仪；GNSS接收机、三轴加速度计、三轴陀螺仪用来同步采集GNSS观测数据、加速度观测数据、角速度观测数据并存储；
所述数据融合单元进一步包括GNSS数据定位处理模块、集成融合测量数据模块和综合滤波器；其中，GNSS数据定位处理模块用来根据GNSS观测数据进行定位解算，获得地心地固坐标系下的GNSS位置信息；集成融合测量数据模块用来将GNSS位置信息、加速度观测数据和角速度观测数据进行融合；综合滤波器用来根据融合后的数据进行解算，获得融合解算后的待求解参数，所述待求解参数包括位置信息、速度信息、加速度信息、姿态角信息和角速度信息。
2.如权利要求1所述的六自由度宽频带一体化地震监测系统，其特征是：
还包括信息输出单元，用来输出传感器集成单元同步采集的原始观测数据、以及根据融合后的位置信息、速度信息、加速度信息、姿态角信息和角速度信息获得的含零频位移时间序列、速度时间序列和旋转性形变时间序列。
3.如权利要求1所述的六自由度宽频带一体化地震监测系统，其特征是：
所述GNSS数据定位处理模块利用相对定位法或精密单点定位法进行定位解算。
4.如权利要求1所述的六自由度宽频带一体化地震监测系统，其特征是：
所述集成融合测量数据模块用来将GNSS位置信息、加速度观测数据和角速度观测数据进行融合，进一步包括：
用来将地心地固坐标系下的GNSS位置信息转换至测站“东北高”坐标系下；
用来对加速度观测数据进行零偏和线性趋势估计；
用来对角速度观测数据进行零偏估计；
用来根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，计算GNSS接收机、加速度计和陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在测站“东北高”坐标系三轴方向的权重。
5.如权利要求1所述的六自由度宽频带一体化地震监测系统，其特征是：
所述综合滤波器用来根据融合后的数据进行解算，进一步包括：
用来根据融合后的数据，结合位置、加速度、角速度观测方程、待求解参数间的递推关系以及GNSS位置信息、加速度观测数据、角速度观测数据的权重，解算出待求解参数；
所述位置观测方程的构建为：利用测站“东北高”坐标系下的GNSS位置信息和GNSS位置误差，来表示待求解位置参数；
所述加速度观测方程的构建为：利用加速度观测数据和加速度计误差，来表示待求解加速度参数和待求解姿态角参数；
所述角速度观测方程的构建为：利用角速度观测数据和陀螺仪误差，来表示待求解角速度参数；
所述待求解参数间的递推关系的构建为：利用第i历元的待求解速度参数和待求解加速度参数表示第i+1历元的待求解速度参数，构建待求解速度参数的递推关系；利用第i历元的待求解位置参数、待求解速度参数和待求解加速度参数表示第i+1历元的待求解位置参数，构建待求解位置参数的递推关系；利用第i历元的待求解姿态角参数和待求解角速度参数表示第i+1历元的待求解姿态角参数，构建待求解姿态角参数的递推关系；
所述权重的获得为：根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，计算GNSS接收机、加速度计和陀螺仪的三个轴向输出数据的标准差或方差，以标准差或方差的倒数作为GNSS位置信息、加速度观测数据、角速度观测数据在测站“东北高”坐标系三轴方向的权重。
6.GNSS接收机、加速度计和陀螺仪观测数据的融合处理方法，其特征是，包括步骤：
S100：对原始观测数据进行预处理，所述原始观测数据包括GNSS接收机、加速度计和陀螺仪同步采集的GNSS观测数据、加速度观测数据和角速度观测数据；
S200：根据震动开始前预设时长静止状态下预处理后的GNSS位置信息、加速度观测数据和角速度观测数据，对GNSS位置信息、加速度观测数据、角速度观测数据进行定权；
S300：构建位置、加速度、角速度观测方程；本步骤进一步包括：
S310：利用预处理后的GNSS位置信息和GNSS位置误差，构建位置观测方程其中，[xe xn xu]T表示预处理后的GNSS位置信息，[ξe ξn ξu]T表示GNSS位置误差，[pe pn pu]T表示待求解的位置参数；
S320：利用预处理后的加速度观测数据和加速度计误差，构建加速度观测方程其中，[Ae'An'Au']T表示预处理后的加速度观测数据，[εe εn εu]T表示加速度计误差，[ae an au]T表示待求解的加速度参数，[φθ ψ]T表示待求解的姿态角参数，g表示重力加速度；
S330：利用预处理后的角速度观测数据和陀螺仪误差，构建角速度观测方程其中，[re rn ru]T表示预处理后的角速度观测数据，[ζe ζn ζu]T表示陀螺仪误差，[ge gn gu]T表示待求解的角速度参数；
S400：构建相邻历元的待求解参数间的递推关系，所述待求解参数包括待求解速度参数、待求解加速度参数、待求解位置参数、待求解姿态角参数和待求解角速度参数；
所构建的待求解参数间的递推关系包括：
其中，pi和pi+1分别表示第i历元和第i+1历元下的待求解位置参数；vi和vi+1分别表示第i历元和第i+1历元下的待求解速度参数；τ表示采样间隔；ai表示第i历元下的待求解加速度参数；[φi θi ψi]T和[φi+1 θi+1 ψi+1]T分别表示第i历元和第i+1历元下的待求解姿态角参数；[ge gn gu]T表示待求解角速度参数；式(3)下标中i和i+1分别表示第i历元和第i+1历元对应的待求解参数；表示3行3列的单位矩阵；表示3行3列的零矩阵；v＝[ve vn vu]T，表示待求解速度参数；
S500：结合位置、加速度、角速度观测方程、待求解参数间的递推关系以及GNSS位置信息、加速度观测数据、角速度观测数据的权重，解算出待求解参数。
7.如权利要求6所述的GNSS接收机、加速度计和陀螺仪观测数据的融合处理方法，其特征是：
所述预处理进一步包括：
S110：根据GNSS观测数据解算地心地固坐标系下的GNSS位置信息，并将GNSS位置信息转换至测站“东北高”坐标系下；
S120：对加速度观测数据进行零偏和线性趋势估计；
S130：对角速度观测数据进行零偏估计。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810580041.0A CN108802805B (zh) | 2018-06-07 | 2018-06-07 | 六自由度宽频带一体化地震监测系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810580041.0A CN108802805B (zh) | 2018-06-07 | 2018-06-07 | 六自由度宽频带一体化地震监测系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108802805A true CN108802805A (zh) | 2018-11-13 |
| CN108802805B CN108802805B (zh) | 2019-10-11 |
Family
ID=64087523
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810580041.0A Active CN108802805B (zh) | 2018-06-07 | 2018-06-07 | 六自由度宽频带一体化地震监测系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108802805B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113126143A (zh) * | 2021-04-22 | 2021-07-16 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种便携式微动及强震动实时监测仪系统 |
| CN113805223A (zh) * | 2021-08-16 | 2021-12-17 | 南京天巡遥感技术研究院有限公司 | 一种地震勘探系统及其采集数据的处理方法 |
| CN114167493A (zh) * | 2021-11-23 | 2022-03-11 | 武汉大学 | Gnss双天线辅助陀螺的地震旋转测量系统及方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1916567A (zh) * | 2006-09-04 | 2007-02-21 | 南京航空航天大学 | 基于自适应闭环h∞滤波器的对北斗双星/捷联惯导组合导航系统进行修正的方法 |
| CN101726752A (zh) * | 2008-10-23 | 2010-06-09 | 鸿富锦精密工业（深圳）有限公司 | 地震监测系统 |
| CN102608625A (zh) * | 2012-03-30 | 2012-07-25 | 武汉大学 | 基于惯性辅助定位接收机的实时形变监测预警系统及方法 |
| CN103235328A (zh) * | 2013-04-19 | 2013-08-07 | 黎湧 | 一种gnss与mems组合导航的方法 |
| CN103969672A (zh) * | 2014-05-14 | 2014-08-06 | 东南大学 | 一种多卫星系统与捷联惯性导航系统紧组合导航方法 |
| CN104181574A (zh) * | 2013-05-25 | 2014-12-03 | 成都国星通信有限公司 | 一种捷联惯导系统/全球导航卫星系统组合导航滤波系统及方法 |
-
2018
- 2018-06-07 CN CN201810580041.0A patent/CN108802805B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1916567A (zh) * | 2006-09-04 | 2007-02-21 | 南京航空航天大学 | 基于自适应闭环h∞滤波器的对北斗双星/捷联惯导组合导航系统进行修正的方法 |
| CN101726752A (zh) * | 2008-10-23 | 2010-06-09 | 鸿富锦精密工业（深圳）有限公司 | 地震监测系统 |
| CN102608625A (zh) * | 2012-03-30 | 2012-07-25 | 武汉大学 | 基于惯性辅助定位接收机的实时形变监测预警系统及方法 |
| CN103235328A (zh) * | 2013-04-19 | 2013-08-07 | 黎湧 | 一种gnss与mems组合导航的方法 |
| CN104181574A (zh) * | 2013-05-25 | 2014-12-03 | 成都国星通信有限公司 | 一种捷联惯导系统/全球导航卫星系统组合导航滤波系统及方法 |
| CN103969672A (zh) * | 2014-05-14 | 2014-08-06 | 东南大学 | 一种多卫星系统与捷联惯性导航系统紧组合导航方法 |
Non-Patent Citations (2)
| Title |
|---|
| JIANGHUI GENG 等: "Recovering coseismic point ground tilts from collocated high-rate GPS and accelerometers", 《GEOPHYSICAL RESEARCH LETTERS》 * |
| 李大威等: "卡尔曼滤波在弹道组合导航仿真中的应用", 《航空计算技术》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113126143A (zh) * | 2021-04-22 | 2021-07-16 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种便携式微动及强震动实时监测仪系统 |
| CN113805223A (zh) * | 2021-08-16 | 2021-12-17 | 南京天巡遥感技术研究院有限公司 | 一种地震勘探系统及其采集数据的处理方法 |
| CN114167493A (zh) * | 2021-11-23 | 2022-03-11 | 武汉大学 | Gnss双天线辅助陀螺的地震旋转测量系统及方法 |
| CN114167493B (zh) * | 2021-11-23 | 2023-08-04 | 武汉大学 | Gnss双天线辅助陀螺的地震旋转测量系统及方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108802805B (zh) | 2019-10-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107270893B (zh) | 面向不动产测量的杆臂、时间不同步误差估计与补偿方法 | |
| CN106405670B (zh) | 一种适用于捷联式海洋重力仪的重力异常数据处理方法 | |
| US10788591B2 (en) | Methods, apparatuses, and computer programs for estimating the heading of an axis of a rigid body | |
| JP5068531B2 (ja) | 測定及び記憶された重力傾度を用いて慣性航法測定値の精度を改善する方法及びシステム | |
| CN103900611B (zh) | 一种惯导天文高精度复合两位置对准及误差标定方法 | |
| CN108802805B (zh) | 六自由度宽频带一体化地震监测系统 | |
| CN103090870B (zh) | 一种基于mems传感器的航天器姿态测量方法 | |
| CN108051866A (zh) | 基于捷联惯性/gps组合辅助水平角运动隔离的重力测量方法 | |
| CN108225308A (zh) | 一种基于四元数的扩展卡尔曼滤波算法的姿态解算方法 | |
| CN110133692B (zh) | 惯导技术辅助的高精度gnss动态倾斜测量系统及方法 | |
| CN113720330B (zh) | 一种亚角秒级的遥感卫星高精度姿态确定设计与实现方法 | |
| JP4316777B2 (ja) | 重力測定装置及び方法 | |
| CN113203418A (zh) | 基于序贯卡尔曼滤波的gnssins视觉融合定位方法及系统 | |
| Lindner et al. | Seafloor ground rotation observations: Potential for improving signal‐to‐noise ratio on horizontal OBS components | |
| CN114518587A (zh) | 一种室内外无缝定位系统及方法 | |
| CN109084755B (zh) | 一种基于重力视速度与参数辨识的加速度计零偏估计方法 | |
| KR20200139613A (ko) | 방향 발견자 | |
| CN110514201A (zh) | 一种惯性导航系统及适用于高转速旋转体的导航方法 | |
| CN113551669A (zh) | 基于短基线的组合导航定位方法及装置 | |
| Veth et al. | Alignment and calibration of optical and inertial sensors using stellar observations | |
| Tu et al. | Tightly integrated processing of high-rate GPS and accelerometer observations by real-time estimation of transient baseline shifts | |
| Bhaskar et al. | Effect of oscillator quality on ultra-tight GPS/INS aided carrier phase tracking | |
| RU2779274C1 (ru) | Способ измерения ошибок начальной выставки инерциальной навигационной системы без привязки к внешним ориентирам | |
| Grejner-Brzezinska et al. | Enhanced gravity compensation for improved inertial navigation accuracy | |
| Noureldin et al. | Improving the performance of alignment processes of inertial measurement units utilizing adaptive pre-filtering methodology |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Geng Jianghui Inventor after: Wen Qiang Inventor before: Geng Jianghui Inventor before: Wen Qiang Inventor before: Li Chenghong Inventor before: Chang Hua |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |