重力梯度儀仿真系統(tǒng)設計

劉 杰, 閆 妍, 申夢嫻, 劉乾坤, 羅偉東

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

0 引 言

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

1 設計方案

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

圖1 旋轉加速度計重力梯度儀GGI 示意圖 Fig.1 Rotary accelerometer gravity gradiometer GGI schematic

圖2 全張量重力梯度儀結構示意圖Fig.2 Schematic diagram of full tensor gravity gradiometer

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

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

2 信號產(chǎn)生

2.1 重力梯度測量原理

空間質量受到的重力是由地球引力和地球自轉引起的慣性離心力構成, 重力位函數(shù)是一個標量函數(shù), 它對各個方向的偏導數(shù)等于重力在對應坐標軸上的分量。 重力位函數(shù)為

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

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

對重力位在各個方向上求二階導數(shù)可得到重力梯度

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

2.2 加速度計測量原理

如果忽略地球自轉角速度重力梯度儀平臺相對于地球坐標系的轉動角速度和圓盤自轉角速度, 根據(jù)擺式加速度計輸出的數(shù)學模型[5-8]可導出單個加速度計的輸出結果為

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

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

3 信號調(diào)理

圖4 坐標系示意圖Fig.4 Coordinate gravity gradiometer simulation system block diagram

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

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

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

圖5 DAC7617 驅動電路Fig.5 The driver circuit of DAC7617

圖7 串口驅動電路Fig.7 Serial port driver circuit

4 信號解調(diào)

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

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

由信號產(chǎn)生部分可知, 當加速度計標度因數(shù)不一致時, 以水平圓盤為例, 單個圓盤的輸出為

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

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

5 仿真結果

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

表1 實驗參數(shù)設置Tab.1 Experimental parameter setting

表2 測試結果Tab.2 Test results

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

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

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

表3 單圓盤半實物仿真測試結果Tab.3 Test results of hardware in the loop simulation of single disk

6 結 語

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