观测系统逼真度特性分析＊

刘序俨 黄声明 梁全强 季颖锋 王 旭

1)福建省地震局,福州 350003

2)厦门地震勘察研究中心,厦门 361021

3)日本神戸大学理学研究科,神戸 657-8501

观测系统逼真度特性分析＊

刘序俨1)黄声明1)梁全强2)季颖锋1,3)王 旭1)

1)福建省地震局,福州 350003

2)厦门地震勘察研究中心,厦门 361021

3)日本神戸大学理学研究科,神戸 657-8501

在对观测系统保真度概念内涵阐述的基础上,对观测系统满足逼真度要求为什么必须为一线性不变系统的原理作了分析,指出观测系统的单位脉冲时间响应函数(权系数)和传递函数与输入的大小及其性质无关,并完全表征了系统的特性,从时域和频域对观测系统的权系数与传递函数的特征进行的分析表明：对于一个物理可实现的观测系统,不可能为一理想线性时不变系统,仅能要求在规定的通频带内和允许的观测精度范围内其传递函数近似为一平坦直线时,其相位滞后近似等于 0或与频率成正比,就能满足观测的逼真度。最后对观测系统的格值及其标定进行讨论,以水平摆为例,从其动力学常系数线性微分方程推导出理论格值,以说明一种优良的观测系统,不但要给出格值人工标定的方法,同时也要给出理论格值公式,并要经受住观测数据的可靠性检验。

观测系统;逼真度;脉冲时间响应函数;传递函数;格值与标定

1 前言

观测系统是指平常所称的观测仪器,观测仪器则是我们人类感觉器官的延伸,例如哈勃太空望远镜就是我们眼睛的延伸,使我们可观测到在宇宙深处超新星爆炸以及在奇点大爆炸后宇宙演化的奇妙景像[1];顶尖的显微镜使我们视觉深入到微观世界去观测原子、分子的结构;固体潮观测系统,例如重力仪、倾斜仪与应变仪使我们能分别发现原来固体地球宛如海潮一样也会在天体起潮力作用下发生潮起潮落现象,致使观测点的重力、地表面倾斜以及应变(线应变、剪应变、面应变与体应变)也会发生潮汐变化。重力仪能观测到由天体起潮力所引起的10-8ms-2(10-8ms-2=0.01μms-2=1微伽)的重力变化,该变化相当于重力的 10-9;石英摆与水管倾斜仪能观测由该力所导致的地面发生 0.05角秒的倾斜变化,相当于 400 m长的距离上地面发生了±0.1 mm的垂直变化;而应变仪则能观测到由该力所致的地面发生 10-8～10-9应变变化[2,3],相当于地面 100 m长的两点距离会发生 1～0.1微米的伸缩变化[2,3],上述这些几何物理量的变化我们用肉眼根本无法察觉,只有凭借仪器才能发现固体地球的潮起潮落现象。在自然界中,唯有海潮与固体潮理论值是可以预先计算出来的,根据观测结果发现,观测曲线简直是理论值曲线放大后的复制品,并且固体潮观测资料的维尼迪柯夫调和分析结果,即所谓的潮汐因子与相位迟后与理论预期值非常接近[4]。众所周知,固体潮主要是由日波、半日波与1/3日波组成,而这 3种波又分别由若干波群组成,而每一个波群又由角频率十分接近的分波叠加而成,为什么上述 3种固体潮观测系统能毫无失真地观测到由这些不同角频率的分波所组成的固体潮波形?这些观测系统到底有什么性质?为什么这些系统能对这些跨越从零频到周期为 1天的分波作出如此高逼真的响应?观测系统的响应值 (输出值)与输入值是不同量纲的,又如何对观测值进行标定,使观测值转换成我们所需要的观测量纲?以上诸问题正是本文所要探讨的。

2 线性时不变系统

所谓系统就是指能获取输入信号并可对其进行操作并产生输出信号的一种装置[5],或者说,凡是反映因果关系的任何一种物理装置或数学运算在数学上都可称之为系统。系统的特性通常由输出信号(系统响应)与输入信号 (系统激励)之间的关系来描述。系统的输入输出可以用数学法则或图表表示,也可用系统的元器件、部件符号图或框图表示,系统按它们具有的性质分类,主要有线性与非线两大类[6]。如果输出与输入成正比,且输出不失真也不延迟的系统,称为线性时不变系统。在地震系统,用以捕捉前兆信号的所有观测系统大多都可归入这种线性时不变系统。实际上,任何一个观测系统都不可能是一种真正的线性时不变系统,对于一个物理可实现的观测系统来说,仅要求在规定的通频带内并在允许的观测误差范围内能满足逼真度要求就可视为一线性时不变系统。

图 1为某一观测系统的输入输出示意图。图 1中 X(t)与 Y(t)分别为观测系统的输入与输出,它们分别为图中所示的输入与输出中相应诸分量 xi(t)与 yi(t)的叠加。

图1 系统输入与输出示意图Fig.1 Sketch of the system input and output

如果一个系统同时满足叠加性和齐次性,则称该系统为线性系统[6]。所谓叠加性是指当若干个激励同时作用在系统时,系统的总响应等于各个激励单独作用时所产生的响应之和,即

式中,yi(t)为相应于 xi(t)的输出。

所谓齐次性,指的是当系统的激励分别增大到α(α为任意常数)倍时,系统相应的响应也随之增加到α倍,α称为系统增量或放大倍数。

综合以上两个性质,一个线性系统须满足:

或

根据其系统线性性质,一个线性系统的响应应当是各个激励单独作用所产生的响应之和。

如果一个线性系统满足

则称该系统为时不变系统[6]。时不变系统表明系统的响应与时间原点无关,系统输出波形仅与输入波形有关。

在固体潮频域内,不管是重力仪、倾斜仪还是应变仪,它们皆为线性不变系统。根据固体潮理论,由杜德森关于起潮力位的展开,固体潮观测值可以写成[2]

式中:t表示自 T起算的以小时表示的时间,T是 t为 0时的时间;Hn是某潮汐波的振幅;ωn是某潮汐波的角频率,量纲为°/小时;φn(T)为在 T时角频率为ωn的某潮汐波的初相位。式中 (ωnt+φn(T))称为角频率为ωn的某潮汐分波的相位。

如果固体潮观测系统为线性时不变系统,则该系统对固体潮中的每一分波都放大了同一倍数,且系统的响应与时间原点无关,即系统输出波形完全为输入波形被放大了同一倍数的波形。因此,一个线性时不变观测系统的输出只可能包含输入中存在的频率,不会有也不可能有新的频率出现,且输出中的每一频率分量都延迟相同的时间,自然它们之和所表示的波形不变,只是发生时间延迟而已。当然,这都是在规定的通频带内以及允许的误差范围内才成立。

3 单位脉冲的响应时间函数

一般来说,一个观测系统或者某种观测仪器,都可看做是一个黑箱,因为其内部的具体结构或性能对观测者来说,都是不甚清楚的,但如何判断该观测系统的特性呢?唯一的办法是以一单位脉冲对该观测系统进行激励并作为系统的输入,然后根据该观测系统对该激励的响应即输出就可清楚地了解该系统的特性。

图 2为某个观测系统对单位脉冲δ(t)的时间响应示意图。图中δ(t)为单位脉冲作为对系统的输入,δ(t)为狄拉克函数,其定义为

式中,h(t)为该观测系统对单位脉冲δ(t)的响应,称之为单位脉冲时间响应函数,也称为权系数以作为系统的输出。不同的观测系统具有不同的权系数,因此权系数表征了观测系统的固有特性。

图 2 某个观测系统对单位脉冲δ(t)的时间响应示意图Fig.2 Ti me response curve of a certain observation system to a unit impulseδ(t)

为了求得某一离散时间线性系统对输入 x(t)的响应,我们把 x(t)表示成等间隔的以输入 x(kΔ)为振幅的一系列脉冲的叠加,如图 3所示,其表达式为

图 3输入 x(t)为脉冲串的示意图Fig.3 Sketch of inputx(t)as a pulse string

图 3中,每一个输入都作为一个振幅不等的脉冲,观测系统对每一个输入脉冲激励所作出的响应不会如输入那样也是一个嘎然而止的脉冲,往往表现为一串逐渐衰减的波形,如图 2所示,这是系统惯性使然,使之产生拖尾效应,使输入 x(t)的一系列前行值 x(t-kΔ)(k=1,2,3,…,n)会对输入的现行值 x(t)分别作出程度不等的贡献,例如前行值 x(t -Δ)的贡献为 x(t-Δ)h(Δ),x(t-2Δ)的贡献为 x (t-2Δ)h(2Δ)…最后系统对输入现行值 x(t)产生的输出响应值为[5,6]:

式中,Δ为采样的时间间隔,“＊”表示褶积。如果系统对单位脉冲的时间响应函数或权系数亦是一个在起始刹那间嘎然而止的脉冲,即除起始权系数h(0)以外,其余的权系数皆为零,此时式 (8)就变为:

但任何一种系统的响应不可能是一种嘎然而止的脉冲,因物理系统总存在惯性,使现行值 x(t)的输出值 y(t)中或多或少总含有一系列前行输入值的影响。如果在式(8)中所有前行值加权和的影响可忽略的话,那么式 (9)就可近似成立。在控制论中,如输出与输入成正比,且输出不失真也不延迟反映输入的环节称为比例或放大环节[7,8],式 (9)即为该环节的动力学方程,式中 h(0)为一常数,称为该环节的放大倍数或增益。如何评价观测系统输出波形的逼真度呢?换句话说,如何验证输出相对于输入没有失真呢?要做到这一点,只要依据式 (8),把对系统的输入时间序列值 x(t)与相应输出值 y(t)的时间序列波形进行对比,就可对输出波形是否与输入波形保持高逼真度作出判断,不过,此种判断方法要求在整个观测时间窗内对输出与输入的波形进行比较,而不能对输出与输入中的每一频率分量的振幅比是否保持为常数,两者的相位滞后是否为零或相位滞后是否与频率成正比作出检验。

4 传递函数

要解决对输出与输入中同一频率分量的振幅比是否为一常数且两者相位滞后是否为零或其与频率之比为一常数,这就要从时间域转换到频率域对输入与输出时间序列的频谱进行分析。

根据傅立叶变换理论,对式 (8)两边取傅立叶变换,有[5,6,9]

式中,X(ω)与 Y(ω)分别为输入与输出时间序列的傅立叶变换,H(ω)为系统的传递函数,H(ω)与 h (t)互为傅立叶变换对[6,9,10],可分别表达为:

式中 H(ω)为一复数,且可表示为

式中,A(ω)称为系统的振幅谱或系统放大倍数, φ(ω)称为系统相位谱或相位滞后。相位滞后是指输出与输入的相位差。由式(10)可得:

单位脉冲的时间响应函数与传递函数互为傅立叶变换对,他们分别从时域与频域决定了观测系统的特性。这是因为单位脉冲的无限小的时间间隔必然具有无限大的频带宽度,且其频谱是平坦的,如图4所示。另一方面,它又暗示在这无限大的频域内,无数多的频率分量必然在时间原点同相,因此,脉冲的相位谱处处为零。如果系统对单位脉冲的响应亦是一个嘎然而止的脉冲,则该脉冲的振幅即为式(8)中的 h(0),而该脉冲的频谱亦如图 4的单位脉冲一样也是平坦的,不过该频谱的振幅为 h(0),此时系统的传递函数即为 h(0)。因此,系统对单位脉冲的嘎然而止的脉冲响应的振幅就是系统的传递函数,在这种情况下,观测当然是完全逼真的,但这仅仅是一种理想情况。实际上时间间隔为无限小的脉冲是不存在的,因此式 (8)不可能简化为式 (9),输出的现行值 y(t)仅能为由式(8)所示的输入的现行值 x(t)与有限项的前行值 x(t-k△)(k=1,2,…, n)的加权和进行表达,并仅能要求在期望的通频带内观测系统的传递函数在精度允许的范围内可视为平坦,且在该通频带内具有零相位谱或相位滞后与频率成正比,则该观测系统就具有高逼真度的观测特性。

图4 单位脉冲的频谱Fig.4 Frequency spectrum of unit-impulse

在式(12)中,如果 A(ω)=常数,φ(ω)=0或φ(ω)/ω为常数,则该系统就能保证输出与输入的相似性,但仅仅要求波形相似是不够的,因为相似仅是指输出与输入的形状相似,而不要求两者大小相等[11]。为了做到这一点,就必须对观测值进行标定,只有经标定取得观测格值后,方能把观测值转换成所要求的物理量纲,以获得可靠的观测数据。

5 格值与标定

传递函数 H(ω)虽然表征了观测系统的固有特征,但它却无法提供有关观测系统物理结构的任何信息,因为在物理学中许多不同系统的传递函数可以是一样的,但它却包括了一些把输入与输出联系起来的单位,从而提供传递函数 H(ω)的振幅谱量纲,该量纲是由输入与输出两者的量纲确定的。例如光杠杆式水平摆倾斜仪,其输入为地面倾斜,其量纲为角秒,输出为水平摆杆上的反光镜随地面倾斜产生的光点移动,其量纲为米,则水平摆的振幅谱A(ω)的量纲为角秒 /米,A(ω)称为水平摆这个观测系统的格值,其倒数称水平摆的灵敏度,格值越小,灵敏度越高,表明水平摆在同样的地倾斜下,其光点移动量越大,表明该摆具有极高的放大能力。目前,石英摆倾斜或水管倾斜,其格值 H(ω)可达到 1 ms/ mm,或者说地面若倾斜 0.001″(1 ms),记录笔在纸介质上移动可达 1 mm,若考虑到倾斜固体潮最大振幅可达 0.05″,那么光点移动可达 50 mm,当然可观测到倾斜固体潮。因此,一个观测系统的格值是由该系统的传递函数的振幅谱给出的。对于一个观测系统而言,确定其格值是非常重要的,因为它牵涉如何由输出值确定所观测某一物理量的大小,其倒数则反映了该系统的灵敏度,即它的放大倍数。

既然格值的确定对于一个线性系统如此重要,但如何确定格值呢?从理论上讲,有两种办法,一种是基于一个逼真度高的观测系统,它必然是一个线性时不变系统,则它必须满足式 (2)、(3)所示的叠加性和齐次性条件,因此输入值 x(t)与输出值 y(t)之间必存在如下线性关系:

式中,η为格值,采用最小二法,由 x(t)与 y(t)的时间序列可确定格值η;另一种方法是基于该观测系统的传递函数的振幅谱在规定的通频带内为一常数,该常数即为格值,由式(9)可得:

因此观测系统观测格值可通过输入值 x(t)与输出值 y(t)的时间序列的傅立叶变换求得。

但可惜的是,我们无法获取输入值 x(t),因为它正是我们所要观测的一个几何物理量,我们无法由上述两种方法来确定格值。一般可行的方法是采用人工标定方法,就是由人工提供被测对象的一系列输入值,相应地我们可得到一系列输出值 (观测值),然后,根据输入与输出值就可由式(14)或(15)求得系统的格值η,例如水平摆倾斜仪采用胀盒标定法,该胀盒放置在水平摆倾斜仪的旁移脚螺丝下,该胀盒与一盛有水银的容器相连通,并悬挂在绕一水平轴作圆周运动的一个圆盘上,随着周期运动,使水银高度产生周期性变化,从而使水银对胀盒施加一个周期性的作用力,使胀盒随之产生一种周期性的伸缩变化,该伸缩变化可由胀盒的力学参数计算出来,从而人工地对水平摆倾斜施加一个地面倾斜输入[3,11];对于水管倾斜仪,是采用标定棒法进行人工格值标定,其主要原理是借助标定棒的上升和下降来增减仪器内的液体体积,从而使液体水位发生变化以转换成电信号来进行仪器标定[11],对于洞体伸缩仪进行标定,是在石英杆基线自由端安置数显电感器,在基线的固定端安置标定装置,人工给出一系列已知的位移值,安置在自由端的数显电感器就可将此位移转换为电信号,由基线长度和位移就可计算出应变,从而可对电信号观测值进行格值标定[12,13]。

由于观测系统的工作状态不是一成不变的,同时台站环境也可能爱到外界的干扰,因此,在平时工作中,要定期对仪器进行标定,以保证观测数据的可靠性。

6 认识与讨论

本文围绕为满足逼真度的要求,对观测系统为什么应该是一个线性时不变系统并对系统输入与输出的关系,单位脉冲时间响应函数与传递函数的物理含义,观测格值标定等问题进行了深入分析并着重指出,一个观测系统特性可由输入与输出来确定,观测系统对单位脉冲激励所作的响应 (输出)即单位脉冲时间响应函数及传递函数决定了该系统的特性,他们与输入的大小及其性质无关。因此,一个线性时不变系统的单位脉冲时间响应函数与传递函数包含了有关该系统的动力学同样多的信息[14],这是不难理解的,因为系统的单位脉冲时间函数与传递函数两者互为傅立叶变换,前者适用于时间域,后者适用于频率域,两种是等价的。传递函数作为观测格值的数学表达式,任何一家仪器研制部门和厂家都应向用户提供格值公式。

下面以水平摆倾斜仪为例,以说明理论格值公式是可以从仪器的动力学方程推导出来的。下式为水平摆动力学常系数线性微分方程①刘序俨,李华,鲍挺.水平摆纵横谈 [J].形变学科通讯,2002, (2)：51-59;2003,(1)：54-63.:

式中:β为在时刻 t时的地面倾斜值,α为水平摆相应偏转角度,α、β皆为微小角度,以弧度为单位;ξ为水平摆阻尼比,为无量纲单位;ω为水平摆无阻尼角速度,以 1/s为单位;m为水平摆锤质量,以 kg为单位;l为水平摆折合摆长,以m为单位;g为重力加速度,以ms-2为单位;I为水平摆转动惯量,单位为kgm2。

设∑为地面倾斜的角频率,Ω为水平摆的角频率,则ω=2π Ω,对式 (16)两边取傅立叶变换,根据傅氏变换定理[6,9],有

由式(10)可得到水平摆的传递函数为:

式中(i+ε)以弧度为单位。

对于某些观测系统,例如承压井水位观测系统,由于它本身就是一个天然体应变计,无法用人工提供地壳体应变对其进行激励,因此也就无法对承压井水位观测系统本身进行格值标定,但井水位既然为一天然体应变计,我们就可利用固体地球的体应变理论值作为激励值,即作为人工的输入值,从而可以采用式 (10)计算井水位的传递函数,不过在这里,传递函数之倒数才是井水位观测的格值[15],把此格值乘以井水位值就可求得体应变值。对井水位观测系统以体应变理论值作为人工激励,是不得已而为之,当然,我们亦可推导出承压井水位观测系统的理论格值表达式[15],但因在该表达式中含有含水层孔隙度这个未知参数,我们还是无法获得事实上的理论格值。对于那些非天然的固体潮观测系统,无论是重力、倾斜、应变 (包括体应变)观测系统都不能这样做,都应给出人工标定的方法,同时亦应给出观测系统的传递函数,以便校核。同时,我们还可利用维尼迪柯夫调和分析方法对上述观测资料进行调和分析以取得潮汐因子与相位滞后,如果各波群的潮汐因子彼此非常接近且与理论预期值十分一致,且相位滞后为负值且接近于零,即就可对观测格值的标定正确性作出判断,且潮汐因子的倒数即是格值。

最后,在这里要着重指出,虽然上述诸种观测系统能清晰地记录到日波至 1/3日波频段上的地球固体潮汐,但对地壳长周期运动进行监测却无能为力,究竟原因,这是由于这些观测系统存在系统误差不确定性所引起的,产生此种不确定性大致有以下两方面的原因:一方面,这些观测系统本身就存在结构或电气方面的漂移,表现在即使安装在同一仪器墩子上的同样两台仪器也有可能出现趋势方向不同的变化,但在各自趋势变化的背景上所记录到的固体潮波形是十分一致的;另一方面,这些地面观测系统易受到地壳红噪声的干扰,地壳“红噪”谱与其频率成倒幂律关系[16],地壳振动频率越低,其噪声越大,在这种“红噪”背景下,当然无法对地壳长周期信号进行监测。存在系统误差不确定性是此类地面观测系统固有的缺陷,如何克服这种缺陷是摆在仪器研制者面前的一项重要课题。

1 徐世忠.哈勃空间望远镜探秘[M].北京：科学出版社, 2000.

2 北京大学地球物理系,武汉测绘学院大地测量系,中国科学技术大学地球物理教研室.重力与固体潮教程[M].北京：地震出版社,1982.

3 梅尔基奥尔(杜品仁,等译).行星地球的固体潮[M].北京：科学出版社,1984.

4 Melchior,等.中国固体潮观测研究 [J].地球物理学报, 1985,28(2)：142-154.

5 Zoher Z Karu.Signals and systems[M].ZiZi Press,Cambridge,MA,2001.

6 林秩盛.信号与线性系统[M].北京：清华大学出版社, 2008.

7 曾燕主编.控制工程基础[M].北京：电子工业出版社, 2007.

8 陆一心主编.现代工程控制理论[M].北京：化学工业出版社,2006.

9 奥本海姆,谢弗著.董大嘉,杨增辉译.茅于海,等校.数字信号处理[M].北京：科学出版社,1983.

10 Sanjit K Mitra.Digital signal processing—A computerbased approachse[M].北京：清华大学出版社,2006.

11 Clapham C.Oxford concise dictionary of mathematics[M]. Shanghai：Shanghai Foreign Language Education Press, 2001.

12 国家地震局科技监测司.地震形变观测技术[M].北京：地震出版社,1995.

13 杜为民,刘序俨.定点形变仪器的标定[J].地壳形变与地震,1989,(3)：85-92.

14 Katsuhiko Ogata.Modern control engineering(fourth edition)[M].北京：电子出版社,2007.

15 刘序俨,等.承压井水位观测系统对体应变响应机制分析[J].地球物理学报,2009,52(12)：3 147-3 157.

16 Perron J T,Kirchner J W and DietrichW E(毕丽思,何宏林译).地形中的特征空间尺度与非分形构造的频谱信号[J].世界地震译丛,2009,(5)：39-54.

FIDEL ITY ANALYSIS OF AN OBSERVATION SYSTEM

Liu Xuyan1),Huang Shengming1),LiangQuanqiang2),Ji Yingfeng1,3)andWang Xu1)

(1)Earthquake Adm inistration of Fujian Province,Fuzhou 350003 2)Xianm en Research Centre of Earthquake Surveging,X iam en 361021 3)Faculty of Science,Kobe University,Rokkotaim achi1-1,N ade,Kobe,Japan,Postcode657-8501)

On the basis of describing the concept of an observation system’s fidelity,it takes a close look on the principle aboutwhy the observational system which meets the requirement of fidelity is sure to be an invariable linear system and points out that for an observation system the response function of unit pulse ti me as well as the transfer function area irrelevantwith the quantity and quality of the input.These functions are capable of totally representing the characteristics of the observation system.Besides,by studying the weight coefficient of observation system from the time-domain and frequency-domain alongwith the characteristics of transfer function,it concludes that a physically realizable observation system for sure will fail to be an ideal linear time-invariant system.Nevertheless,given an appointed transmissional band and a certain requirement forprecision,the transfer function curves will turn similarly into a level line with a phase lag approximately to zero or directly related to the frequency,which can largelymeet the demands observation.Finally,it discusses the scale unit and its calibration of observation system,by analyzing a example of horizontal pendulum that derives a theoretical scale unit from the dynamic constantcoefficient linear differential equations,to demonstrate that an excellent observation system should not only has a manualmethod for scale unit calibration,but also a theoretical formula for the scale unit deriving and it is dependable in any reliability examinations ofwhatever observation data.

observation system;fidelity;response function of unit impulse ti me;transfer function;scale unit and calibration

1671-5942(2010)Supp.(Ⅰ)-0001-07

2010-05-12

中国地震局老专家科研基金

刘序俨,男,研究员,长期从事固体潮与地形变研究.E-mail：xuyanliu@126.com

TH76

A