一種基于扭秤的重力梯度動(dòng)態(tài)測(cè)量方法

葉子薇, 王 飛, 包 福, 黃安貽,余 燁, 胡紅波, 王若琳

(1.武漢理工大學(xué) 機(jī)電工程學(xué)院，湖北 武漢 430070；2.湖北省計(jì)量測(cè)試技術(shù)研究院，湖北 武漢 430223；3.武漢光谷航天三江激光產(chǎn)業(yè)技術(shù)研究院有限公司，湖北 武漢 430073)

1 引 言

本文針對(duì)重力梯度儀在動(dòng)態(tài)測(cè)量研究中因結(jié)構(gòu)限制只能測(cè)得重力水平梯度Wxz和Wyz這一現(xiàn)狀,將Z字型秤桿設(shè)計(jì)引入重力梯度動(dòng)態(tài)測(cè)量,這種扭秤關(guān)于中心對(duì)稱,結(jié)構(gòu)簡(jiǎn)單,能夠通過動(dòng)態(tài)測(cè)量同時(shí)獲取重力水平梯度和重力曲率。相比厄缶扭秤通過細(xì)絲懸掛檢驗(yàn)質(zhì)量的方式,直接將豎桿與秤桿粘貼,大大減小懸掛系統(tǒng)的穩(wěn)定時(shí)間,同時(shí)也避免了檢驗(yàn)質(zhì)量微小擺動(dòng)對(duì)重力梯度測(cè)量的影響[18]。此外,由于APES算法處理動(dòng)態(tài)測(cè)量信號(hào)的精度仍存在改進(jìn)的空間,在動(dòng)態(tài)測(cè)量信號(hào)后處理的過程中采用遺傳算法,實(shí)現(xiàn)更高準(zhǔn)確度的動(dòng)態(tài)測(cè)量信號(hào)處理。

2 扭秤法測(cè)量原理

扭秤是一種經(jīng)典的弱力測(cè)量工具,由一根長(zhǎng)細(xì)絲懸掛,在水平方向自由旋轉(zhuǎn)。由于對(duì)水平方向的外力極其敏感,因此可實(shí)現(xiàn)扭秤沿水平方向扭矩的測(cè)量。本文設(shè)計(jì)的扭秤式重力梯度測(cè)量裝置,其Z字型扭秤結(jié)構(gòu)如圖1所示。

圖1 扭秤基本結(jié)構(gòu)示意圖

以秤桿質(zhì)心為原點(diǎn)建立直角坐標(biāo)系,xyz坐標(biāo)系為地固坐標(biāo)系,x軸與y軸在水平面上,z軸沿細(xì)絲方向垂直向下。ξηz坐標(biāo)系為固定在扭秤上的扭秤坐標(biāo)系,ξ沿著秤桿長(zhǎng)方向,η軸為垂直秤桿方向[19]。秤桿與x軸之間的夾角為α,需測(cè)量的量為多個(gè)方位角α處懸絲受重力梯度影響時(shí)的扭秤偏轉(zhuǎn)角度θ。

設(shè)扭秤系統(tǒng)的某1個(gè)質(zhì)量單元為dm,對(duì)應(yīng)的力矩平衡方程為[20]:

-f2Wxzsinα+f2Wyzcosα

(1)

由于地球上不存在重力場(chǎng)均勻的點(diǎn),因此θ0為未知量,Wyz、Wxz、Wxy、WΔ這4個(gè)重力梯度也為未知量,若想解出這5個(gè)未知量,則需要5個(gè)方程。因此,需要在同一點(diǎn)設(shè)置不同的方向上測(cè)量5次,即可得到方程求解未知量。

動(dòng)態(tài)測(cè)量重力梯度時(shí),將扭秤放置在勻速轉(zhuǎn)動(dòng)的平臺(tái)上,在初始測(cè)量狀態(tài)下,扭秤坐標(biāo)系與地固坐標(biāo)系重合。啟動(dòng)勻速轉(zhuǎn)臺(tái),整體裝置與勻速轉(zhuǎn)臺(tái)同步以ω的速度旋轉(zhuǎn),此時(shí),扭秤系統(tǒng)的方位角α(t)=ωt,扭秤偏轉(zhuǎn)角度θ也是隨時(shí)間變化的,用θ(t)表示。待信號(hào)平穩(wěn)時(shí),扭秤標(biāo)準(zhǔn)平衡方程為:

(2)

由非線性振動(dòng)理論可知,在扭秤偏轉(zhuǎn)角度信號(hào)θ(t)中,除了基頻ω和二倍頻2ω的角速度運(yùn)動(dòng)信號(hào),還包含了扭秤的自由振蕩運(yùn)動(dòng)信號(hào)[21]。因此,扭秤偏轉(zhuǎn)角度信號(hào)模型為:

θ(t)=A0(t)cos[ω0(t)t+φ0]

(3)

式中:A0(t)、ω0、φ0分別為自由振蕩信號(hào)的時(shí)變振幅、時(shí)變頻率、初相位;Ak為信號(hào)k倍頻分量的幅值;ω為信號(hào)的基頻頻率;φk為k倍頻分量的初相位;n(t)為高斯白噪聲。

+a1cos(ωt)+b1sin(ωt)+a2cos(2ωt)

+b2sin(2ωt)+n(t)

(4)

式中:A0為自由振蕩初始振幅;β為振幅衰減因子;為自由振蕩平均周期;T為數(shù)據(jù)采集時(shí)間;ε為自由振蕩周期變化范圍。

扭秤偏轉(zhuǎn)角度信號(hào)數(shù)據(jù)處理的目標(biāo)量是基頻ω和二倍頻2ω的余弦分量幅值a1、a2及正弦分量幅值b1、b2。結(jié)合式(2)和式(4),可得到a1、b1、a2、b2與待測(cè)梯度值Wyz、Wxz、Wxy、WΔ滿足式(5)所示關(guān)系[21]。

(5)

扭秤結(jié)構(gòu)設(shè)計(jì)為Z字型時(shí),f1和f2均不為零,通過測(cè)量θ(t),由式(5)可以計(jì)算出重力水平梯度Wyz、Wxz和重力曲率Wxy、WΔ。

3 數(shù)據(jù)處理方法

3.1 旋轉(zhuǎn)濾波

實(shí)驗(yàn)時(shí)由于溫度、濕度、噪聲等環(huán)境的環(huán)境因素影響,扭秤偏轉(zhuǎn)角度信號(hào)中的自由振蕩信號(hào)為頻率時(shí)變信號(hào)。若直接對(duì)θ(t)采用多重回歸分析處理,頻率時(shí)變的自由振蕩信號(hào)的存在,會(huì)導(dǎo)致ω0(t)中的參數(shù)無法準(zhǔn)確求解。此外,自由振蕩信號(hào)還會(huì)影響對(duì)基頻和二倍頻的參數(shù)提取。因此,需要剔除自由振蕩信號(hào),使θ(t)滿足多重回歸分析的要求,提高參數(shù)提取的準(zhǔn)確度。

旋轉(zhuǎn)濾波可去除單一頻率信號(hào)[22],因此使用旋轉(zhuǎn)濾波對(duì)自由振蕩信號(hào)進(jìn)行剔除。以自由振蕩平均周期T0=2π/ω0為時(shí)間間隔對(duì)信號(hào)進(jìn)行180°的翻轉(zhuǎn),即有翻轉(zhuǎn)后的信號(hào)為:

(6)

式(6)中T0可通過Quinn頻率估計(jì)法得到,將式(3)代入式(6),并通過三角函數(shù)性質(zhì)化簡(jiǎn)可得:

(7)

式中:Aj為信號(hào)j倍頻分量的振幅;jω為信號(hào)j倍頻分量的頻率;φj為信號(hào)j倍頻分量的初相位。

由式(7)可知,若能準(zhǔn)確評(píng)估自由振蕩平均周期,則θ(t)中自由振蕩信號(hào)能被完全剔除。

3.2 Quinn頻率估計(jì)方法

采用Quinn頻率估計(jì)方法提取整段自由振蕩信號(hào)的平均周期T0。該方法的主要優(yōu)點(diǎn)是算法簡(jiǎn)單,計(jì)算量小,可估計(jì)整段目標(biāo)信號(hào)的平均周期,滿足提取T0的要求。

3.3 APES算法

使用APES算法獲取相位和幅值的步驟如下:

1) 通過旋轉(zhuǎn)濾波法,去除θ(t)的主要干擾部分即自由振蕩信號(hào);

2) 利用Matlab進(jìn)行基于APES算法的實(shí)信號(hào)幅值和相位提取過程的數(shù)值計(jì)算實(shí)驗(yàn)。

3) 在2～100調(diào)節(jié)濾波器階數(shù)參量,得到APES算法提取4個(gè)目標(biāo)梯度值的誤差值,對(duì)比分析,選擇使梯度值誤差最小的濾波器階數(shù)為APES的最佳濾波器階數(shù)。

3.4 遺傳算法

θ(t)中基頻和二倍頻相隔較近,相互影響,所以在遺傳算法處理過程中采用如下先構(gòu)造并去除干擾量后擬合目標(biāo)量的方式:

1) 通過旋轉(zhuǎn)濾波,去除θ(t)的主要干擾部分即自由振蕩信號(hào)。

2) 在求基頻的幅值和相位時(shí),二倍頻為干擾量。先以A2sin(2ωt+φ2)為目標(biāo)函數(shù)構(gòu)造二倍頻分量,自適應(yīng)函數(shù)為扣除自由振蕩信號(hào)后的角度信號(hào)θ(t)減去目標(biāo)函數(shù)后的幅值譜在2ω頻率附近信號(hào)幅值的平均值。通過遺傳算法,得到使自適應(yīng)函數(shù)最小的參數(shù),將參數(shù)帶入目標(biāo)函數(shù)模型,構(gòu)成二倍頻分量,并將其從原信號(hào)中減去,用來消除二倍頻對(duì)基頻的影響。然后以A1sin(ωt+φ1)為目標(biāo)函數(shù)擬合基頻分量,評(píng)價(jià)函數(shù)為扣除自由振蕩信號(hào)和二倍頻信號(hào)后的角度信號(hào)θ(t)減去目標(biāo)函數(shù)后的幅值譜在ω頻率附近信號(hào)幅值的平均值。通過遺傳算法,得到使評(píng)價(jià)函數(shù)最小的參數(shù),這兩個(gè)參數(shù)即為所求的基頻幅值和相位。

3) 在求二倍頻的幅值和相位時(shí),與步驟2)同理,先消除基頻分量的影響。然后通過遺傳算法擬合二倍頻分量,從而得到所求的二倍頻幅值和相位。

4) 得到目標(biāo)量基頻與二倍頻的幅值和相位,由式(5)可得到4個(gè)待測(cè)梯度值,并分別計(jì)算誤差。

4 仿真實(shí)驗(yàn)

4.1 動(dòng)態(tài)仿真模型

根據(jù)扭秤法動(dòng)態(tài)測(cè)量原理進(jìn)行仿真模型設(shè)計(jì),以MATLAB為工具進(jìn)行仿真,通過數(shù)值計(jì)算實(shí)驗(yàn)比較不同算法的數(shù)據(jù)處理結(jié)果。設(shè)計(jì)系統(tǒng)的整體示意圖如圖2所示,系統(tǒng)主要由二維導(dǎo)引、鎢絲、Z字型秤桿、石英窗口、光電自準(zhǔn)直儀、旋轉(zhuǎn)轉(zhuǎn)臺(tái)和密度均勻質(zhì)量塊組成。二維導(dǎo)引與旋轉(zhuǎn)轉(zhuǎn)臺(tái)共軸,轉(zhuǎn)臺(tái)以ω的角速度旋轉(zhuǎn),扭秤相對(duì)于轉(zhuǎn)臺(tái)的旋轉(zhuǎn)角度為θ(t)。

圖2 整體裝置示意圖

測(cè)量時(shí)為了消除動(dòng)態(tài)測(cè)量中慣性對(duì)于轉(zhuǎn)動(dòng)的影響,旋轉(zhuǎn)速度越小越好。但考慮到轉(zhuǎn)臺(tái)旋轉(zhuǎn)速度決定數(shù)據(jù)采集時(shí)間,選擇0.3°/s的旋轉(zhuǎn)速度,即轉(zhuǎn)臺(tái)旋轉(zhuǎn)周期為1 200 s。

根據(jù)第2節(jié)分析的動(dòng)態(tài)測(cè)量的扭秤偏轉(zhuǎn)角度數(shù)據(jù)模型,如式(4)所示。鎢絲直徑取50 μm,則秤桿轉(zhuǎn)動(dòng)慣量為1.277×10-3kg·m2,可以得出自由振蕩周期為T0=303.5 s[23]。設(shè)置自由振蕩信號(hào)幅值的幅值衰減系數(shù)為10-6,自由振蕩信號(hào)初始幅值為2 000角秒,初相位為0.1 rad?？紤]到T0在 162 000 s 內(nèi)的時(shí)變范圍大約為2 s,設(shè)置數(shù)據(jù)采集時(shí)間T=162 000 s,ε=2 s。

在仿真實(shí)驗(yàn)中,假設(shè)扭秤處的重力梯度全部來源于右側(cè)質(zhì)量塊所產(chǎn)生的人工梯度場(chǎng),則其理論梯度可通過計(jì)算得到。設(shè)置質(zhì)量塊重120 kg,長(zhǎng)寬高分別為420 mm、190 mm、290 mm。如圖2所示,定義質(zhì)量塊長(zhǎng)方向?yàn)閤軸,寬方向?yàn)閥軸,高方向?yàn)閦軸建立地固坐標(biāo)系。設(shè)秤桿中心相對(duì)于質(zhì)量塊質(zhì)心坐標(biāo)為40.0 mm,454.0 mm,-159.0 mm。得到人工梯度場(chǎng)的理論梯度值為Wyz=53 E,Wxz=4 E,Wxy=-11 E,WΔ=-156 E。根據(jù)式(5)得到,理論上目標(biāo)量a1=52 μrad、b1=-4 μrad、a2=-26 μrad、b2=-181 μrad。仿真模型參數(shù)表如表1所示。

表1 仿真模型參數(shù)表

建立仿真數(shù)據(jù)的數(shù)學(xué)模型:

+52cos(ωt)-4sin(ωt)-26cos(2ωt)

-181sin(2ωt)+n(t)

(8)

式中n(t)的方差為0.25。

根據(jù)式(8),在matlab中以時(shí)間間隔1 s生成1組162 000 s的仿真數(shù)據(jù)。

4.2 自由振蕩信號(hào)剔除效果

對(duì)仿真數(shù)據(jù)進(jìn)行Quinn頻率估計(jì),可得自由振蕩信號(hào)的平均周期T0=304.8 s。

對(duì)仿真數(shù)據(jù)進(jìn)行旋轉(zhuǎn)濾波,濾波前后對(duì)比圖如圖3所示?？梢钥吹酵ㄟ^旋轉(zhuǎn)濾波,自由振蕩信號(hào)的幅值減小46.71 dB。

圖3 旋轉(zhuǎn)濾波前后頻譜對(duì)比圖

4.3 APES算法數(shù)據(jù)處理

由于APES算法的提取效果與其構(gòu)造的有限長(zhǎng)單位沖激響應(yīng)(finite impulse response,FIR)濾波器的階數(shù)M有關(guān)。當(dāng)APES算法的M在2～100的范圍變化,所得4個(gè)梯度的誤差值如圖4所示。提取梯度值的最大誤差最小時(shí),所用的M為最佳濾波器階數(shù)。

圖4 不同濾波器階數(shù)下梯度值的誤差

取最佳濾波器階數(shù),進(jìn)行100次蒙特卡洛實(shí)驗(yàn),結(jié)果如圖5所示,可以看出APES算法穩(wěn)定性較好,但參數(shù)提取效果較差,APES提取梯度值的最大誤差為-0.82 E。

圖5 最佳濾波器階數(shù)下APES提取梯度值的誤差

4.4 遺傳算法數(shù)據(jù)處理

使用第3.4節(jié)中的仿真原理及步驟,對(duì)仿真模型進(jìn)行遺傳算法處理,設(shè)置蒙特卡洛循環(huán)次數(shù)為100次,結(jié)果如圖6所示。可以看出遺傳算法因求解過程帶有一定的隨機(jī)性,所以穩(wěn)定性一般,但是提取梯度值的最大誤差為0.013 3 E,準(zhǔn)確度比APES算法高出近2個(gè)數(shù)量級(jí)。

圖6 遺傳算法提取梯度值的誤差

4.5 數(shù)據(jù)結(jié)果對(duì)比分析

使用APES算法和遺傳算法提取4個(gè)目標(biāo)梯度值Wyz、Wxz、Wxy、WΔ的誤差如表2所示。APES算法提取梯度值的最大誤差為-0.82 E,遺傳算法提取梯度值的最大誤差為0.013 3 E。因此,在該重力梯度儀動(dòng)態(tài)測(cè)量模式的信號(hào)處理上,使用遺傳算法,可以顯著提高梯度值的測(cè)量準(zhǔn)確度,數(shù)據(jù)處理效果優(yōu)于APES算法。

表2 APES算法和遺傳算法仿真梯度值誤差對(duì)比

5 不確定度分析

5.1 扭秤結(jié)構(gòu)誤差引入的不確定度

實(shí)驗(yàn)使用的Z字型秤桿的結(jié)構(gòu)系數(shù)f1=1.273×10-3kg·m2,f2=5.372×10-4kg·m2。由于扭秤的加工誤差、材料的非均勻性以及組裝時(shí)扭秤姿態(tài)的非理想性等問題,其扭秤結(jié)構(gòu)參數(shù)f1和f2會(huì)存在誤差,為最終的測(cè)量結(jié)果引入不確定度。在此對(duì)f1、f2引入誤差的因素進(jìn)行評(píng)估,主要有以下幾項(xiàng):

1) 加工誤差:扭秤尺寸存在10 μm的加工誤差時(shí)引入的誤差。

2) 中心定位誤差:鎢絲懸掛在小擺件的質(zhì)心位置,小擺件質(zhì)心位置存在1 mm位置偏差時(shí)引入的誤差。

3) 密度誤差:小擺件豎桿材料密度的非均勻性引入的誤差。

4) 重物傾斜姿態(tài)誤差:組裝時(shí)的非理想性導(dǎo)致重物傾斜姿態(tài)引入的誤差。

5) 整體傾斜姿態(tài)誤差:小擺件由鎢絲懸掛時(shí)整體裝置存在傾斜引入的誤差。

表3所示為影響f1、f2大小的誤差項(xiàng)及它們對(duì)f1、f2的標(biāo)準(zhǔn)差。

表3 f1、 f2的誤差項(xiàng)及誤差標(biāo)準(zhǔn)差

由表3可知f1、f2的主要誤差項(xiàng)為加工誤差,計(jì)算加工誤差為1 μm和10 μm時(shí),f1、f2的誤差引入的不確定度,結(jié)果如表4所示。

表4 加工誤差1 μm和10 μm時(shí),f1、 f2誤差引入的不確定度

扭秤結(jié)構(gòu)誤差引入的不確定度u1,取加工誤差為1 μm時(shí)的不確定度。

5.2 數(shù)據(jù)處理過程引入的不確定度

表5 用APES和GA參數(shù)提取數(shù)據(jù)處理過程引入的不確定度

5.3 合成不確定度

被測(cè)量的各輸入量之間不相關(guān),因此系統(tǒng)的測(cè)量合成不確定度uc(W)為[24]:

(9)

針對(duì)4個(gè)待測(cè)梯度值,按式(9)計(jì)算合成不確定度。采用APES算法進(jìn)行參數(shù)提取時(shí)系統(tǒng)的測(cè)量不確定度和采用遺傳算法進(jìn)行參數(shù)提取時(shí)系統(tǒng)的測(cè)量不確定度如表6所示。

表6 用APES和GA參數(shù)提取時(shí)系統(tǒng)的測(cè)量不確定度

6 總 結(jié)

設(shè)計(jì)了一款Z字型扭秤式重力梯度測(cè)試裝置,能夠通過動(dòng)態(tài)測(cè)量同時(shí)獲取重力水平梯度和重力曲率。針對(duì)其動(dòng)態(tài)測(cè)量模式下的角度信號(hào),提出用遺傳算法對(duì)其進(jìn)行信號(hào)后處理。通過仿真分析,在剔除自由振蕩信號(hào)的條件下,證明在這種Z字型扭秤式重力梯度儀中遺傳算法的準(zhǔn)確度比APES算法高出近2個(gè)數(shù)量級(jí)。

為評(píng)估系統(tǒng)的測(cè)量不確定度,除了對(duì)數(shù)據(jù)處理算法引入的不確定度進(jìn)行分析,還對(duì)扭秤的加工誤差、材料的非均勻性以及組裝時(shí)扭秤姿態(tài)的非理想性等問題引入的不確定度進(jìn)行分析。在使用APES算法進(jìn)行數(shù)據(jù)處理和使用遺傳算法進(jìn)行數(shù)據(jù)處理的情況下,系統(tǒng)的最大測(cè)量不確定度分別為0.58 E和0.058 2 E。由此得出,在扭秤式重力梯度儀動(dòng)態(tài)測(cè)量實(shí)驗(yàn)的數(shù)據(jù)處理中,遺傳算法優(yōu)于APES算法。論證了遺傳算法在扭秤式重力梯度儀動(dòng)態(tài)測(cè)量模式的信號(hào)處理上的優(yōu)越性。所提出的重力梯度測(cè)量系統(tǒng)將Z字型秤桿結(jié)構(gòu)與動(dòng)態(tài)測(cè)量、遺傳算法相結(jié)合,從理論上實(shí)現(xiàn)了高精度的重力水平梯度和重力曲率同時(shí)測(cè)量。

未來有望優(yōu)化實(shí)驗(yàn)裝置,提高系統(tǒng)的穩(wěn)定性和準(zhǔn)確度。主要措施為:提高扭秤加工精度,減小扭秤結(jié)構(gòu)引入的不確定度分量;增添隔振裝置,降低轉(zhuǎn)臺(tái)盤面跳動(dòng)的影響;實(shí)時(shí)測(cè)量溫度、濕度、電場(chǎng)以及磁場(chǎng)等外界干擾因素,保持重力梯度在測(cè)量過程中對(duì)外界響應(yīng)一致;提高整體裝置維度,實(shí)現(xiàn)重力梯度全張量測(cè)量。由此便于在尋找油氣田[25]、航空梯度測(cè)量[26]等領(lǐng)域擴(kuò)大檢測(cè)范圍,提高測(cè)量精度,作出重要貢獻(xiàn)。