高解析度衛星影像之定位精度研究 (未完)

高解析度衛星影像之定位精度研究 (林義乾 ,指導教授：趙建哲老師)

一、 研究動機：
近年來遙測技術的發展快速，空間解析度從稍早以前SPOT1~4全色態影像的10米空間解析度，至今QUCIKBIRD已可達到次米級（全色態0.61~0.72m）之高解析度，雖然相較於航空攝影測量之空間解析度仍有一段差距，但衛星影像具有取得快速，且衛星行徑有一定的規律，並可由地面基地台直接控制衛星感測器掃描範圍，不需傳統航測之航線規劃過程，對於同樣大小之區域，衛星感測器之姿態變化相對於空載像機之變化為小，此外高解析度衛星影像皆採用推掃式（push-broom）掃描，僅於y方向會產生幾何變形。縱上所論，高解析度衛星影像（High Resolution Satellite Imagery, HRSI）具有較穩定和效率較高等優勢，對於大比例尺製圖方面的應用理論上較傳統航測便利，因此透過分析HRSI的定位精度，即可探討其可用於何種比例尺的製圖。
二、 研究目的：
根據目前已有之方法，分析各高解析度衛星可達之最高定位精度。並預期能找出提升定位精度之方法。

三、 高解析度衛星介紹：
近年來，許多國家陸續發射高解析度衛星，自1999年9月24日IKONOS商用高解析度衛星發射升空之後，高解析度衛星影像（High Resolution Satellite Imagery, HRSI）邁入一米解像力的紀元，此後2000年QUICKBIRDⅡ更達到次米級之解析度，而今年我國也於5月21日發射了華衛二號，號稱可達兩米之解析度（全色態影像），另外較早發射的SPOT 5、EROS等，也都有2.5米以上之解析度（SPOT 5於super-mode下可達2.5米解析度，EROS在兩米左右），由於目前對HRSI並無清楚的定義，且因各種應用而會對要求的解像力有所不同，本研究著重於定位精度，因此選定上述幾顆衛星加以探討，以下為各衛星的基本資料及用於定位研究所需之相關資料介紹。

1. 高解析度衛星基本資料：
衛星 SPOT 5
（HRG） ROCSATⅡ EROS A1 IKONOS QUICKBIRD
焦距（m) 1.082 2.898 3.4 10 8.839
衛星高度（km） 822 891 480 681 450
成像方式
推掃式(push-broom) 推掃式(push-broom) 推掃式(push-broom) 推掃式(push-broom) 推掃式(push-broom)
掃描模式 同步取樣 同步取樣 非同步取樣 同步取樣 同步取樣
獲得立體像對方式 同軌道或跨軌道 跨軌道 同軌道或跨軌道 同軌道或跨軌道 同軌道或跨軌道
感測器間最大夾角/（deg） 同軌道：27  同軌道：45 同軌道：26
跨軌道：45 同軌道：25
掃描寬度 2×60 km 24 km  12.5 km  11 km 16.5 km
空間解析度（m） XS：10
PAN：2.5 or 5 XS：8
PAN：2 PAN：1.8 XS：4
PAN：1 XS：2.44~2.88
PAN：0.61~0.72
輻射解析力 PAN：0.49~0.69
B1：0.50~0.59
B2：0.61~0.68
B3：0.78~0.89
SWIR：1.58~1.75 PAN：0.45~0.90
XS： 0.45~0.52
    0.52~0.60
    0.63~0.69
    0.76~0.90 PAN：0.5~0.9 PAN：0.45~0.90
XS： 0.45~0.52
    0.52~0.60
    0.63~0.69
    0.76~0.90 PAN：0.45~0.90
XS： 0.45~0.52
    0.52~0.60
    0.63~0.69
    0.76~0.90
資料型態 8 Bits  11 Bits 11 Bits 11 Bits
衛星軌道種類 太陽同步衛星 太陽同步衛星 太陽同步衛星 太陽同步衛星 太陽同步衛星
衛星軌道與赤道交角 98.7度  97.3度 98.1度 98度
衛星軌道參數 提供 預計不提供 提供 目前不提供 提供
衛星軌道週期 101.4 mins
26 days  15.3 圈/天 98 mins
14.6圈/天 93.5 mins
飛行速度 7.4 km/sec   6.79 km/sec
表1－高解析度衛星基本資料比較表
其中EROS A1之掃描模式為非同步取樣（Asynchronous），其掃描速率約為750 lines/sec，而同步取樣與非同步取樣之差異如下圖：

圖1－同步取樣與非同步取樣
另外IKONOS衛星採用三線式立體成像系統，擁有3排CCD陣列感測器，分別可向前（fore）、向下（nadir）、向後（aft）掃瞄成像。當衛星在軌道上行進時，如下圖2，前視，底視及後視之掃描器分別對地面掃描成像，由此類似航帶攝影之方式，可獲得前後重疊之影像已組成立體模型，以提供後續之應用[李茂園,2001]。

圖2－三線式影像成像方法

圖3－三線式掃描器立體成像原理及座標系統定義
2. 各衛星之影像處理等級：
  SPOT 5：以中央大學太空遙測中心公佈之資料為主。
處理等級 說明
原始影像(等級1A、1B-相對於原等級1、4) 『原始影像』是指將最原始所接收到的資料進行輻射方面之相對及絕對校正，其幾何部份則保持原始掃瞄紀錄之形式。由於衛星感測器通常以一定數量之電荷耦合裝置(Charged Coupled Device, CCD)來拍攝地表，各個CCD之間的輻射反應函數不盡相同，因此需要進行相對的輻射校正，未經此校正之影像在均調區域將呈現不同灰階之線條(Striping)現象。此外，所需之校正，則是為了修正衛星感測器在不斷的運轉下，隨著運作時間的增加而降低其敏感度所產生所謂的劣化(Degrade)現象。進行整體的輻射校正後，方能使所獲取之影像維持固定的輻射品質。
系統改正影像(等級2A-相對於原等級8) 『系統改正影像』之輸出範圍位置大小，可依使用者需求選取範圍。此類產品在幾何糾正上，利用軌道載體所紀錄之衛星位置與姿態等參數來預估影像四個角點之地面座標，再進行影像糾正，同時將成果投影在所需要之地圖座標系統上。不同的衛星，因為所用以探測及紀錄載體參數之儀器精度不同，使得此系統改正影像之幾何 絕對精度有所不同。
精密幾何改正影像(等級2B、3-相對於原等級9、10) 在製作『精密幾何改正影像』時，首先利用本站所建立之全台灣地面控制點資料庫來點選相對應之影像控制點，求解得到精確之方位。再利用林務局農航所製作之全台灣40 公尺網格數值地形模型(Digital Terrain Model, DTM)資料進行位移修正，並將成果投影在設定之地圖座標系統上成為等級3產品，以產生正射影像(Ortho-Image)。若為等級2B產品則未使用數 值地形模型，而改以控制點平均高度來進行幾何糾正。此產品之影像範圍及位置，可依使用者需求給定，不限定為全幅影像之輸出範圍。針對台灣及澎湖地區所提供之地圖投影系統可以為UTM、TWD67及TWD97三種，但上述等級2B及等級3之標準產品以TWD67之地圖投影系統為主，若有使用者需要UTM或 TWD97地圖投影系統之產品，必須另外指定。
表2－SPOT 5衛星影像各處理等級及說明[太遙中心使用者手冊,2003]
  EROS A1：中央大學太空遙測中心亦有接收，故也以其公佈資料為準。
處理等級 改正項目
原始影像－等級1A 影像範圍約13.5公里x 13.5公里，影像空間解析度隨拍攝角度而改變。在近地點之空間解析度為1.8公尺。影像格式為Raw Format，輻射解析度為11 Bits/Pixel，儲存時為16 Bits/Pixel。
系統改正影像－等級1B 此等級影像與SPOT等級2A相同，影像已經過系統幾何改正，影像上方為地理北方。影像格式為GeoTiff，輻射解析度為11 Bits/Pixel，儲存時為16 Bits/Pixel，採用WGS84-UTM座標系統。
精密幾何改正影像－等級3 僅限台灣地區或使用者提供地面控制點與數值地形模型之區域。最大範圍與原始影像範圍一致，台灣地區採用TWD67-TM2座標系統。影像格式為Raw Format，輻射解析度為11 Bits/Pixel，儲存時為16 Bits/Pixel。
表3－EROS A1衛星影像各處理等級及說明[太遙中心使用者手冊,2003]
  IKONOS：
CARTERRA Geo CARTERRA Reference CARTERRA Precision and CARTERRA Precision plus
已經過標準幾何改正之影像，以軌道參數及投影方法將影像糾正至參考橢球上，約有50公尺的誤差（CE90）。
此類影像已經過幾合糾正及地圖投影過程，但未加入地面資料故不能稱之為正射糾正。
嚴格地說，此影像已無法滿足製圖應用需求中之最原始影像幾何，及不適合進行光束法平差分析。 以軌道參數及投影方法將影像糾正至橢球上，但加上DEM資料，因此可使誤差縮小至25公尺（CE90），在攝像仰角小於60度內皆可接受。 利用DTM資料、控制點資料及軌道參數來糾正影像，即一般所謂之正射影像，平面精度最高可達到2公尺（CE90），而高程可達到3公尺（CE90），但美國本土以外影像使用者需自行提供DTM資料及控制點資料，或需政府單位才能購得。
表4－IKONOS衛星影像各處理等級及說明[李茂園,2001]
  QUICKBIRD：
Basic處理等級 Standard處理等級 Ortho Ready Standard
處理等級
影像處理過程，包含：
輻射校正
衛星系統校正
空間解析度
全色態影像：61~72公分
多光譜影像：2.44~2.88公尺
資料型態（Dynamic Range）8 Bits 或 16Bits 影像處理過程，包含：
輻射校正（radiometric distortion）
衛星系統校正（sensor distortion）
地圖投影（map projection）
空間解析度
全色態影像：70公分
多光譜影像：2.8公尺
資料型態（Dynamic Range）8 Bits或16Bits 影像處理過程，包含：
輻射校正（radiometric distortion）
衛星系統校正（sensor distortion）
地圖投影（map projection）
空間解析度
全色態影像：70公分
多光譜影像：2.8公尺
資料型態（Dynamic Range）8 Bits或16Bits
表5－QUCIKBIRD衛星影像各處理等級及說明[RITI,2004]
四、 方法介紹：
1. 高解析度衛星成像幾何：
上述介紹的高解析度衛星影像皆採用推掃式（Push-broom）之掃描方式，因推掃式較橫掃式（Whisk-broom）掃描具有較好的空間與解析力，更好的幾何性（偵檢器元件間有固定的關係），較輕小的設備，需求較少能量等優點[黃灝雄,2004]。
由於線性陣列（linear arrays）較區域陣列（area arrays）容易建構且不須機械掃描程式，通常被用在固定姿態的感測器上。推掃式是沿著感測器前進之方向掃描，所得之影像會與航向相垂直，且每一列（line）影像由於載台的移動會有不同的位置和方位，使得推掃式影像之幾何強度受限，在傳統航測應用上，通常採用GPS及INS來改善上述情形，載台的移動速度決定感測器的取樣速率，但若感測器瞄準邊緣時，對地速度的改變會造成取樣過度或是不足的情形[Edward,2001]。
高解析度衛星影像皆是採推掃式掃描，因此一張完整的衛星影像乃由數千至數萬"列"影像所組成，理論上每一列掃描線皆有其對應之外方位參數，若吾人欲探討HRSI之定位問題，應逐一解算各列影像之外方位參數，但在實務上無法取得如此大量的控制資料來解算全部的外方位元參數（lines×6）。此外亦有部分衛星之軌道參數不公開，更使得利用共線式（Collinearity equations）來進行影像的幾何改正有諸多的困難，因此部分學者提出了不考慮影像的物理意義的幾何校正模式：多項式糾正（Polynomial Rectification）、透視投影轉換（Projective Transformation）、有理函數轉換（Rational Function Model ,RMF）等等，數學模式較簡單且計算速度較快，不需影像之方位參數，因此較不嚴密；但若有衛星的軌道參數，則可使用較嚴謹的物理校正模式進行幾何改正，可考慮成像時造成成像變形的原因，並利用共線式配合DTM資料來進行影像糾正，糾正效果較佳[李茂園,2001]。
2. HRSI之成像變形因素：
Guennadi[2004]提到衛星影像在成像時，會受到感測器（視角變化）、衛星載台（姿態、高度、速度）、地球本體（自轉、地球曲率）和拍攝時的情形（大氣折光、高差移位）等產生幾何變形，由於產生這些變形之原因已知，且目前對其影響程度已有相當的瞭解，故上述的系統性變形是可被改正的。
  視角傾斜：
當感測器的像主軸（principle axis）或光軸（optical axis）傾斜時，會使影像幾何產生明顯的變形，此一情形可以數學模式描述並完善的改正至垂直面上，但對掃描地表的認知不足可能會使改正存有殘差，乃因地形高度與傾角有關。
另一個因視角傾斜會產生的現象，亦即一個像元相對的涵蓋範圍會隨著傾角越大而增加，或可說是對地有效像元變大，引入瞬間視場角（IFOV）的概念，當感測器在nadir方向時其像元大小為pix，傾斜α角後為

其中H為載台高度。
  地球曲率：

圖3－地球曲率所產生的移位
如上圖3，地球表面點P投影至卡氏座標系（Cartesian coordinate system）上會落於在P"點，其在影像上及切平面上的相對移位分別為 和 ，當載台高度H、地球半徑R、焦距f和垂直影像上點位元與像主點間之距離 ，則可計算 ，在利用比例關係計算 ，如下

若要採取更精確的改正，可由地球為一橢球體之概念來計算

計算D可由光線SP與地球表面的交點P之位置決定

其中 ，並以橢球座標方程式表示P點座標如下

其中 和 分別表示P點之經度及緯度，e為地球之離心率，N為曲率半徑，則X可表示為

代入前述公式及可求得位移量 。
  地球自轉：
在掃描的過程中，地球同時呈現西向東的自轉，會使取得之影像像元產生平移（shift）的現象，因此必須將影像線列向西方等比例的補償校正以消除此一位移。根據上圖3，在P點上之地球自轉速度為 ，其中 為平均地球自轉角速度。
而衛星之角速度為 ，其中 。將衛星之速度向量 轉換至地表，表示為

為地心座標系之經度，而 乃由 及 所組成，衛星軌道之方位角κ，可由軌道傾角θ計算，其關係為 。
因此地球自轉與衛星之相對速度（線性）即為

  大氣折光：
從太空中感測地表資料，光波在通過大氣層時會產生折射現象，會影響視角及點位的真實位置，而HRSI乃由太空中掃描，其大氣折光校正方式與空載航測像片不同，必須經由以天文學為基礎之校正方式。
  衛星參數差異：
建構一個精密模式來描述微變的衛星運動是非常困難的工作，乃基於精密的太空力學演算法和精確的軌道元件，衛星的高度、姿態與速度參數的變化皆會使影像產生相當程度的位移。
3. HRSI影像校正方式：
  仿射轉換（Affine Transformation）：
仿射轉換屬於幾何校正模式中多項式轉換中的低階轉換，由於未知參數較少，在解算上最簡單且迅速。首先探討平面仿射轉換，由於衛星高度遠大於地表起伏，故將HRSI與相對應之地球表面視為兩個平面，不考慮點位高程變化所產生的高差移位，其考慮到兩座標系之間的平移、旋轉、兩軸尺度不一致及兩軸不正交等因素，但根據李茂園[2001]指出平面仿射轉換因未考慮地形起伏所導致的高差移位影響，故僅適用於平坦地區之影像糾正。
引入地面點位元的高程值，並將影像座標以 的三維模式呈現（由於影像座標系統為一平面，高程值自然為0），即可進行三維仿射轉換，其考慮座標系間原點的平移、三個方向不一致的尺度變化及三軸的旋轉量。
其數學模式

以矩陣方式表示為，經過簡化之後可表示成

三維仿射轉換共有八個未知參數，理論上只需量測四個點即可求解，但為了消除量測的偶然誤差，會量測四個以上的點做為多餘觀測，並利用最小二乘法進行解算，以求得最佳之轉換參數。

圖4－中心透視投影與仿射投影[C.S. Fraser, 2004]
如上圖所示，仿射投影從像空間投影至物空間時，並非投影於水平面上，而是基於平行投影（parallel projection）的概念，因此兩平面間有一傾角α，Zhang and Zhang（2002）提出嚴僅仿射幾何轉換模型（strict affine geometric model），可校正因傾角所造成的偏差，又由於其為平行投影，故可將物空間的傾斜情形移致像空間，此傾角會使像點之y座標產生偏移（因其為推掃式，僅有y方向有變形）如圖5，必須進行透視投影與仿射投影之間的轉換，因此前述之三維仿射轉換應修正為

其中
圖5－物點P之影像座標偏移[C.S. Fraser, 2004]
此修正模式可應用於涵蓋大範圍區域的衛星影像及地形起伏較大的山地區域，且維持未知參數為八個，在解算上仍保有簡單迅速的特性。但三維仿射轉換（嚴謹模式亦同）仍將物空間視為一個三維平面，雖然有學者指出在觀測量足夠的情形下，其定位精度與有理函數模式相當接近[羅秋月,2002]，不過其模式明顯忽略地形起伏變化，對於仿射轉換之定位精度將在往後以實驗驗證。
  有理函數模式（Rational Function Model, RFM）：
C.V. Tao及T. Thierry等人於2000年提出以有理函數模式來校正IKONOS衛星影像[羅秋月,2002]。相較於一般多項式，RFM位較複雜的形式，因其考慮到高程的變化，故理論上較一般多項式更能有效表達地形的變化，且可應用在各種形式的感測器上，如框幅式（frame）、推掃式（push-broom）、橫掃式（Whiskbroom）和SAR影像[李茂園,2001]。
屬於幾何模式的校正，使用者不需知道感測器的相關資訊，在處理HRSI時，可用於衛星軌道資訊無法得知或是影像已經過處理的情形。
RFM之影像幾何模式是利用2個到4個多項式的比例關係來計算影像的行、列座標，每個多項式皆是地面座標 的函數，最多共有20項（3階RFM函數），共有78個未知係數，如下式：

其中r, c是正常化（Normalized）後的影像座標， 為正常化後的地面座標， 為多項式的係數，正常化是指利用平移（Offset）及尺度（Scaled）調整，使座標值的範圍在-1.0到+1.0之間。由C.V. Tao[2000]等人的研究指出，函數的一階項可消除光學投影產生的扭曲誤差，二階項可消除因地球曲率、大氣折射及透鏡扭曲等之變形，而其他未知因素所產生的變形則可由三階項來吸收。
當使用三階的RFM時，共要求解78個係數，至少需要39個控制點才能求解。利用最小二乘法（Least Squares Method）解算時，由於上式非線性方程式，必須先將其線性化，可先將上式以下面型式表示

其誤差方程式如下式中
 
表達成矩陣形式如下

其中
  令則 組成法方程式 未知數向量I可由下是求得

其中
  光束法平差：
當考慮到成像時的物理因素，如成像方法、方位參數等資訊，其糾正精度理論上應較幾何模式高。前面所介紹高解析度衛星皆採用推掃式掃描，因此每一條掃描線即屬於一透視投影，亦即每條掃描線皆有一組外方位參數（位置參數及姿態參數），而HRSI乃由數千數萬條掃描線所組成，無法利用大量的控制資料來直接求解如此大量的未知參數，必須予以簡化，一般是採用時間參數來描述衛星的方位參數，位置參數採用一階及二階表示，姿態參數以三階多項式表示。
根據前述的推掃式成像原理，在飛行方向上為平行投影，垂直飛行方向上為透視投影，與一般航測像片的純透視投影原理不同[李茂園,2001]，故將共線方程式略作修改如下

其中
   焦距
各掃描線相應影像中心線之掃描時間
附加尺度參數
第i點之影像座標
第i點之地面座標
投影中心之地面座標（為時間t之函數）

旋轉矩陣元素（為時間t之函數）

其中投影中心之姿態參數（為時間t之函數）：

A. 量測影像座標可得觀測方程式
 
再對各參數偏微分以線性化，則得誤差方程式，將未知參數分成地面座標及其他參數（包含外方位及附加尺度參數）兩群，以矩陣形式表示如下

簡化表示如下
 
其中 表示外方位參數及附加尺度參數 表示地面座標
B. 外方位參數（多項式係數）虛擬觀測方程式

其中參數符號右上方加"0"者代表該參數的近似值；而加"00"者代表該參數之觀測值，簡化如下：

C. 地面控制點虛擬觀測方程式

簡化如下

假設共有m張影像，每張影像有n點，則共可列得三類觀測方程式如下：

簡寫成矩陣表示式或 其相應的權矩陣為

其中 ：影像座標量測之權
  ：外方位參數之權
  ：地面點座標之權
故法方程式（Normal Equation）為

或簡寫為
式中

4. 後續工作－影像重新取樣流程：

圖6－影像重新取樣之流程[李茂園,2001]
五、 結論與未來方向：
1. 分析各成像變形因素對影像幾何的影響及探討正確改正方式。
2. 實際操作實驗以驗證幾何校正模式之正確性，並評估其定位精度。
3. 探討推掃式（push-broom）影像之幾何特性，分析是否有影像特性可加入平差模式作為約制條件。
4. 探討太陽同步（sun-synchronous）衛星具有的軌道特性，分析是否可作為平差之約制條件。
5. 根據太陽同步衛星之軌道特性，探討是否有優於時間參數的衛星位置參數與姿態參數的描述方式。
6. 探討影像幾何校正後之各項後續工作及其對定位精度的影響，如影像重新取樣方式。

六、 參考文獻：
  Ahmed Shaker, "The Line Based Transformation Model (LBTM): A New Approach to The Rectification of High-resolution Satellite Imagery"。
  Ball Company," http://www.ballaerospace.com/quickbird.html"。
  C.S. Fraser*, T. Yamakawa ,2003,"Insights into the affine model for high-resolution satellite sensor orientation", ISPRS Journal of Photogrammetry & Remote Sensing 58(2004) 275-288。
  Edward M. Mikhail,2001,"Introduction To Modern Photogrammetry" 。
  Guennadi A. Guienko,2004,"Geometric Accuracy of Ikonos: Zoom In", IEEE TRANSACTION ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO.1,209-214。
  Liang-Chien CHEN and Tee-Ann TEO,"RIGOROUS GENERATION OF DIGITAL ORTHOPHOTOS FROM EROS A HIGH RESOLUTION SATELLITE IMAGES"。
  RITI Technology,"http://www.quickbird.com.tw/"。
  SPOT 5," http://spot5.cnes.fr/gb/index2.htm"。
  李茂園,2001,"高解析度衛星影像之幾何處理與定位精度分析",國立台灣大學土木工程學研究所碩士論文。
  羅秋月,2002,"IKONOS衛星影像正射改正之研究",國立中央大學土木工程研究所碩士論文。
  黃灝雄,2004,"遙感探測課程投影片"。
  國立中央大學太空遙測中心," http://www.csrsr.ncu.edu.tw/"。
國家太空計畫室，" http://www.nspo.org.tw/c60/home/"。