以光達資料之控制直線求解單張像片外方位參數之模式探討與可行性評估

(1)

Volume 10, 1, March 2005, pp. 89-102

以光達資料之控制直線求解單張像片外方位參數

之模式探討與可行性評估

趙鍵哲

1

彭念豪

2

摘要

空載光達(Airborne LIDAR)及地面光達(Ground-based LIDAR)系統能迅速地收集空間三維資料，提供豐

富的物空間幾何訊息；而像片資料則具有較豐富的輻射資訊，便於影像判讀，但其空間三維幾何必須藉助方位求解及交會定點才能確定。因此結合光達資料與攝影測量方法應用之一即為如何藉助從光達資料中萃取的物空間特徵(如點、線、面特徵)作為像片求解外方位參數之控制資料。本研究利用三維的光達資料，從中萃取物空間直線特徵(本文稱為控制直線)作為控制資料，利用直線特徵在物、像空間與投影中心之幾何關聯求解像片外方位參數，並在數學模式中引入控制直線之誤差估計量，藉以合理化平差程序獲致較佳的外方位參數解算成果。本研究中實例之分析除了說明光達控制直線用於求解單張像片外方位參數之實用性，並闡述平差計算後之誤差特性。

關鍵詞：控制直線、單張像片外方位參數求解

1. 前言

傳統攝影測量空中三角作業主要以量測點為主，例如連結點(Tie Points)及控制點(Control Points)，此類特徵必須在影像呈現足夠的訊號，例如屋角點及佈標點，但也可以是一些不具幾何意義的位置，例如屋頂的斑駁痕跡或地面的石塊角落。相對於此，線特徵具有較清楚的幾何及輻射訊息，例如幾何結構線條或在輻射特徵表現的邊緣線 (Edge Lines)，也因此提供偵測及匹配上的便利(彭念豪、趙鍵哲，2004)。在電腦視覺(Computer Vision) 領域，直線特徵用於相機內方位率定、遁點偵測 (Vanishing Point Detection)及物件模型之應用相當

普遍，甚至著有專書論述(Faugeras, 1993; Hartley and Zisserman, 2000)。在攝影測量界約莫自 1980 年代起，部份之研究者開始涉獵於直線特徵用於航攝及近景像片之應用研究。線形特徵用於方位求解 問題最早是由 Masry(1981)所提出，Liu et al.(1990) 利用直線控制特徵以兩步驟方式求解像片外方位參數，Habib(1998)結合點及線形特徵進行空中三角解算，Zalmanson(2000)利用三維參數化之曲線特徵求解像片外方位參數，李志宏(2002)利用共面直線特徵物求解相對方位參數及物型重建， Schenk(2004)利用連結線(Tie Lines)聯結像片，並利用少數的控制直線進行空中三角解算。光達資料因具有空間三維訊息，地、物表之斷 1_{國立臺灣大學土木工程學系助理教授} 2_{國立臺灣大學土木工程學系碩士班研究生} 收到日期:民國 93 年 09 月 30 日修改日期:民國 94 年 02 月 02 日接受日期:民國 94 年 02 月 16 日

(2)

線及結構線可藉由影像處理或結構化 (Structuralization)方式推演產生，所差別的是因處理之程序造成誤差的等級不同。除了光達本身固有的誤差外，以影像處理之邊緣線偵測方式獲致直線特徵必須依賴網格式資料，而建構網格資料之平滑處理常伴隨降低直線特徵空間精度，但其優勢是方法較為普及簡單，因此本研究工作依此簡易程序，將光達離散的點雲資料利用內插方式組成規則網格點，並以 Canny 邊緣線偵測元決定物空間中的直線特徵。本研究假設物空間場景具豐富的三維直線資料(例如都會地區)，相對於利用雙像或多重影像，以單像求解外方位參數不須考量像空間共軛像對之匹配對應(雖然其呈現較強之求解幾何)，此種作業方式尤其助益於排除近景攝影呈現顯著遮蔽之共軛特徵對應。在以光達資料作為控制資料求解的數學模式中，考慮到物空間資料的隨機性，在平差模式中將控制直線視為觀測量處理，而控制直線之物、像及投影中心幾何關聯採取以直線為基礎之共線式，方便於物、像觀測量平差後改正數之分析。

2. 研究方法

2.1 直線特徵於物空間之參數化

以控制直線求解像片外方位參數，必須收集直線特徵的像片觀測量及物空間控制資料。描述三維物空間直線之幾何，最簡單的方式為利用兩點構成一直線，即使用六參數來描述一條直線，但經過兩個條件的約制:(1).沿著直線方向與 X-Y、Y-Z 或 X-Z 平面之交點，稱為穿刺點(Penetration Point)，則描述穿刺點的一座標分量 0，因此以穿刺點為參考點只需使用 2 個參數，即式(1)之

p

與

q

;(2).將方向向量除以三分量中之其中一分量，此時方向向量中有一分量固定為 1，因此只需使用 2 個參數，即式(1) 之

a

與

b

，如此三維直線可以四個獨立參數 (

a

,

b

,

p

,

q

)來組成三維直線參數式(Habib et al., 2000)，如式(2)及示意圖 1(a)。 0 0 0 0 1 L L L X X dX Y Y t dY Z Z dZ p a q z b ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ₌⎢ ⎥_{+ ×}⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ =_{⎢ ⎥}+ ×_{⎢ ⎥} ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ (1) L L L

X

p

z

a

Y

q

z

b

Z

z

=

+

×

⎧

⎪

₌

₊

_×

⎨

⎪

₌

⎩

(2) 其中

X Y Z

₀

,

₀

,

₀和

X Y Z

_L

,

_L

,

_L為物空間直線上兩點；

dZ

dY

dX

,

為方向向量；

q

p

b

a

,

為描述三維直線的四參數；

,

t z

為三維直線上之變數；式(2)並非一通式，例如當物空間中直線平行於 X-Y 平面，則穿刺點應在 X-Z 或 Y-Z 平面上。除此之外，空間三維直線亦可能在 X-Y、Y-Z 及 X-Z 平面上均有交點，所以在決定此直線參數式時，須先找尋合適的平面穿刺點。穿刺點既為三維直線與前述三平面之交點，考量三維直線資料固有的誤差特性，從交會幾何的觀點來看，三維直線與所選擇之平面交角越接近垂直，則獲致之穿刺點參數

(

p

,

q

)

精度越高，故本研究在擇取穿刺點所採用的方式為判斷三維直線和 X-Y、Y-Z 及 X-Z 三個平面的法向量交角，選擇交角最小的平面為穿刺點所在平面，交角越小者表示該三維直線和此相應平面越接近正交，則其穿刺點參數精度越高。

(3)

(a) (b) 圖 1、(a)為直線參數式之示意圖，(b)為直線共線式之示意圖

2.2 以直線為基礎之共線式

(Line‐based Collinearity

Equations)

以直線為基礎的外方位參數求解方式分為兩種方法，一為利用共面的條件，另一法為將直線方程式代入共線式求解，兩者為等效的解算模式(彭念豪、趙鍵哲，2004)，但在實際應用上，後者具有可擴展直線特徵至三維曲線特徵、適合於自動作業流程及可直接獲取改正數訊息等優勢(Schenk, 2004)，故本研究採用以直線方程式代入共線式的解算模式。以直線為基礎之共線模式和傳統以點為基礎之共線模式相似，是描述投影中心、像片量測點及物空間對應點三者共線的關係，惟物空間直線上點位以三維直線參數式描述，讓平差模式自行解算像片量測點對應在物空間直線上的點軌跡。將直線參數式(2)代入傳統之共線方程式，產生以直線為基礎之共線方程式，如式(3)及示意如圖 1(b)。 (3) 其中

X Y Z

c

,

c

,

c為投影中心物空間座標；

x y

,

為量測點框標座標；

x

0

, y

0為像主點框標座標；

f

為相機焦距； ij

r

為旋轉矩陣元素，旋轉元素定義為從物空間轉至像空間，旋轉順序為

ω

,

ϕ

,

κ

；相較於以點為基礎之共線式，式(3)每一像片觀測點多了一個未知參數 z，利用此 z 值可計算所觀測之像點在物空間直線上之對應位置。在此須強調的是，對平差具有約制效力的是像空間直線與物空間直線的幾何對應，因此平差前像空間直線觀測量與物空間控制直線並無點對點的必要對應條件，此模式下求算的 z 值是平差後反映的物、像點對應狀況。考慮最少對應關係(亦為本研究工作目前的作業方式)，式(3)每條控制直線包含兩個線段端點的像片點位，同時引入兩個 z 未知數，一般情形在內方位已知的條件下，假使有 m 條控制直線，則單像求解外方位參數之未知參數總個數為 (2m+6)，而每條控制直線可列 4 組觀測方程式，考慮求解條件

4 m

≥

2 m

+

6

，此模式至少需要 3 條有效的控制直線才能求解。以幾何方式解析以控制直線進行單像後方交會，若將通過控制直線物、像的平面視為投影中心軌跡，則決定空間中投影中心的最小條件必須至少有三個相對應的投影中心軌跡平面且此三平面必須交於一點(投影中心)，因此當控制直線分佈無法造就此一交點或造就無限交點，例如在物空間交於屋角點的三直線(其三軌跡平面之交集為通過物、像屋角點之空間直線)，將形成奇異解。式(3)線性化後之觀測方程式表為式(4)，其中係數矩陣

A

₁如式(5)。 1 1 1 1 4 1 4 (2 6) (2 6) 1 4 1

,

m m m m m

l

A

ξ

e

× × + + × ×

=

+

2 1 1

~ (0,

1 0 1

)

e

Σ =

σ

P

− (4) 11 12 13 0 31 32 33 21 22 23 0 31 32 33 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) c c c c c c c c c c c c p z a X r q z b Y r z Z r x x f p z a X r q z b Y r z Z r p z a X r q z b Y r z Z r y y f p z a X r q z b Y r z Z r + × − × + + × − × + − × ⎧ − =− × ⎪ _{+ × − × + + × − × + − ×} ⎪ ⎨ _{+ × − × + + × − × + − ×} ⎪ − =− × ⎪ + × − × + + × − × + − × ⎩

(4)

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0 c c c c c c c c c c c c x x x x x x x X Y Z z y y y y y y y X Y Z z A x x x x x x x X Y Z z y y y y y y y X Y Z z ω ϕ κ ω ϕ κ ω ϕ κ ω ϕ κ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎡ ⎤ ⎢_∂ _∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎢ ⎥ ⎢_∂ _∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎢ ⎥ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎢_∂ _∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎢ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎢∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎢ ⎢⎣ ⎦ L L L L M M M O ⎥ ⎥ ⎥ ⎥ (5) 簡化式(3)的表示式，令 0

,

0

U

V

x

f

y

f

W

−

= − ×

−

= − ×

，則式(5)中偏微分之各項如下：依照最小二乘法(Least-squares Method)並進行漸近解算，其收斂解如下： 1 1 1 1 1 1 1 1 1 1 1 1 ˆ ₍ T _; T ₎ N U N A P A U A Pl ξ ₌ − ₌ ₌ (6) 2 1 1 1 0

ˆ

4 (2

6)

T

e Pe

m

σ

=

−

+

%

(7) 1 1 2 ˆ 0 1

ˆ

_{ˆ N}

ξ

σ

−

Σ =

(8) 其中

l

₁為約化後的觀測量向量

A

₁為係數矩陣； 1

[

0

,

0

,

0

, , , , ,

1 2

,...,

2

]

T m

X Y Z

z z

z

ξ

=

ω ϕ κ

為未知參數向量；

ξ

ˆ

₁為未知參數估值向量；

e

1為誤差向量；

e

%

1為改正數向量；

σ

02為先驗單位權方差； 2 0

ˆ

σ

為後驗單位權方差； 1

P

為權矩陣；

m

為控制直線個數； ˆ

ˆ

ξ

Σ

為後驗參數方差-協方差矩陣；

2.3 非固定式約制

(Soft Constraint)

一般而言，單片後方交會是將像片量測值作為平差時的觀測值，而視地面控制資料為完全無誤差 (即固定其座標值)。但考慮到控制資料的隨機特性，為了合理的誤差分析，必須將已知的控制資料依照其精度施以非固定式的約制，視為另一類之觀測值，加入平差模式(吳怡燊，2003)。描述控制直線之參數為

( , , , )

p q a b

，其觀測方程式如式(9)，而約化後之觀測方程式如式(10)。 0 0 2 0 2 0 i i i i p i i i q _i _i _i i a i i i i i b

l

_p

l

q

e

a

l

a

b

l

⎡ ⎤

_{⎡ ⎤}

_⎡

_{+ ∆}

_⎤

⎢ ⎥

_{⎢ ⎥}

_⎢

_⎥

+ ∆

⎢ ⎥

₌

_{⎢ ⎥}

_{+ =}

_⎢

_⎥

₊

⎢ ⎥

_{⎢ ⎥}

_⎢

_{+ ∆}

_⎥

⎢ ⎥

_{⎢ ⎥}

_⎢

_⎥

⎢ ⎥

⎣ ⎦

⎢

⎣

+ ∆

⎥

⎦

⎣ ⎦

(9) 1 1 0 1 0 1 2 0 2 2 4 1 0 1 1 1 1 1 1 1 1 i i i i i i i p i i q i m a i m i b _m m m p q a l p _p b l q q l e e a l a p b l b _q a b φ φ × ∆ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ _{⎥ ∆} ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ∆ ⎡ − ⎤ _⎡_∆ _⎤ _⎢ _⎥_⎢ _⎥ ⎢ ⎥ _⎢ _⎥ _∆ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ _⎢∆ _⎥ _⎢ _⎥_⎢ _⎥ =⎢ _{⎥ ⎢ ⎥}= _∆ + = + ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ _⎢ _⎥ _⎢ _⎥_⎢_∆ _⎥ ⎢ − ⎥ ⎢⎣∆ ⎥⎦ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ∆ ⎣ ⎦ _⎢ _⎥_⎢ _⎥ ∆ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ _{⎥ ∆} ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ L L (10) 2 2 2

,

A

ξ

e

=

+

e

2

~ (0,

2

σ

02

P

21

)

−

Σ =

)) 1 ( 1 ( )) ( ) ( ) ( ) ( ) ( ) (( ))) cos sin ( ) ( ) sin (sin ) ( ) (cos ) (( ) cos cos cos ( ) ( ) sin cos (cos ) ( ) sin cos ( ) (( ))) ( ) ( ) ( ) (( ) ( ) ( ) ( ) (( ) ( ) ( ) ( 33 32 31 13 12 11 23 22 21 32 33 12 13 33 13 32 12 31 11 r r b r a W U r r b r a W f z x r Z Z r Y Y r X X W f x Z Z Y Y X X W U Z Z Y Y X X W f x r Z Z r Y Y W U r Z Z r Y Y W f x r W U r W f Z x r W U r W f Y x r W U r W f X x c L c L c L c L c L c L c L c L c L c L c L c L c L c c c × + × + × − × + × + × × − = ∂ ∂ × − + × − + × − × − = ∂ ∂ − × − + × − + × − − − × − + × − + − × − × − = ∂ ∂ × − + − × − − × − + − × − × − = ∂ ∂ − × = ∂ ∂ − × = ∂ ∂ − × = ∂ ∂ κ ω ϕ ω ϕ ϕ ω ϕ κ ω ϕ κ ϕ κ ϕ ω )) 1 ( 1 ( )) ( ) ( ) ( ) ( ) ( ) (( ))) cos sin ( ) ( ) sin (sin ) ( ) (cos ) (( ) cos cos (sin ) ( ) sin cos sin ( ) ( ) sin (sin ) (( ))) ( ) ( ) ( ) (( ) ( ) ( ) ( ) (( ) ( ) ( ) ( 33 32 31 23 22 21 13 12 11 32 33 22 23 33 23 32 22 31 21 r r b r a W V r r b r a W f z y r Z Z r Y Y r X X W f y Z Z Y Y X X W V Z Z Y Y X X W f y r Z Z r Y Y W V r Z Z r Y Y W f y r W V r W f Z y r W V r W f Y y r W V r W f X y c L c L c L c L c L c L c L c L c L c L c L c L c L c c c × + × + × − × + × + × × − = ∂ ∂ − × − + − × − + − × − × − = ∂ ∂ − × − + × − + × − − × − + − × − + × − × − = ∂ ∂ × − + − × − − × − + − × − × − = ∂ ∂ − × = ∂ ∂ − × = ∂ ∂ − × = ∂ ∂ κ ω ϕ ω ϕ ϕ ω ϕ κ ω ϕ κ ϕ κ ϕ ω

(5)

其中

p

_i、

q

_i、

a

_i、

b

_i為求解參數； 0 i

p

、

q

_i0、

a

_i0、

b

_i0為已知之控制直線參數之疊代近似值； i

p

∆

、

∆

q

_i、

∆

a

_i、

∆

b

_i為求解之參數增量；平差始於令 0 i i p

p

=

l

、 0 i i q

q

=

l

、 0 i i a

a

=

l

、 0 i i b

b

=

l

，但隨著平差模式的漸近計算，

p

_i0、 0 i

q

、

a

_i0、

b

_i0近似值變動，約化後之觀測值向量

l

₂也會隨著改變。結合像點觀測量

l

1及控制直線參數觀測量

l

2 兩類之觀測方程式後的平差模式中的 A 矩陣將擴充為式(11)。 1 3 4 (2 6) 4 4 8 (6 6) 2 4 4

0

m m m m m m m m

A

× + × × + ×

⎡

⎤

⎢

⎥

= ⎢

⎥

⎢

⎥

⎣

⎦

(11) 式(11)中

A

₁與式(4)中

A

₁相同，為式(3)對 6 個外方位參數(X Y Z_C, _C, _C, , , )

ω ϕ κ

及 z 偏微分所得之係數矩陣，

A

₃ 為式 (3) 對控制直線參數

( , , , )

p q a b

偏微分所得之係數矩陣，如式(12)，

A

₂ 為式(10)中之

A

₂，為一單位矩陣。從式(11)中計算多餘觀測數為(2m-6)，此數目與將控制直線視為固定之系統多餘觀測數是一致的。 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 x x x x p q a b y y y y p q a b x x x x A p q a b y y y y p q a b ⎡∂ ∂ ∂ ∂ ⎤ ⎢_∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎢∂ ∂ ∂ ∂ ⎥ ⎢_∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎢ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ = ∂ ∂ ∂ ∂ ⎢ ⎥ ⎢ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ∂ ∂ ∂ ∂ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ L L L L O M M M M M M M M (12) 式(12)中偏微分之各項如下： 31 11 ( U r ) x f r p W W × ∂ ₌− _× ₋ ∂ , (12 32) U r x f r q W W × ∂ ₌− _× ₋ ∂ 31 11 ( U r ) x f z r a W W × ∂ ₌− × _× ₋ ∂ ， (12 32) U r x f z r b W W × ∂ ₌− × _× ₋ ∂ 31 21 ( V r ) y f r p W W × ∂ ₌− _× ₋ ∂ ， 32 22 ( V r ) y f r q W W × ∂ ₌− _× ₋ ∂ 31 21 ( V r ) y f z r a W W × ∂ ₌− × _× ₋ ∂ ， 32 22 ( V r ) y f z r b W W × ∂ ₌− × _× ₋ ∂

2.4 控制直線參數觀測量的方

差-協方差矩陣

控制直線參數觀測量

(

p q a b

_i

, , , )

_i _i _i 之方差-協方差矩陣

Σ

₂須視來源而定，在本研究中利用已知物空間線段兩端點方式推演的直線參數觀測量，採用一條穿刺點位於 X-Y 平面上的三維直線為例(當所選之穿刺點位於不同平面時，推導出之式(13)、(14)中之參數將不相同)，從式(1)可推得如式(13)，其相應之方差-協方差矩陣可由針對各項觀測值偏微分所得之 B 矩陣經誤差傳播 (Error Propagation)方式推得，如式(15)所示。 1 1 2 1 2 1 1 1 2 1 2 1 2 1 2 1 2 1 2 1

(

) /(

)

(

) /(

)

(

) /(

)

(

) /(

)

p

X

Z

X

Z

q

Y

Z

Y

Z

a

X

Z

b

Y

Z

=

−

×

−

⎧

⎪ = − ×

−

⎪

⎨ =

−

⎪

⎪ =

−

⎩

(13) 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 p p p p p p X Y Z X Y Z q q q q q q X Y Z X Y Z B a a a a a a X Y Z X Y Z b b b b b b X Y Z X Y Z ∂ ∂ ∂ ∂ ∂ ∂ ⎡ ⎤ ⎢_∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ⎢ ⎥ ⎢_∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ = ⎢ ∂ ∂ ∂ ∂ ∂ ∂ ⎥ ⎢_∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎢ ∂ ∂ ∂ ∂ ∂ ∂ ⎥ ⎢_∂ _∂ _∂ _∂ _∂ _∂ ⎥ ⎢ ⎥ ⎣ ⎦ 2 2 2 1 1 1 2 1 2 2 2 1 2 1 2 1 2 1 2 2 2 1 1 1 2 1 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 1 ( ) ( ) 0 0 ( ) ( ) ( ) ( ) 0 0 ( ) ( ) 1 ( ) 1 ( ) 0 0 ( ) ( ) 1 ( ) 1 ( ) 0 0 ( ) ( ) Z Z X X Z Z X X Z Z Z Z Z Z Z Z Z Z Y Y Z Z Y Y Z Z Z Z Z Z Z Z X X X X Z Z Z Z Z Z Z Z Y Y Y Y Z Z Z Z Z Z Z Z − − − − ⎡ ⎢ ₋ ₋ ₋ ₋ − − − − − − − − = − − − − − − − − − − − − − − − − ⎣ ⎤ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎦ (14)

(6)

1 1 1 1 1 1 2 1 2 1 2 1 1 1 1 2 1 2 1 2 1 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ( , , , ) 2 2 4 4 2 2 2 x X XY XZ X X XY XZ Y YZ YX YY YZ p pq pa pb Z ZX ZY ZZ qp q qa qb pqab X X Y X Z ap aq a ab Y YZ bp bq ba b Z B sym σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ × ⎡ ⎤ ⎢ ⎥ ⎡ _{⎤ ⎢} _⎥ ⎢ _{⎥ ⎢} _⎥ ⎢ ⎥ Σ =_⎢ = ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ _{⎥ ⎢} _⎥ ⎢ ⎥ ⎣ _{⎦ ⎢} _⎥ ⎢ ⎥ ⎣ ⎦ T B (15) 假設共有 m 條控制直線且所有控制直線之間不相關，則其相應之方差-協方差矩陣

Σ

₂可表為如式(16)。 1 1 1 1 ( , , , ) 2 1 2 0 2 ( m,m, m,m) p q a b p q a b P φ σ φ − ⎡Σ ⎤ ⎢ ⎥ Σ =_⎢ _⎥= ⎢ _Σ ⎥ ⎣ ⎦ L (16)

2.5 考量控制直線誤差之函式

及隨機模式

考慮像點量測誤差及控制直線誤差之整體觀測方程式可以式(17)表示。 1 1 3 1 2

0

2 2

l

A

l

e

l

A

ξ

⎡ ⎤ ⎡

⎤ ⎡ ⎤

=

_{⎢ ⎥ ⎢}

=

_{⎥ ⎢ ⎥}

+

⎣ ⎦ ⎣

⎦ ⎣ ⎦

A

ξ

e

=

+

1 2 1 0 2

0 , ~ (0,

)

0 e

Σ =

⎡

_⎢

Σ

⎤

_⎥

=

σ

P

−

Σ

⎣

⎦

(17) 其中，

[

0 0

, , , , , , , , , , , , , , , , , ,

0 1 2 2 1 1 1 1

]

T m m m m m

X Y Z

z z

z p q a b

p q a b

ξ

=

ωϕκ

_L

，依最小二乘法漸近解算，其收斂解如下： 1

ˆ

₍

T

_;

T

₎

N U

N

A PA U

A Pl

ξ

₌

−

₌

(18) 2 0

ˆ

2

6

T

e Pe

m

σ

=

−

%

(19) 2 1 ˆ 0

ˆ

_{ˆ N}

ξ

σ

−

Σ =

(20)

3. 實驗成果與分析

為說明及分析利用光達資料萃取之控制直線 用於單片後方交會求解之成效，本研究工作進行以下兩組實驗，分別以空載(實驗一)及地面(實驗二) 光達資料用於航空影像及近景影像之外方位參數解算及精度分析。作業方式為將離散光達點雲資料經由 Surfer 軟體以克利金法(Kriging Interpolator)內插規則網格，其網格大小約與點雲密度相同，再以 Canny 邊緣線偵測元偵測出線特徵，藉由擬合或連接兩端點的方式得到控制直線，而影像上之控制直線現階段則以人工方式量測。

3.1 實驗一：空載光達

實驗一採用利用空載光達所量測的新竹地區三維點雲資料(史天元等，2002)及對應的航照影像 (圖 2(a)、(b)、(c))。空載光達的點位原始高程精度預估為 15cm，點位平面(X,Y)精度各約 40cm。如 2-4 節所述方法，控制直線的精度可經由誤差傳播推算，但在實際作業中發現，控制直線因大多位於高差明顯之處，受網格化平滑的影響顯著，所以需考量給予較差的點位精度。實驗一的各項配置如表 1，控制直線分佈參見圖 2(b)、(c)。表 1、實驗一的配置像片比例尺

≈

1/5400 像元地面解析度

≈

25 cm

×

25 cm

像片觀測量先驗中誤差 1 pixel =

0.05mm

光達網格密度

1 m

×

1 m

光達(X,Y,Z)先驗中誤差 (

0.5m

,

0.5m

,

0.2m

) 採用控制直線數目左、右像各 12 條

(7)

(a) 0 0 0 0 0 0 0 0 0 0 (b)○：控制直線分佈區堿 (c)○：控制直線分佈區堿

圖 2、(a)空載光達之 Range Image，(b)航照影像左像，(c) 航照影像右像

3.1.1 單片後方交會計算分析

物、像空間之直線觀測量其隨機特性在單片外方位參數求解的平差過程中，可藉由改正數見其端倪，進而評估合理性，並作為偵錯依據。圖 3 展示利用實驗一計算出之外方位參數將控制直線投影回影像上之情形，圖中藍色線條為像片直線的觀測量，紅色線條為初始之控制直線(即未經過平差計算之值)，綠色線條為經平差模式改正後之控制直線投影回影像的位置，圖 3 中(a)為左像像片觀測量改正數最大處；(b)為左像直線參數改正數最大處；(c)為右像像片觀測量改正數最大處；(d)為右像直線參數改正數最大處，所有控制直線物、像空間觀測量平差成果之改正數統計量如表 2。 3340 3345 3350 3355 3360 3365 3370 3375 1235 1240 1245 1250 1255 1260 1265 1270 (a) 4530 4535 4540 4545 4550 4555 4560 4565 670 675 680 685 690 695 700 705 (b) 2040 2045 2050 2055 2060 2065 2070 2075 1245 1250 1255 1260 1265 1270 1275 1280 (c) 2430 2435 2440 2445 2450 2455 2460 2465 435 440 445 450 455 460 465 (d) 圖 3、平差前後之控制直線

(8)

本實驗中兩張影像之個別單片後方交會平差結果在 95%信心區間內以

χ

2分佈統計檢定零假說：

σ

ˆ

₀

=

σ

₀，其結果均顯示後驗單位權中誤差與先驗單位權中誤無顯著區別(參見表 2)，而像點改正量不及原先預估之誤差量，其現象與彭念豪及趙鍵哲(2004)之模擬測試結果相呼應，物點改正量亦小於原先預估之誤差量。另外，值得注意的是景深方向的改正量相對不敏感。

3.1.2 前交共軛光束之相對一致性

檢查

本研究工作利用控制直線以單片後方交會方式解算外方位參數(A)，相較於傳統空中三角作業採用連結點以共軛光束加強相鄰像片間外方位參數(B)一致性之方式，針對左、右像共軛光束在物空間交會的誤差量，可在相對模式下利用 y 視差量作為評估，本實驗採用以影像匹配方式選取兩張影像內高精度、可靠之 20 個均勻分佈之共軛點(圖 4)，並利用上述 A、B 之方位參數計算 y 視差量，其結果如表 3 所示，y 視差量統計如圖 5 所示。表 3、y 視差量評估最大值最小值平均值標準差 A 之 y 視差量 (pixel) 2.70 0.02 1.03 0.74 B 之 y 視差量 (pixel) 0.47 0.01 0.18 0.14 左像右像後驗單位權中誤差

σ

ˆ

₀ 1.01 1.10 2

χ

分佈測試(在 95%信心區間內)

[

]

12.17 ∈

8.23, 31.53

14.41 ∈

[

8.23, 31.53

]

最大值 (絕對值) 最小值 (絕對值) 標準差最大值 (絕對值) 最小值 (絕對值) 標準差像片觀測量改正數(pixel) 0.74 0.02 0.29 0.79 0.01 0.32

p

(m) 12.014 0.200 4.684 9.051 0.351 3.167

q

(m) 0.605 0.003 0.212 0.310 0.001 0.123

a

0.051 0.003 0.030 0.056 0.001 0.029 直線參數改正數

b

0.003 0.000 0.001 0.003 0.000 0.001 X 方向 0.657 0.015 0.285 0.738 0.005 0.306 Y 方向 0.624 0.021 0.272 0.622 0.006 0.300 物點改正量(m) Z 方向 0.030 0.000 0.012 0.025 0.000 0.011 表 2、實驗一改正數統計

(9)

圖 4、右像共軛點分佈 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 pixel (a) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 pixel (b)

圖 5、y 視差統計圖，(a)對照 A，(b)對照 B 如所預期，單像求解外方位參數因未如傳統空中三角作業加諸共軛光束之約制，另加上控制直線本身之精度等級，其 y 視差量明顯大於空中三角成果之 y 視差量。但控制直線精度之提昇有助於改善此一現象，本研究工作利用空中三角成果萃取控制直線，將這些控制直線資料引入本工作之單像外方位參數求解中，並據此評估 y 視差量，如表 4 所示。經由較高精度的控制直線所計算出的外方位參數，其共軛光束之一致性確實可獲得顯著改善。表 4、y 視差量評估(利用精度較高之控制直線) 最大值最小值平均值標準差 y 視差量 (pixel) 0.45 0.00 0.21 0.14

3.1.3 外部精度檢查

前述兩種程序檢視平差成果之合理性及內部精度之一致性，而在外部精度評估的考量上，由於此實驗光達資料相對於與航測資料的控制基準存有一待查明之系統性問題或光達系統未率除完善而殘留顯著的系統性誤差(周富晨，2004)，因此不適合進行外部精度評估。

3.2 實驗二：地面光達

實驗二中所採用的資料為成功大學圖書館側面的近景影像及地面光達資料(劉燈烈，2004)，所使用之地面雷射掃描儀為加拿大 Optech 公司所生產之 ILRIS-3D ，而影像由數位相機 FujiFilm FinePix S1 Pro 搭配 Nikon 28mm f/2.8s 鏡頭取得，實驗區光達資料及左、右像兩張近景影像參見圖 6(a) 、 (b) 、 (c) 。地面光達的定位精度約為 8mm-12mm(劉燈烈，2004)。實驗二之配置如表 5。表 5、實驗二之配置像片比例尺

≈

1/1400 像元地面解析度

≈

1 cm

×

1 cm

像片觀測量先驗中誤差 1 pixel =

0.007mm

光達網格密度

3 cm

×

3 cm

光達(X,Y,Z)先驗中誤差 (

0.008m

,

0.008m

,

0.008m

) 採用控制直線數目 12 條 18 條

(10)

-20 -15 -10 -5 0 5 10 15 -10 -5 0 5 10 (a) (b)紅線為所採用之控制直線，黃圈為共軛點，藍點為地面檢核點 (c)紅線為所採用之控制直線，黃圈為共軛點，藍點為地面檢核點

圖 6、(a)實驗二地面光達之 Range Image，地面近景影像左像(b)及右像(c)

3.2.1 單片後方交會計算分析

本實驗區之物空間座標系統使用地面光達測站座標系統(劉燈烈，2004)。表 6 呈現實驗二所有控制直線物、像空間觀測量平差成果之改正數統計量。零假說：

σ

ˆ

₀

=

σ

₀，在 95%信心區間內，顯示後驗單位權中誤差與先驗單位權中誤無顯著區別；物、像點改正量仍顯示出不及預估之誤差量；另外，景深方向(左像:約為 X 方向; 右像:Y 方向) 改正量亦相對不敏感。

3.2.2 前交共軛光束之相對一致性

檢查

利用左右像重疊區選定之 9 個共軛點(參見圖 6)，以前方交會來檢驗光束相交情形，成果統計如表 7。本實驗區因未有其他空中三角資料可資比對，因此只評估單像求解外方位參數前交之 y 視差量。類如實驗一之結果，實驗二之控制直線精度直接影響 y 視差量，若能提昇控制直線的品質(例如經由結構化程序的特徵萃取)，將有助於降低 y 視差量。

3.2.3 外部精度檢查

本實驗例於影像重疊區域中僅提供有 3 個地面量測之檢核點(圖 6)可供檢驗外方位參數之物空間定位品質。為使光達資料與地面檢核點具有一致的物空間三維座標系統，本研究援用劉燈烈(2004) 測站連結獨立模型法之成果以七參數實施物空間三維座標轉換後以本文方法進行外方位參數求解。外方位參數物空間前交的定位精度檢核成果如表 8。受限於實驗資料中所提供之檢核點數太少，其均方根誤差值(R.M.S.E.)雖不具統計意義，但可見其方位參數呈現公分級之物空間定位誤差，而此等誤差來源除控制直線物、像空間量測誤差之影響外，尚受限於三維座標轉換的七參數估值之精度。

(11)

表 6、實驗二改正數統計左像右像後驗單位權中誤差

σ

ˆ

₀ 1.26 1.10 2

χ

分佈測試(在 95% 信心區間內)

28.43 ∈

[

8.23, 31.53

]

36.60 ∈

[

16.79, 46.98

]

最大值 (絕對值) 最小值 (絕對值) 標準差最大值 (絕對值) 最小值 (絕對值) 標準差像片觀測量改正數 (pixel) 2.682 0.0044 0.628 1.884 0.000 0.525

p

(m) 0.034 0.000 0.012 0.363 0.000 0.124

q

(m) 0.211 0.001 0.072 0.478 0.000 0.155

a

0.004 0.000 0.001 0.028 0.000 0.010 直線參數改正數

b

0.025 0.000 0.009 0.012 0.000 0.006 X 方向 0.002 0.000 0.001 0.015 0.000 0.005 Y 方向 0.011 0.000 0.005 0.003 0.000 0.001 物點改正量(m) Z 方向 0.007 0.000 0.003 0.011 0.000 0.004 表 7、y 視差量評估最大值最小值平均值標準差 y 視差量(pixel) 3.32 0.03 1.15 1.19 表 8、實驗二成果經前方交會後與地面檢核點 (視為真值)比較之定位成果 E N H R.M.S.E.(m) 0.070 0.055 0.074

4. 結論與建議

本研究探討利用光達資料之控制直線解算單張影像外方位參數之數學模式，並以空載及地面資料實例分析定位精度，其內部精度檢核結果顯示融合光達資料與影像外方位求解之可行性。本研究屬先期階段，諸多考量不近完善，若干工作仍留後續深入研究： 1. 分析網格化處理對直線特徵幾何精度的影響，據以合理提供控制資料相應的隨機誤差。 2. 除了分析以影像處理方式萃取直線特徵之效益外，亦應朝向利用更精確的偵測方式 (例如結構化的特徵偵測)產生較高精度的光達資料控制直線，有助於提昇外方位參數品質。 3. 發展自動化的物、像控制直線量測及匹配，提昇作業效率。 4. 在平差模式中引入方差分量模式(Variance Component Model)，反推較合理之光達資料原始精度。

(12)

5. 在光達資料及航測資料具共同的基準底下並經適當程序濾除光達系統性誤差，利用擁有足夠地真資料的實驗場，進行外部精度檢核、評估。 6. 將控制直線引入衛照影像之外方位參數求解。

致謝

本研究感謝農委會提供空載光達資料、交通大學史天元教授提供之光達資料處理經驗、及成功大學曾義星教授研究群提供之地面光達、近景影像資料及地面控制與檢核資料，並承蒙行政院國家科學委員會支持之研究計劃 (NSC 92-2211-E-002-095 , NSC 93 – 2211 – E-002-054) 經費支援。兩位審查委員提供精細的修正意見增加本文的完整性，謹表謝忱。

參考文獻

史天元、彭淼祥、徐偉城，2002。應用空載雷射掃描儀進行地震災區形變研究，行政院農業委員會九十一年度，試驗研究計畫研究報告，計畫編號 91 農科-5.1.1-林 -R1(8)。李志宏，2002。應用線特徵物求解攝影測量方位參數與物型重建，國立台灣大學土木研究所碩士論文。吳怡燊，2003。影像控制區塊自動化量測及方位求解作業之探討，國立台灣大學土木研究所碩士論文。周富晨，2004。適應性點雲過濾演算法於空載光達資料產生數值高程模型之研究，國立成功大學測量及空間資訊研究所碩士論文。彭念豪、趙鍵哲，2004。以控制直線求解單張像片外方位參數模式之探討，第二十三屆測量學術及應用研討會論文集， pp.197-204。劉燈烈，2004。地面光達點雲資料的平差結合與影像敷貼，國立成功大學測量及空間資訊研究所碩士論文。

Faugeras, O., 1993. Three-dimensional computer vision, The MIT Press, Cambridge, Massachusetts, London, England.

Habib, A. F., 1998. Aerial Triangulation Using Point and Linear Features, International Archives of Photogrammetry and Remote Sensing, Munich/Germany, Vol.XXXIII, Part 3-2W5, pp.137—141.

Habib, A., A. Asmamaw, D. Kelly, and M. May, 2000. Linear Features in Photogrammetry, Report No. 450, Department of Civil & Env. Eng. & Geodetic Science, The Ohio State University, Columbus, OH.

Hartley, R. and A. Zisserman, 2000. Multiple view geometry in computer vision, Cambridge University Press, Cambridge. Liu, Y., T.S. Huang, and O. Faugeras, 1990.

Determination of Camera Location from 2-D to 3-D Line and Point Correspondences, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol.12, No.1, pp.28—37.

Masry, S.E., 1981. Digital Mapping Using Entities: A New Concept, Photogrammetric Engineering and Remote Sensing, Vol. 48, No.11, pp.1561—1599.

(13)

Schenk, T., 2004. From point-based to feature-based aerial triangulation, ISPRS Journal of Photogrammetry & Remote Sensing, Vol. 58, pp.315—329.

Zalmanson, G. H., 2000. Hierarchical Recovery of Exterior Orientation from Parametric and Natural 3-D Curves, International Archives of Photogrammetry and Remote Sensing, Vol. XXXIII, Part B2, pp.610—617.

(14)

LIDAR-based Control Lines for Single

Photo Orientation

Jen-Jer Jaw

1

Nei-Hao Perng

2

ABSTRACT

Airborne and ground-based LIDAR systems are able to quickly support the tasks of three dimensional spatial data acquisition and provide users with abundant geo-spatial information. Optical images, on the other hand, recorded with adequate radiometric detail offer advantages for better interpretation of the scene, three dimensional coordinates of object points can only be determined after orientation and intersection steps, however. Thus, one of the apparent gains in fusing LIDAR data sets with photogrammetry is to solve the photo orientation by using control information resulting from extracting features from LIDAR data sets.

This study aims at employing control lines extracted from LIDAR data sets solving for single photo orientation. The functional model describing the geometric relationships among perspective center, line observations in image space, and that in object space, the random model characterizing the control line measuring errors both in image space and object space are established and tested in this work. The experimental results suggest the application potential of fusing LIDAR data sets with images for the task of single photo orientation.

Key Words：

Control Lines, Single Photo Orientation

1

Assistant Professor, Department of Civil Engineering, National Taiwan University

2_{Graduate student, Department of Civil Engineering, National Taiwan University}

Received Date: Sep. 30, 2004 Revised Date: Feb 02, 2005 Accepted Date: Feb 16, 2005