以現有正射影像與 DTM 資料解算外方位參數

(1)

Volume 9, 4, December 2004, pp. 13-20

以現有正射影像與 DTM 資料解算外方位參數

楊靈峰

¹

廖揚清

²

摘要

航照像片之外方位參數求定一直是攝影測量中相當基本的課題，隨著電腦科技的進步，配合影像處理技術，許多數值航測的工作可以朝著自動化的方向發展。對於較小的區域範圍，若能利用現有的資料，將區域內已經變遷的地貌局部更新，而非以一般數值航測的方式進行，勢必可以節省時間與人力。對於有部分變遷的航攝影像，在本研究中，利用了區域匹配與蒙地卡羅方法，來設計一套半自動求解該航攝影像外方位之作業流程，並可在個人電腦上操作。藉由此方式，求定了該航攝影像的外方位最或是值之後，再配合DTM 資料，來進行該地區的正射影像更新，並且其結果可達到與一般數值航測作業相同之之精度要求。

關鍵字：正射影像、數值地形模型、蒙地卡羅法、局部更新

1. 前言

在一張空照影像的涵蓋範圍內，若能利用現有的資料，將範圍內部分已經變遷的地貌更新(Local Regeneration)，而非以一般數值航測區域空中三角平差的方式進行，勢必可以節省時間與人力，國外已有的相關文獻如OEEPE 的研究計畫中，

Paszotta(1999)、Shan(1999)和 Höhle(1999)分別提出不同的更新方法。

考量到不損失精度以及電腦計算速度已逐漸提升，發展了一套方法，以「現有正射影像」與「數值地形模型」為基礎，在個人電腦上進行一般數值航測的作業，利用「蒙地卡羅法」之演算以半自動方式求解單張航攝影像外方位，並可對該變遷地區的正射影像進行更新。

2. 蒙地卡羅法於單張航攝影像外方位求定

在此引入一個「影像距離」－d的概念，兩共軛影像之影像距離d其定義為

r

d = 1 −

(1) 其中：

r

為兩影像之相關係數值

此概念類似於兩點間之實際距離，當兩點越接近則其距離越短，若兩點重合則距離為零。同樣的，若是兩影像區塊越相似，則其影像距離應該要越短，若兩張影像完全相同時，此時其影像距離就為零。

並可以將影像距離函數

d ( C ,' C " )

定義成：

"

1 ,'

)

"

,'

(C C r_C _C

d = − (2) 其中：

r

_C_,'C_"－

C '

與

C "

兩影像之相關係數值

1國立成功大學測量與空間資訊學系碩士

2 國立成功大學測量與空間資訊學系教授

收到日期:民國 92 年 09 月 02 日修改日期:民國 93 年 05 月 21 日接受日期:民國 93 年 05 月 26 日

(2)

'

C

－現有正射影像

"

C

－經正射糾正後的新飛航攝影像若是現有正射影像與新拍航攝影像所製作出來新的正射影像兩者間無任何差異的話，則共軛影像其影像距離函數d(C ,'C")=0。

但所要進行更新的影像有地貌上的改變，所要進行求解的目標函數為：

min )

"

, '

( C C =

d

(3)

而要在大範圍內尋找目標函數的最小值是相當不容易的，所以，將全區域的匹配方式改變成只選取需要更新的範圍之中地貌沒有變化的部分來進行匹配工作。

因影像距離函數為一非線性多變數的不連續函數，所以選定以蒙地卡羅法做為演算法。其主要的概念為，藉由給定不同的外方位初始值、配合 DTM 資料所製作而成的正射影像應該和現有正射影像具最大相關性。以此原則來求解目標函數

)

"

,' ( C C

d

的最小值。

依最小值搜尋演算法，其流程可分成：

(一)初始值給定、(二)亂數範圍、(三)權值給定及(四)中止條件，四個部分說明。(徐業良，1995) (一) 初始值給定

為了有多餘觀測,並求得未知數之中誤差，以四個控制點空中後交會求得新飛航攝影像之外方位參數，做為後續應用蒙地卡羅法計算之外方位初始值。

(二) 亂數範圍

因為所使用的演算法為一隨機方法，所以要預先訂定一亂數範圍，並於此範圍內尋找目標函數的最小值(Baran，1983)：

i

i = t ⋅m

2

ε α (4) 由於外方位僅由四個點位決定，其自由度為 2，在 90%的信心水準下，外方位元素會落在 2.92

倍中誤差之內。

(三) 權值給定

每個影像區塊是人為選定，其區塊之影響力受到其相關係數值高低與區塊範圍大小兩個變量而有所不同：

∑

=

⋅

=

^N

i

d

p d

1

，

p

_i

= s

_i

⋅ r

_i² (5)

其中：

N－區塊總數

pi－第i 個區塊的權值 si－第i 個區塊的大小 ri－第i 個區塊的相關係數值 (四) 中止條件

1. 化量指標V_c

以新舊兩個時期之正射影像進行比較在兩次迭代間其影像區塊灰值之變化程度。影像變化量指標 Vc之定義為：

指標 (6) 其中：

N－影像區塊的總數

ri－第i 個影像區塊的相關係數值

K

－第

K

次的迭代點 2. 移量指標

V

_p

兩次迭代間像元位置之變動量應在容許範圍內。像元位移量指標

V

_p之定義為：

其中：

N－影像區塊的總數

) ,

( px

_ij

py

_ij －第i 個影像區塊其第 j 個角點

K

－第

K

次的迭代點

[ ]

N

py py px

px Vp

N

i j

K ij K ij K

ij K

∑ ∑

ij

= =

−

− + −

−

= ¹

4 1

2 ) 1 ( ) ( 2 ) 1 ( )

( ) ( )

(

N r r V

N

i

K i K i c

∑

=

− −

= ¹

2 ) 1 ( )

( )

(

指標 (7)

(3)

3. 統計量

另外，也可由統計量上預估取樣數多寡(Paszotta，

1999)：

) 1 ( ) 1

( − ε

ⁿ ^M

≤ − γ

(8) 其中：

M

－取樣數(number of sampling) n－變數維度(variable dimension) ε－精確性指標(the measure of precision) γ－機率(probability)

3. 實驗結果與分析

以對局部地區進行正射影像更新為前提之下，本研究使用之實驗資料與檢核資料如下：

(一)現有正射影像

1995 年台南市成功大學 1/5000 正射影像。

(二)數值地形模型資料

DTM 間隔大小為 50mX50m。

(三)新拍航攝影像

同一地區新拍的航照相片，像比例尺約為 1/5000。

(四)檢核資料

在解析製圖儀上，以一般數值航測作業流程解算得外方位參數值並製作正射影像供後續比較，並且立體觀測15 個平面點，做為點位檢核資料。

以上各資料所涵蓋的範圍如圖1 所示：

整體作業的架構上可分成兩個部分進行說明：

(1)資料前置處理(2)蒙地卡羅法主體架構

3.1 資料前置處理

前置處理是將資料進行預處理，其流程如圖 3-2，處理完成後可得到以下的中繼資料：

(1) 外方位初始值及各外方位元素之中誤差

(2) 選的影像區塊及其四個角點的地面坐標 (3) 轉換後的 DTM 格點資料

(4) 仿射轉換參數

圖 2 資料前置處理

3.2 蒙地卡羅法主體架構

圖 1 各資料涵蓋範圍示意圖

Areial Photo Existing orthoimage

Reconstruct orthoimage DTM

(4)

圖 3 蒙地卡羅法主體架構

蒙地卡羅主體架構如圖3，其主要原則為共軛影像區塊應具最大相關值，也就是，當某個取樣點上，造成了影像距離函數值下降時，表示該次的外方位參數值所製成的正射影像區塊，與現有正射影像中的影像區塊，兩者呈現較高的相關係數值，則該次的外方位參數較前次為佳。

程式的流程中再配合門檻值與指標的檢驗來決定是否合乎要求。依本實驗之情況，各門檻值之設定為，影像距離d<0.3*影像區塊總數，影像變化量指標

V

_c<0.5% ，像元位移量指標

V

_p ＜3 pixel；若只單純的以取樣數

M

來決定中止條件，

則迭代次數<147 次。（楊靈峰，2003）

3.3 實驗平台與設備

本研究中所使用之作業系統為Windows XP。

個人電腦的硬體設備為CPU：AMD Athlon Processor 1.60GHz

記憶體：512MB

程式主體以Matlab 6.5 撰寫。

3.4 實驗結果與分析

以其中一次的實驗成果，來說明程式的迭代過程和其兩指標值的變化：

圖 4 目標函數迭代圖(橫軸為迭代次數，縱軸為影像距離)

表 1 迭代次數與指標次數 d 值

V

_C V_P

1 4.126 0 0 3 2.2025 1347.3 0.56721 5 1.8015 10.154 0.00953 11 1.6454 7.6032 0.00500 98 1.6396 2.7467 0.00053

在目標函數迭代圖4 中，下方實線是當函數值降低時被選定視為目標函數的最小值，而上方虛線是表示在這之中的取樣點的函數值高於前次的最小值，屬於無效的計算。

表格1 之指標值為，在 d 值小於門檻值後，若函數值有下降時，則記錄下該次迭代之指標值，並進行判斷，直到符合中止條止後才停止。

3.5 總結

0 10 20 30 40 50 60 70 80 90 100

1.5 2 2.5 3 3.5 4 4.5

2.2025

1.8015

1.6454

1.6396

(5)

以蒙地卡羅方法解算所得40 組外方位參數實驗結果與空三解算之成果作一比較，見表2。

表 2 40 組實驗解算得之外方位統計數據檢核資

料 8839.738 4338.512 1606.995 -2.2588 -2.5799 87.2643 0

X (m) Y0(m) Z0(m)

ω

(deg)

ϕ

(deg)

κ

(deg) Mean 8839.753 4338.385 1606.671 -2.2546 -2.5798 87.2597 STD 0.174 0.212 0.256 0.007 0.007 0.007

0

∆X (m) ∆Y0(m) ∆Z0(m) ∆

ω

(deg) ∆

ϕ

^(deg)

∆ κ

(deg) Max 0.754 0.559 0.867 0.0301 0.0261 0.0214 Min -0.516 -0.896 -0.998 -0.017 -0.021 -0.026

並以解得之外方位最或是值，配合 DTM 資料，以前交會方式求解 15 個檢核點之地面點坐標，可與實際檢核點坐標資料做一比較，圖5 即為

各檢核點之平面位移量(以公尺為單位)。(楊靈峰，

2003)

另外，在全部40 組的實驗結果中，15 個檢核點平面位移量之統計資料，記錄於表3。雖然每組的外方位最或是值有些許的差異，但是由表3 中的

統計數據可看出，整體上而言，是能夠達到1/5000 正射像片圖的精度要求。

1.684 1.686 1.688 1.69 1.692 1.694

x 10⁵ 2.5436

2.5438 2.544 2.5442 2.5444 2.5446 2.5448

2.545x 10⁶

1 0.174

2 0.0667 3 0.0527 4

0.106

5 0.157

6 0.153 7 0.0967

8 0.036

9 0.116

10 0.0238

11 0.239

12 0.117

13 0.122

14 0.188

15 0.252

圖 5 平面點位偏移圖(縱橫坐標以公尺為單位)

(6)

X坐標(m) Y坐標(m) 點

位 Max Min Mean Std. Max Min Mean Std.

1 0.393 -0.002 0.077 0.149 -0.394 -0.010 0.084 0.178 2 0.235 -0.001 0.052 0.10 -0.263 -0.000 0.030 0.110 3 0.209 -0.003 0.047 0.097 0.201 -0.003 0.001 0.099 4 -0.271 -0.004 0.064 0.131 0.298 0.005 0.036 0.115 5 -0.264 -0.000 0.027 0.132 0.363 0.001 0.102 0.161 6 0.257 -0.001 0.037 0.111 -0.410 0.001 0.094 0.168 7 0.203 0.001 0.029 0.089 -0.344 -0.010 0.067 0.130 8 0.128 -0.006 0.014 0.068 -0.271 0.003 0.032 0.100 9 -0.180 -0.001 -0.012 0.080 0.253 0.006 0.061 0.138 10 0.133 0.003 0.006 0.062 0.235 0.001 0.010 0.094 11 -0.314 -0.0078 -0.017 0.173 -0.486 -0.021 0.129 0.193 12 -0.214 -0.012 -0.016 0.104 -0.343 -0.016 0.065 0.120 13 -0.231 -0.004 -0.044 0.079 0.251 -0.011 0.031 0.136 14 -0.295 -0.002 -0.077 0.104 0.328 -0.006 0.066 0.194 15 -0.350 -0.002 -0.108 0.126 0.629 -0.019 0.110 0.267

並再與解析製圖上所製作的正射影像相比較，由圖3-6 中可以看出，色彩變化較顯著的地方

均是在建築物的邊緣部分，原因為使用的DTM 資料並不包含建物的高程，所以仍存在高差移位。

圖 6 正射影像灰值變化圖

表 3 以 40 組實驗結果計算得之檢核點坐標其統計數據

(7)

4. 結論與建議

1. 使用本方法時可以只選擇四個地面控制點來進行空中後交會，若是存在粗差，也可以藉由影像區塊匹配的相似性來識別，偵錯動作較為容易。

2. 經糾正後的正射影像在人造建物或是植被覆蓋的地方，仍會有高差移位存在，所以較適合於裸露地表的更新。但若能取得DSM 資料，

本方法應用於建物密集的都會區也是可行的。

3. 以 40 組實驗數據所得之外方位最或是值以前交會求得15 個檢核點地面坐標並與實際的 15 個檢核點位坐標資料進行比較，在平面坐標位移量上所有點位平均位移在20±10cm以內。

4. 一般數值航測作業所需時間遠超過本法的平均時間，而本法可以達到相同數值製圖之精度要求，加上在個人電腦上執行之方便性，所以相當適合於小範圍地物變遷之正射影像更新。

5. 若是在小比例尺的正射影像中，其高差移位更不明顯。但是地控點精度會受到影響，那麼由後交會所得到的外方位初始值就可由蒙地卡羅法得到良好的改善。

參考文獻

徐業良，「工程最佳化設計」，國立編譯館主編，宏明書局印行，華泰書局總經銷，1995。

楊靈峰，2003，「以現有正射影像與DTM資料解算外方位參數」，國立成功大學測量系碩士論文。

Baran, W., 1983. Teoretyczne podstwy opracowania wynikow pomiarow geodezyjnych. PWN, Warszawa.

Höhle, J., 1996, Experiences with the production of digital orthophotos. Photogrammetric

Engineering & Remote Sensing, 62 (10): pp.

1189-1194.

Paszotta, Z., 1999. Method of Exterior Orientation of Aerial Image by Matching Orthoimages, pp.

9-15. ISBN 83-88343-79-3.

Shan J., 1999.Automtic Exterior Image Orientation with Orthoimage and DTM. OEEPE Publication No 36, pp. 151-157.

(8)

Automatic Exterior Image Orientation Using Existing Orthoimage and DTM

Leng-Feng Yang

¹

Yang-Cheng Liao

²

ABSTRACT

Determination of aerial photo orientation is very fundamental issue in photogrammetry. Along with the developments of computer science and join in image-processing technology, many works in digital photogrammetry can be automatically accomplished.

As for smaller area, we can update changed terrain features by using the existing data not by using AT on analytical plotter that should be saving much time and human resources.

Concerning some partial changes in aerial images, we used Monte Carlo method with area-based matching technique and design a semiautomatic operation procedure which can perform on PC. By determination of aerial image orientation and DTM, we can update the existing orthoimage. And the accuracy of experience result is equal to the accuracy of digital photogrammetry.

Key Words：Orthoimage, DTM, Monte Carlo method, Local Regeneration

Received Date: Sep. 02, 2003 Revised Date: May 21, 2004 Accepted Date: May 26, 2004

1 Graduate Student with M.S. degree

2 Professor Department of Surveying Engineering National Cheng-Kung University

以現有正射影像與 DTM 資料解算外方位參數