空載重力測量系統架構及演算法之設計

在一個完整的空載重力測量的施測過程中，觀測資料共分為兩個種類，一為 GPS 觀測資料，另一為重力儀觀測資料。然而一個完整的空載重力測量系統也可 分為兩個部分，一為 GPS 動態定位系統，用於提供飛機的位置，另一為重力值約 化系統，由飛機的位置計算飛機的速度量與加速度量，並結合重力儀觀測資料，

求解空間中該點的重力值或重力異常值。在本章中將針對這兩個系統的各別解算 流程並整合兩個系統形成一個完整的空載重力測量系統，最後並說明空載重力測 量時之 GPS 作業注意事項與空載重力測量作業程序。

4-1 GPS 動態定位系統

相對於靜態定位而言，動態定位的意思就是希望能夠在待定點為運動狀態下 求解其位置，對空載重力測量而言，吾人必須知道每一個觀測時刻載具之位置，

進一步求解載具之速度及加速度，因此必須利用 GPS 採取動態定位的方式求解載 具的位置，其中最關鍵的部分，即在於快速的求解 GPS 載波相位觀測量中的整數 週波未定值，吾人選擇 OTF 週波未定值演算法做為 GPS 動態定位系統中求解整數 週波未定值的方法。GPS 動態定位系統架構則如圖 4-1 所示。

4-1-1 OTF 週波未定值演算法

在要求較高定位精度的情況下，就必須採用相對定位(Relative)的方法，同時 採精度較高的載波相位觀測量。相對定位是指移動站坐標的決定乃以相對於某一 已知坐標之主站的方式。至於載波相位觀測量因具有非常小的偶然誤差，非常適 合用於精確的定位，然而卻具有載波相位未定值(Phase Ambiguity or 相位模糊度) 的問題(楊，1997)。因此，載波相位未定值的求解速度及正確性將是影響動態定位 的關鍵。在此使用 OTF 演算法做為我們載波相位未定值求解之方法。

將四個二次差分方程式組成如下之矩陣形式：

17

ˆ 60 ^{1 2}

2

K

N K −

=

N^ˆ₁ =N^ˆ₂+K₁ (4-5)

由(4-5)式可知，如果能精確觀測K 、₁ K ，即誤差項為零，便可以輕易求得整₂ 數解N 、₁ N 。若誤差項不為零，則任一組整數₂ K 、₁ K 代入不一定可求得整數的₂ N 、1 N ，因此可允許₂ K 的誤差範圍在₁ ±2cycle、K 的誤差範圍在₂ ±6cycle，在 這麼多整數的K 、₁ K 的組合當中，可求得僅有一組₂ K 、₁ K 可求解出整數的₂ N 、₁

N 。 2

建立搜尋空間快速求得所有衛星的相位未定值，再代入所有衛星觀測量進行 最小二乘平差，求得相位觀測量的改正值和後驗變方因子，然後選擇最大和最小 的的後驗變方因子進行統計的驗證來決定正確的相位未定值(楊，1997)。

經由以上演算法計算後，我們便能快速求解得主站和移動站的 GPS 天線間的 基線向量，進而快速求得飛機的位置向量，在此我們可得到飛機的緯度量，經度 量與橢球高。

圖 4-1 GPS 動態定位系統架構圖 No

移動站觀測資料 (xxxxx.0xo)

計算 GPS 衛星軌道

起始值的獲得 (利用電碼解) 建立搜尋空間

驗證與剔除準則

比率測試與合理 性測試

確定相位未定值

載具之 WGS84 坐標

OTF 演算法 地面站觀測資料

(xxxxx.0xo)

廣播星曆(xxxxx.0xn)或 精密星曆(xxxxx.sp3)

No

Yes

Yes

4-2 重力儀觀測值解算流程

f 為重力儀於空間中該點的觀測值 _Z

0 Eotvos local local

z

(1) 位置修正

利用程式語言 Fortran 的 IMSL 資料庫所提供的微分副程式(DERIV/DDERIV) 加以撰寫改編，使之可利用於我們所需的用途。此微分程式的原理為將資料組成

(3) 坐標轉換

利用由前一個步驟所得之飛機之三維速度量(x&^,y&^,z&^{)、加速度量(}x&&^,y&&^,z&&^)，透 過坐標系統轉換的方式轉換至飛機的區域坐標系下，即可得在飛機區域坐標系下 沿 n,e,u 方向的速度量、加速度量。在我們將飛機之三維速度量(x&^,y&^,z&^{)、加速度}

量(x&&^,y&&^,z&&)經過坐標轉換方式得到飛機區域坐標系下之三維速度(n&^,e&^,u&^)及三維

加速度(n&&^,e&&^,u&&)，其轉換的過程意涵等同於從地固坐標系統(Conventional Terrestrial

System，CTS)轉換至區域坐標系統的過程。因此我們可以得到如下之轉換公式：

在求解得到由地固坐標系統(CTS, Conventional Terrestrial System)轉換至地區 坐標系統的轉換矩陣後，我們將飛機之三維速度量(x&^,y&^,z&^{)、加速度量(}x&&^,y&&^,z&&⁾

( ) ∑ ( )

∑

圖 4-2 重力值約化系統架構圖 天線位置

t ，a φ_a，λ_a，h _a

輔助檔案

d，S ，_d h， f ，_z0 g 0

重力儀位置(φ,λ,h) t ，a φ ，_g λ ，_g h_g

重力儀位置(X,Y,Z) t ，a X_g，Y_g，Z_g

在 (X,Y,Z)坐標系統下的 速度量與加速度量

t ，a X& ，_g Y& ，_g Z& ， _g

X&& ，g Y&& ，_g Z&&_g 位置修正

坐標轉換和內插

微分處理

在(N,E,U)坐標系統下的速 度量、加速度量與科氏力

t ，a φ ，_g λ ，_g h_g，_n&g，_e&g， u&g，n&&_g，e&&_g，u&&_g，δg_Eto&&vo&&s 坐標轉換、濾波處

理和計算科氏力 改正項

重力儀觀測資料 tg， f _z

Time shift 修正後 的重力儀觀測資料

t′g， f _z

(

)

由u&&_g與 f 的關係_z 找出 time shift

空中的重力值

( )

f

g = _z + − &&_g +δ _Eto&&vo&&s− _z +

圖 4-3 模擬飛機飛行時之側視圖

4-4 GPS 作業注意事項

動態 GPS 定位測量之程序不同於一般 GPS 靜態測量之程序，因此探討動態 GPS 定位測量之主站位置與移動站設置的方式，以做為日後實際作業之參考。

(1)主站位置之選擇

主站必須是已知坐標的控制點，周圍環境必須對空通視良好，及仰角 10 度以 上無遮蔽物。此次台灣地區空載重力測量所選擇的主站之位置如表 4-3 所示。

(2)飛機 GPS 天線位置之選擇

飛機上 GPS 移動站之天線必須裝置在對空通視良好之處，而我們使用的飛機 是將天線安裝在機背上，所以並無通視上的問題。天線的安裝必須在飛機上挖洞，

由於考量飛機結構安全問題，因此天線的位置必須經過飛機結構工程師確認沒有 問題之後，才能在該位置進行安裝。由於我們所使用的飛機僅安裝了一個天線，

因此安裝分訊器(splitter)將天線接收訊號分離至兩部 GPS 接收儀，為用來檢驗定位 精度的一可行方法且能防止儀器損壞造成無觀測資料。

4-5 空載重力作業程序

完整的空載重力測量作業程序如下所列之步驟：

1. 於機場建立一重力基準點，此重力基準點可由最近之重力點引測。

2. 以相對重力儀(如 GRAVITON-EG)觀測重力基準點與飛機上 LCR 空載重力 儀之間的相對重力，此步驟僅需施測一次，目的為求得飛機停放位置與重 力基準點間的相對重力。

3. 飛機起飛後，觀測飛行經過地區與起飛點之相對重力。

4. 降落後，在飛機停放的位置上，此時並不把空載重力儀關機，而是使其繼 續觀測達 30 分鐘以上，以求得其最穩定的讀數。

5. 以 GPS 資料後處理計算飛機之位置、速度、加速度而改正步驟 3 測得之相 對重力。

6. 以步驟 3,4,5 求得之空中的重力值與現有重力(地面重力經向上延續後)比較 以偵測粗差。

7. 依步驟 3∼6 完成所有航線，再進行交叉點平差分析。

8. 向下延續重力到適當的地面高度以為後續之應用。

表 4-3 GPS 動態定位主站列表

名稱 編號 地心地固坐標系統下之坐標

水湳機場站 7803 (-2968166.2744,5008620.4250,2596236.7407)

陽明山站 ymsm (-3024807.5842,4921747.1613,2696033.6580)

北港站 Pkgm (-2951343.5090,5049506.4590,2535725.4330)

墾丁站 kdnm (-3028999.6338,5084820.6365,2369241.4009)

鳳林站 flnm (-3048091.3511,4983134.4170,2552645.0155)

太麻里站 tman (-3034565.1975,5048871.0503,2437550.6962)