劉 杰, 閆 妍, 申夢(mèng)嫻, 劉乾坤, 羅偉東

(1.吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院, 長(zhǎng)春130061; 2.中國(guó)地質(zhì)調(diào)查局 廣州海洋地質(zhì)調(diào)查局, 廣州510000)

0 引 言

重力梯度儀是用于測(cè)量重力梯度張量的儀器, 是目前探測(cè)地球的先進(jìn)手段之一, 在能源勘探、 航空、地形恢復(fù)等方面發(fā)揮著重要的作用[1]。 重力梯度在空間上有9 個(gè)分量, 其中5 個(gè)分量是獨(dú)立的, 利用這5 個(gè)分量的數(shù)據(jù)及其變化即可以對(duì)地球上物質(zhì)發(fā)生的變化進(jìn)行預(yù)測(cè)。 目前國(guó)內(nèi)外常見的重力梯度儀主要有4 種[2]: 旋轉(zhuǎn)加速度計(jì)重力梯度儀、 靜電加速度計(jì)重力梯度儀、 超導(dǎo)重力梯度儀和冷原子重力梯度儀。筆者所研究的全張量旋轉(zhuǎn)加速度計(jì)重力梯度儀在常溫下工作, 非常適合應(yīng)用于地下資源勘探、 海洋技術(shù)以及慣性導(dǎo)航等領(lǐng)域[3-6]。 為了模擬旋轉(zhuǎn)加速度計(jì)重力梯度儀的工作原理, 研究重力梯度儀信號(hào)處理的方法, 筆者分別對(duì)全張量重力梯度儀和單圓盤重力梯度儀進(jìn)行了仿真。

1 設(shè)計(jì)方案

旋轉(zhuǎn)加速度計(jì)重力梯度儀的4 個(gè)加速度計(jì)對(duì)稱、 正交地安裝在一個(gè)低速轉(zhuǎn)臺(tái)上( GGI: Gravity Gradient Instrument), 加速度計(jì)的敏感軸沿轉(zhuǎn)臺(tái)的切線方向, 而且相對(duì)的兩個(gè)加速度計(jì)敏感軸方向相反,如圖1 所示[7]。 工作時(shí), GGI 轉(zhuǎn)臺(tái)以特定的角速率旋轉(zhuǎn), 單個(gè)GGI 轉(zhuǎn)臺(tái)能測(cè)量部分張量的重力梯度, 將3 個(gè)相同的GGI 轉(zhuǎn)臺(tái)在空間正交安裝, 構(gòu)成全張量重力梯度儀, 能實(shí)現(xiàn)全張量的重力梯度測(cè)量, 如圖2所示[6]。

圖1 旋轉(zhuǎn)加速度計(jì)重力梯度儀GGI 示意圖 Fig.1 Rotary accelerometer gravity gradiometer GGI schematic

圖2 全張量重力梯度儀結(jié)構(gòu)示意圖Fig.2 Schematic diagram of full tensor gravity gradiometer

旋轉(zhuǎn)加速度計(jì)重力梯度儀的系統(tǒng)仿真流程圖如圖3 所示。 以單圓盤為例, 計(jì)算機(jī)A 模擬正交放置的4 個(gè)加速度計(jì)輸出數(shù)據(jù), 通過D/ A 轉(zhuǎn)換為模擬信號(hào), 并對(duì)4 路加速度計(jì)信號(hào)進(jìn)行求和求差運(yùn)算, 再通過A/ D 將調(diào)理過的信號(hào)經(jīng)過串口發(fā)送到計(jì)算機(jī)B 中, 由計(jì)算機(jī)B 對(duì)數(shù)據(jù)進(jìn)行濾波和解調(diào)。

圖3 重力梯度儀仿真系統(tǒng)框圖Fig.3 Gravity gradiometer simulation system block diagram

2 信號(hào)產(chǎn)生

2.1 重力梯度測(cè)量原理

空間質(zhì)量受到的重力是由地球引力和地球自轉(zhuǎn)引起的慣性離心力構(gòu)成, 重力位函數(shù)是一個(gè)標(biāo)量函數(shù), 它對(duì)各個(gè)方向的偏導(dǎo)數(shù)等于重力在對(duì)應(yīng)坐標(biāo)軸上的分量。 重力位函數(shù)為

其中ωe是地球自轉(zhuǎn)角速度,r是空間質(zhì)量到地心的距離,x,y分別是空間質(zhì)量在固地坐標(biāo)系相應(yīng)坐標(biāo)軸上的分量,M是地球的質(zhì)量,G是萬有引力常數(shù)。

重力加速度是重力位沿著重力方向的一階偏導(dǎo)數(shù), 重力梯度是重力位沿重力方向的二階偏導(dǎo)數(shù)。 各個(gè)方向的重力加速度為

對(duì)重力位在各個(gè)方向上求二階導(dǎo)數(shù)可得到重力梯度

其中Γxy=Γyx,Γxz=Γzx,Γyz=Γzy,Γxx+Γyy+Γzz=0。

2.2 加速度計(jì)測(cè)量原理

如果忽略地球自轉(zhuǎn)角速度重力梯度儀平臺(tái)相對(duì)于地球坐標(biāo)系的轉(zhuǎn)動(dòng)角速度和圓盤自轉(zhuǎn)角速度, 根據(jù)擺式加速度計(jì)輸出的數(shù)學(xué)模型[5-8]可導(dǎo)出單個(gè)加速度計(jì)的輸出結(jié)果為

水平圓盤的輸出為E=(E1+E2)-(E3+E4), 在理想條件下, 當(dāng)4 個(gè)加速度計(jì)的標(biāo)度因數(shù)一致時(shí), 圓盤輸出為E=4KRcos(2ωt)(-Γxy)+2KRsin(2ωt)(Γxx-Γyy)。

下面推導(dǎo)全張量重力梯度儀的加速度計(jì)輸出模型[3,7-9]。 全張量重力梯度儀由3 個(gè)GGI 組成, 安裝在由陀螺儀穩(wěn)定的平臺(tái)上, 3 個(gè)GGI 的軸線互相垂直。 為了簡(jiǎn)化計(jì)算, 以3 個(gè)圓盤的正交旋轉(zhuǎn)軸建立平臺(tái)坐標(biāo)系OXYZ, 再以旋轉(zhuǎn)軸為法線在每個(gè)圓盤上建立測(cè)量坐標(biāo)系OiXiYiZi(i=1,2,3), 如圖4 所示。 根據(jù)3 個(gè)圓盤之間的位置關(guān)系, 已知1 個(gè)質(zhì)量塊在1 個(gè)圓盤坐標(biāo)系位置的情況下可推導(dǎo)出相同的質(zhì)量塊與另1 個(gè)圓盤的位置關(guān)系。 把質(zhì)量塊看作一個(gè)質(zhì)點(diǎn)P, 已知質(zhì)點(diǎn)P在平臺(tái)坐標(biāo)系下的坐標(biāo)為(x,y,z), 平臺(tái)坐標(biāo)系的重力加速度為g=(gx,gy,gz), 則將平臺(tái)坐標(biāo)系轉(zhuǎn)換到每個(gè)GGI 的測(cè)量坐標(biāo)系, 圓盤1 中心的重力加速度為g1=(g1x,g1y,g1z), 圓盤2 中心的重力加速度為g2=(g2x,g2y,g2z), 圓盤3 中心的重力加速度為g3=(g3x,g3y,g3z)。 因此可得到每個(gè)GGI 上的加速度計(jì)輸出。

3 信號(hào)調(diào)理

圖4 坐標(biāo)系示意圖Fig.4 Coordinate gravity gradiometer simulation system block diagram

信號(hào)調(diào)理的目的是將4 個(gè)加速度計(jì)的測(cè)量信號(hào)相加減, 消除了動(dòng)基座引起的共模角加速度和線速度對(duì)重力梯度儀輸出的影響并放大了梯度信號(hào)[4]。 該部分功能通過硬件電路實(shí)現(xiàn), 對(duì)重力梯度儀進(jìn)行半實(shí)物仿真。 半實(shí)物仿真將某些非線性度較高的關(guān)鍵部件實(shí)物引入仿真回路, 避免全數(shù)字仿真中由于建立非線性部件數(shù)學(xué)模型的誤差, 從而可以提高仿真的可信度。

為使單個(gè)圓盤的4 個(gè)加速度計(jì)信號(hào)同時(shí)輸出, 選用DAC7617 4 路串行輸入、 12 位電壓輸出的數(shù)模轉(zhuǎn)換器將計(jì)算機(jī)產(chǎn)生的加速度計(jì)數(shù)據(jù)轉(zhuǎn)換為加速度計(jì)信號(hào), 電路如圖5 所示。 為了提高信噪比, 選用了LT1814 4 運(yùn)放芯片搭建求和求差放大電路, 電路圖如圖6 所示。 4 路信號(hào)經(jīng)過求和求差放大電路后,利用單片機(jī)STM32F103ZET6 內(nèi)部集成的AD 通道對(duì)調(diào)理后的3 路信號(hào)(E1+E2),(E3+E4),(E1+E1)-(E3+E4)進(jìn)行采集[10], 并通過串口發(fā)送到計(jì)算機(jī)中進(jìn)行后續(xù)的數(shù)據(jù)解調(diào)處理。 串口電路如圖7 所示。

圖6 信號(hào)調(diào)理電路Fig.6 Signal conditioning circuit

圖5 DAC7617 驅(qū)動(dòng)電路Fig.5 The driver circuit of DAC7617

圖7 串口驅(qū)動(dòng)電路Fig.7 Serial port driver circuit

4 信號(hào)解調(diào)

信號(hào)解調(diào)[11-16]與反饋模塊由計(jì)算機(jī)B 完成, 該部分主要包括濾波、 數(shù)據(jù)解調(diào)以及加速度計(jì)標(biāo)度因數(shù)不一致信息反饋。 整體流程圖如圖8 所示, 主要過程為: 下位機(jī)通過串口向計(jì)算機(jī)B 傳輸數(shù)據(jù), 并將數(shù)據(jù)存儲(chǔ)在讀取緩沖區(qū)中, 不斷進(jìn)行讀取, 將數(shù)據(jù)存儲(chǔ)在數(shù)組中; 為了提高信號(hào)的信噪比和解調(diào)精度, 在解調(diào)前通過中值濾波器對(duì)數(shù)據(jù)進(jìn)行濾波; 信號(hào)解調(diào)在濾波后進(jìn)行, 該部分能將加速度計(jì)標(biāo)度因數(shù)不一致量及重力梯度值進(jìn)行解調(diào), 并利用PID(Proportion Integral Derivative)控制器將標(biāo)度因數(shù)不一致信息反饋給計(jì)算機(jī)A, 計(jì)算機(jī)A 對(duì)標(biāo)度因數(shù)調(diào)整后重新發(fā)送數(shù)據(jù), 計(jì)算機(jī)B 進(jìn)行處理, 直到解調(diào)出的重力梯度滿足精度為止。

圖8 信號(hào)解調(diào)框圖Fig.8 Signal demodulation block diagram

由信號(hào)產(chǎn)生部分可知, 當(dāng)加速度計(jì)標(biāo)度因數(shù)不一致時(shí), 以水平圓盤為例, 單個(gè)圓盤的輸出為

從式(5)可以看出, 在輸出信號(hào)的一倍頻處解調(diào)[3-7], 可得到加速度計(jì)1 和加速度計(jì)2 的標(biāo)度因數(shù)不一致量及加速度計(jì)3 和加速度計(jì)4 的標(biāo)度因數(shù)不一致量, 通過串口, 將標(biāo)度因數(shù)不一致量反饋給計(jì)算機(jī)A,經(jīng)過調(diào)整使兩兩加速度計(jì)的標(biāo)度因數(shù)達(dá)到一致。 通過對(duì)平臺(tái)輸出信號(hào)進(jìn)行積分解調(diào), 即輸出信號(hào)與相乘后積分, 可得重力梯度值。

通過對(duì)兩個(gè)相對(duì)的加速度計(jì)的標(biāo)度因數(shù)實(shí)時(shí)反饋調(diào)節(jié), 可以有效的抑制平臺(tái)的平動(dòng)加速度對(duì)輸出信號(hào)的影響。 筆者利用增量式PID 控制器, 該控制器在進(jìn)行控制時(shí)需要3 個(gè)時(shí)刻的輸出結(jié)果。 通過設(shè)置PID 的3 個(gè)參數(shù)(Kp,Ki,Kd)調(diào)整加速度計(jì)標(biāo)度因數(shù), 這3 個(gè)參數(shù)決定了控制過程的過渡時(shí)間和穩(wěn)定性。經(jīng)過計(jì)算得到合適的PID 參數(shù)后, 將經(jīng)過PID 調(diào)節(jié)后的結(jié)果反饋給計(jì)算機(jī)A。

5 仿真結(jié)果

假設(shè)地下異常體為理想的球體, 以3 個(gè)圓盤的正交旋轉(zhuǎn)軸建立平臺(tái)坐標(biāo)系OXYZ, 參數(shù)設(shè)置如表1所示。 重力梯度儀的參數(shù)設(shè)置也如表1 所示。 根據(jù)加速度計(jì)的工作原理, 運(yùn)用LabVIEW 仿真全張量旋轉(zhuǎn)加速度計(jì)重力梯度儀產(chǎn)生12 個(gè)加速度計(jì)數(shù)據(jù), 并對(duì)數(shù)據(jù)進(jìn)行解調(diào), 得到全張量的重力梯度值。 正演仿真結(jié)果和解調(diào)結(jié)果如表2 所示。 在實(shí)驗(yàn)中, 多次對(duì)表1 中的參數(shù)設(shè)置進(jìn)行更改, 可以在表2 看到對(duì)解調(diào)數(shù)據(jù)精度的影響。 經(jīng)過多次測(cè)試, 發(fā)現(xiàn)圓盤半徑為1 m, 采樣頻率為16 Hz 時(shí), 數(shù)據(jù)解調(diào)的精度最高, 驗(yàn)證了筆者所使用的數(shù)據(jù)解調(diào)方法是有效可行的。

表1 實(shí)驗(yàn)參數(shù)設(shè)置Tab.1 Experimental parameter setting

表2 測(cè)試結(jié)果Tab.2 Test results

下面以單圓盤為例, 引入部分硬件電路, 對(duì)重力梯度儀進(jìn)行半實(shí)物仿真, 通過對(duì)Γxx-Γyy和Γxy進(jìn)行解調(diào), 進(jìn)一步驗(yàn)證該數(shù)據(jù)解調(diào)方法的可行性和可靠性。 圓盤的參數(shù)設(shè)置如表1 所示, 調(diào)整加速度計(jì)的標(biāo)度因數(shù), 通過硬件對(duì)信號(hào)進(jìn)行運(yùn)算, 再由計(jì)算機(jī)B 對(duì)信號(hào)進(jìn)行解調(diào)處理, 由于硬件電路引入部分隨機(jī)噪聲, 因此相當(dāng)于4 個(gè)加速度計(jì)的標(biāo)度因數(shù)存在不一致量。 實(shí)際上, 在每個(gè)GGI 上的加速度計(jì)的擺放都不可能保證完全的正交, 通過調(diào)整加速度計(jì)的標(biāo)度因數(shù)可以消除這種誤差。 在仿真過程中由計(jì)算機(jī)B 將標(biāo)度因數(shù)調(diào)整量和解調(diào)結(jié)果反饋給計(jì)算機(jī)A, 對(duì)輸出數(shù)據(jù)進(jìn)行調(diào)整, 使解調(diào)出的梯度值精度滿足要求。 受到硬件電路的限制, 所選用的ADC 和DAC 都是12 bit 且3.3 V 供電, 因此測(cè)量的重力梯度值的范圍在-5 000 E ～+5 000 E, 精度可達(dá)到100 E。

單圓盤的仿真結(jié)果如表3 所示。 首先, 設(shè)置采樣點(diǎn)數(shù)為64, 每隔125 ms 從計(jì)算機(jī)A 發(fā)送一組加速度計(jì)數(shù)據(jù)到下位機(jī)。 發(fā)送格式為“0num64N16fsaAbBcCd”, 其中, 0 為發(fā)送數(shù)據(jù)的序號(hào), 作為接收數(shù)據(jù)的校驗(yàn), 64 為采樣點(diǎn)數(shù), 16 為采樣頻率, “a”,“b”,“c”, “d” 分別表示4 個(gè)加速度計(jì)的輸出數(shù)據(jù)。 接收格式為 “ek1Eek2At1Bt2”, 其中 “ek1” 為加速度計(jì)1 的標(biāo)度因數(shù)調(diào)整量, “ek2” 為加速度計(jì)3 的標(biāo)度因數(shù)調(diào)整量, “t1” 為 Γxx-Γyy的解調(diào)結(jié)果, “t2” 為 Γxy的解調(diào)結(jié)果。 計(jì)算機(jī)A 將數(shù)據(jù)發(fā)送到硬件電路進(jìn)行處理后發(fā)送給計(jì)算機(jī)B, 計(jì)算機(jī)B 的接收格式為 “0num64N16fss1As2Bs3”, 其中 “s1” 和 “s2” 分別為圓盤上相對(duì)兩加速度計(jì)相加后的數(shù)據(jù), “s3” 為 “s1-s2” 后的數(shù)據(jù)。 計(jì)算機(jī)B 解調(diào)后將調(diào)整量和解調(diào)值反饋給計(jì)算機(jī)A, 發(fā)送格式與計(jì)算機(jī)A 的接受格式相同。

通過對(duì)表3 中的數(shù)據(jù)進(jìn)行分析可知, 在單圓盤重力梯度儀半實(shí)物仿真的精度在100 E 內(nèi)。 其誤差主要來源于以下幾個(gè)方面: 1) 數(shù)模轉(zhuǎn)換的精度低, 在轉(zhuǎn)換過程中丟失部分?jǐn)?shù)據(jù); 2)單片機(jī)內(nèi)部集成的ADC精度低且在0 ～0.5 V 之間存在非線性區(qū); 3) 信號(hào)經(jīng)過硬件電路后含有噪聲和工頻干擾。 因此, 在后續(xù)的研究工作中, 為了提高重力梯度儀的精度, 一方面可以選用精度更高的D/ A 和A/ D, 盡可能多地保留數(shù)據(jù)的精度位數(shù), 另一方面改進(jìn)硬件電路, 比如采用差分輸入差分輸出的電路形式抑制共模干擾, 電源處做好電磁屏蔽等減少信號(hào)通過硬件電路產(chǎn)生的噪聲。

表3 單圓盤半實(shí)物仿真測(cè)試結(jié)果Tab.3 Test results of hardware in the loop simulation of single disk

6 結(jié) 語

筆者依據(jù)旋轉(zhuǎn)加速度計(jì)重力梯度儀的測(cè)量原理, 對(duì)地下某異常體引起的重力梯度異常值進(jìn)行了計(jì)算, 簡(jiǎn)化梯度模型后, 得到了全張量的重力梯度。 通過模擬加速度計(jì)數(shù)據(jù), 對(duì)全張量重力梯度儀進(jìn)行了仿真、 運(yùn)算和數(shù)字解調(diào)。 以單圓盤為例, 對(duì)重力梯度儀進(jìn)行了半實(shí)物仿真, 研究了旋轉(zhuǎn)加速度計(jì)重力梯度儀的信號(hào)產(chǎn)生原理, 加速度計(jì)信號(hào)解調(diào)重力梯度值的方法以及反饋標(biāo)度因數(shù)不一致的調(diào)整方法, 為重力梯度儀硬件設(shè)計(jì)提供理論和技術(shù)支持。