循序式模擬退火系統於圖形偵測與震測圖形識別之應用

(1)

國

立

交

通

大

學

多媒體工程研究所

碩

士

論

文

循序式模擬退火系統於

圖形偵測與震測圖形識別之應用

Th

e

S

e

qu

e

nt

i

al

Si

m

ul

a

te

d

A

n

ne

a

li

n

g

S

ys

t

em

fo

r

P

a

tt

e

rn

De

t

ec

t

io

n

a

n

d

S

ei

s

mi

c

A

p

pl

i

ca

t

io

n

s

研

究

生：周映良

指導教授：黃國源

教授

(2)

循序式模擬退火系統於

圖形偵測與震測圖形識別之應用

Th

e

S

e

qu

e

nt

i

al

Si

m

ul

a

te

d

A

n

ne

a

li

n

g

S

ys

t

em

fo

r

P

a

tt

e

rn

De

t

ec

t

io

n

a

n

d

S

ei

s

mi

c

A

p

pl

i

ca

t

io

n

s

研

究

生：周映良 Student：Ying-Liang Chou

指導教授：黃國源

教授

Advisor：Dr.

Kou-Yuan Huang

國

立

交

通

大

學

多媒體工程研究所

碩

士

論

文

A Thesis

Submitted to Institute of Multimedia Engineering

College of Computer Science

National Chiao Tung University

in partial Fulfillment of the Requirements

for the Degree of

Master

in

Computer Science

June 2008

Hsinchu, Taiwan

中華民國九十七年六月

(3)

循序式模擬退火系統於圖形偵測與震測圖形識別之應用研究生：周映良指導教授：黃國源教授國立交通大學多媒體工程研究所碩士班

摘要

Hough transform (HT) 通常用來偵測影像中的直線、圓、橢圓、和雙曲線。但此方法需要大量的計算，且耗費大量記憶體空間。我們提出一個運用模擬退火演算法的循序式偵測系統來偵測影像中的直線、圓、橢圓、和雙曲線之參數。循序式模擬退火圖形偵測系統可搜尋到一組圖形參數，使得影像上的點到這組圖形的距離為最小。在偵測影像中的圖形時，我們不同時偵測所有圖形總類，而是以循序式的觀念，一次偵測一個種類的圖形，並將偵測到的圖形所包含的點移除，再對剩餘的資料重複偵測。我們先將循序式模擬退火圖形偵測系統應用於模擬的影像資料上，偵測模擬影像成功之後，再將之應用於偵測模擬的與真實的單炸點震測圖形上，偵測直接波圖形 (直線) 與反射波圖形 (雙曲線) ，對於瞭解地層特性及後續的研究有很大的幫助。我們也將循序式偵測系統應用於偵測 CDP 圖形中的雙曲線，並利用偵測出來的雙曲線參數，推算出地層中每層波的旅行速度 (interval velocity)，此結果可使速度分析達到自動化的效果。關鍵字：模擬退火、霍夫演算法、震測圖形。

(4)

The Sequential Simulated Annealing System for Pattern

Detection and Seismic Applications

Student：Ying-Liang Chou Advisor：Dr.

Kou-Yuan Huang

Institute of Multimedia Engineering

National Chiao Tung University

Abstract

Usually Hough transform (HT) was used to detect lines, circles, ellipses, and hyperbolas in an image. But it needed huge calculation and consumed the memory.

Using simulated annealing (SA), we propose the sequential detection system to detect the parameters of lines, circle, ellipse, and hyperbola in images. The proposed sequential simulated annealing pattern detection system has the capability of searching a set of parameter vectors with global minimal error. In pattern detection, instead of detection the patterns of all types at one step, we use a sequential concept to detect the patterns of a type at one step. Then we remove the detected patterns. The reaming data will be used in the SA to detect the patterns of other type.

Using sequential simulated annealing pattern detection system to detect simulated images can get good detections. The method is also applied to detect the parameters of direct wave (line) and reflected wave pattern (hyperbola) in simulated and real one-shot seismograms, the results can improve seismic interpretation and further seismic data processing. We also use the system to detect the hyperbolic patterns in the CDP image. The detected parameters can be used to calculate the interval velocity in each layer. The results makes velocity analysis automated.

(5)

誌

謝

本文承蒙黃國源教授之悉心指點、指正引導方能順利完成，在此致上最誠摯的謝意。黃教授不僅教學態度嚴謹，做研究更是實事求是、講求科學精神，是我在求學上以及將來就業上的典範。同時感謝參加口試的口試委員，劉長遠教授、蘇豐文教授、以及王榮華教授給予學生修改論文上的諸多寶貴建議。感謝實驗室的同學及好友，在論文研究期間所給予我精神上莫大的鼓勵；最後，謹以此文獻給我摯愛的雙親。

(6)

中文摘要

………

…

………

i

英文摘要

………

…

………

ii

誌謝

………

…

………

iii

………

…

………

iv

1 介紹 … … …

1

2 模擬退火演算法 ……….

₂

3 循序式的模擬退火圖形偵測系統……….

₅

3.1

_{循序式圖形偵測系統 … … … . .}

₅

3.1.1

_{移除已偵測到的圖形之觀念...}

₅

3.1.2

_{系統架構...}

₅

3.2

_{定義圖形參數…….……….………}

₈

3.3

_{定義系統誤差值.……….………}

₉

3.3.1

_{一個點到一個圖形的距離………..………}

₉

3.3.2

_{一個點到所有圖形的距離………..………}

₉

3.3.2

_{所有點到所有圖形的距離………..………}

₁₀

3.4

_{使用模擬退火演算法尋找圖形參數………..……}

₁₁

3.4.1

_{使用模擬退火演算法尋找直線參數………..……}

₁₂

3.4.2

_{使用模擬退火演算法尋找雙曲線參數…………..……}

₁₄

3.4.3

_{使用模擬退火演算法尋找圓與橢圓參數………….…}

₁₇

3.4.4

_{使用模擬退火演算法尋找圖形參數之範例……….…}

₂₀

3.5

_{循序式震測圖形偵測系統……….………..……}

₂₄

3.5.1

_{系統架構………..………….…}

₂₄

3.5.2

_{使用模擬退火演算法尋找雙曲線的漸進線參數….…}

₂₇

4 實驗結果………

30

4.1

_{模擬參數圖形偵測……….}

₃₀

4.2

_{真實影像中的圖形偵測………}

₆₇

4.3

_{模擬單炸點震測訊號圖形偵測………}

₉₂

4.3.1

_{偵測模擬水平地層單炸點震測訊號圖形………}

₉₂

4.3.2

_{偵測模擬傾斜地層單炸點震測訊號圖形………}

₁₀₀

4.4

_{真實震測訊號圖形偵測……….………}

₁₀₈

5 偵測 C

D

P

g

a

th

e

r之圖形並反推地層原貌.………

117

5.1

_{介紹……….}

₁₁₇

5.2

_{實驗步驟……….}

₁₁₈

5.2.1

_{模擬地層圖………}

₁₁₈

5.2.2

_{T h e C D P ga t he r … … … …}

₁₂₂

5.2.3

_{利用循序式圖形偵測系統做速度分析………}

₁₂₈

5.2.4

_{N M O … … … . . … … … …}

₁₃₆

(7)

5.2.5

_{還原地層原貌………..………}

₁₃₉

5.3

_{使用 VSPEC 做速度分析……….}

₁₄₀

5.4

_{討論……….}

₁₄₄

6 結論與討論 ………

145

6.1

_{結論……….}

₁₄₅

6.2

_{討論……….}

₁₄₆

References

………

149

(8)

圖說明

圖 2.1 描述一個函數解的能量上升與下降方向，虛線箭頭為能量上升方向，實線箭頭為能量下降方向 … … … 3 圖 2.2 模擬退火演算法的流程圖 … … … 4 圖 3.1 一個使用循序式圖形偵測的例子………5 圖 3.2 系統偵測圖形種類的順序………5 圖 3.3 系統架構圖………6 圖 3.4 循序式圖形偵測系統的流程圖………7 圖 3.5 使用模擬退火演算法來尋找直線圖形之參數向量的流程圖………13 圖 3.6 使用模擬退火演算法來尋找雙曲線圖形之參數向量的流程圖 … … … 16 圖 3.7 使用模擬退火演算法來尋找雙曲線圖形之參數向量的流程圖 … … … 19 圖 3.8 調整參數 a 和 b。(a)原始參數， (b)新的參數， (c)不接受新的參數………20 圖 3.9 調整參數 c 。 (a)原始參數， (b)新的參數， (c)接受新的參 數 … … … 21 圖 3.10 調整中心參數(mx, my)。(a)原始參數，(b)新的參數，(c)接受新的參數………21 圖 3.11 調整形狀參數 a 和 b。(a)原始參數，(b)新的參數，(c)接受新的參數………22 圖 3.12 整旋轉角度參數 θ 。 (a)原始參數， (b)新的參數， (c)接受新的參數………23 圖 3.13 調整圖形大小參數 f。(a)原始參數，(b)新的參數，(c)不接受新的參數………24 圖 3.14 偵測震測圖形種類的順序 … … … 24 圖 3.15 循序式震測圖形偵測系統的流程圖 … … … 26 圖 3.16 使用模擬退火演算法來尋找雙曲線之漸進線圖形之參數向量的流程圖………29 圖 4.1 欲偵測的影像 … … … 32 圖 4.2 同步式偵測系統偵測的結果………33 圖 4.3 循序式偵測系統第一次偵測的結果………33 圖 4.4 第一次偵測的 error下降圖 … … … 34 圖 4.5 移除已偵測到的圖形所包含的點………34 圖 4.6 剩餘的資料………35 圖 4.7 循序式偵測系統第二次偵測的結果………35 圖 4.8 第二次偵測的 error下降圖 … … … 36 圖 4.9 移除已偵測到的圖形所包含的點………36 圖 4.10 剩餘的資料………37 圖 4.11 循序式偵測系統第三次偵測的結果………37 圖 4.12 第三次偵測的 error下降圖 … … … 38 圖 4.13 移除已偵測到的圖形所包含的點………38 圖 4.14 最後偵測的結果………39

(9)

圖 4.15 欲偵測的影像 … … … 40 圖 4.16 同步式偵測系統偵測的結果………41 圖 4.17 循序式偵測系統第一次偵測的結果………41 圖 4.18 第一次偵測的 error下降圖 … … … 42 圖 4.19 移除已偵測到的圖形所包含的點………42 圖 4.20 剩餘的資料………43 圖 4.21 循序式偵測系統第二次偵測的結果………43 圖 4.22 第二次偵測的 error下降圖 … … … 44 圖 4.23 移除已偵測到的圖形所包含的點………44 圖 4.24 剩餘的資料………45 圖 4.25 循序式偵測系統第四次偵測的結果………45 圖 4.26 第四次偵測的 error下降圖 … … … 46 圖 4.27 移除已偵測到的圖形所包含的點………46 圖 4.28 最後偵測的結果………47 圖 4.29 欲偵測的影像 … … … 48 圖 4.30 同步式偵測系統偵測的結果………48 圖 4.31 循序式偵測系統第三次偵測的結果………49 圖 4.32 第三次偵測的 error下降圖 … … … 49 圖 4.33 移除已偵測到的圖形所包含的點………50 圖 4.34 剩餘的資料………50 圖 4.35 循序式偵測系統第四次偵測的結果………51 圖 4.36 第四次偵測的 error下降圖 … … … 51 圖 4.37 移除已偵測到的圖形所包含的點………52 圖 4.38 最後偵測的結果………52 圖 4.39 欲偵測的影像 … … … 54 圖 4.40 同步式偵測系統偵測的結果………54 圖 4.41 循序式偵測系統第一次偵測的結果………55 圖 4.42 第一次偵測的 error下降圖 … … … 55 圖 4.43 移除已偵測到的圖形所包含的點………56 圖 4.44 剩餘的資料………56 圖 4.45 循序式偵測系統第二次偵測的結果………57 圖 4.46 第二次偵測的 error下降圖 … … … 57 圖 4.47 移除已偵測到的圖形所包含的點………58 圖 4.48 剩餘的資料………58 圖 4.49 循序式偵測系統第三次偵測的結果………59 圖 4.50 第三次偵測的 error下降圖 … … … 59 圖 4.51 移除已偵測到的圖形所包含的點………60 圖 4.52 剩餘的資料………60 圖 4.53 循序式偵測系統第五次偵測的結果………61 圖 4.54 第五次偵測的 error下降圖 … … … 61 圖 4.55 移除已偵測到的圖形所包含的點………62

(10)

圖 4.60 剩餘的資料………64 圖 4.61 循序式偵測系統第九次偵測的結果………65 圖 4.62 第九次偵測的 error下降圖 … … … 65 圖 4.63 移除已偵測到的圖形所包含的點………66 圖 4.64 最後偵測的結果………66 圖 4.65 硬幣圖………68 圖 4.66 RGB轉 graylevelvalue………68 圖 4.67 取 threshold把影像變成 binarycode………69 圖 4.68 Sobel取邊界點………69 圖 4.69 輸入的資料………70 圖 4.70 同步式偵測系統偵測的結果………70 圖 4.71 循序式偵測系統第一次偵測的結果………71 圖 4.72 第一次偵測的 error下降圖 … … … 71 圖 4.73 移除已偵測到的圖形所包含的點………72 圖 4.74 最後偵測的結果………72 圖 4.75 真實硬幣偵測結果………73 圖 4.76 真實壁磚圖 … … … 74 圖 4.77 RGB 轉 graycode………75 圖 4.78 使用 imfill將壁磚內部雜訊去除………75 圖 4.79 Sobel取邊界點………76 圖 4.80 輸入的資料 … … … 76 圖 4.81 同步式偵測系統偵測的結果………77 圖 4.82 循序式偵測系統第一次偵測的結果………77 圖 4.83 第一次偵測的 error下降圖 … … … 78 圖 4.84 移除已偵測到的圖形所包含的點………78 圖 4.85 剩餘的資料………79 圖 4.86 循序式偵測系統第二次偵測的結果………79 圖 4.87 第二次偵測的 error下降圖 … … … 80 圖 4.88 移除已偵測到的圖形所包含的點………80 圖 4.89 剩餘的資料………81 圖 4.90 循序式偵測系統第三次偵測的結果………81 圖 4.91 第三次偵測的 error下降圖 … … … 82 圖 4.92 移除已偵測到的圖形所包含的點………82 圖 4.93 剩餘的資料………83 圖 4.94 循序式偵測系統第四次偵測的結果………83 圖 4.95 第四次偵測的 error下降圖 … … … 84 圖 4.96 移除已偵測到的圖形所包含的點………84 圖 4.97 剩餘的資料………85 圖 4.98 循序式偵測系統第六次偵測的結果………85 圖 4.99 第六次偵測的 error下降圖 … … … 86 圖 4.100 移除已偵測到的圖形所包含的點………86 圖 4.101 剩餘的資料………87 圖 4.102 循序式偵測系統第八次偵測的結果………87 圖 4.103 第八次偵測的 error下降圖 … … … 88 圖 4.104 移除已偵測到的圖形所包含的點………88

(11)

圖 4.105 剩餘的資料………89 圖 4.106 循序式偵測系統第十一次偵測的結果………89 圖 4.107 第十一次偵測的 error下降圖………90 圖 4.108 移除已偵測到的圖形所包含的點………90 圖 4.109 最後的偵測結果………91 圖 4.110 真實壁磚圖的偵測結果………91 圖 4.111 水平地層構造圖………92 圖 4.112 模擬水平地層的單炸點震策訊號圖………94 圖 4.113 模擬單炸點震策訊號經過 Envelop處理………94 圖 4.114 經過 threshold後的模擬水平地層震測訊號………95 圖 4.115 經過 peak之後的震測訊號點………95 圖 4.116 循序式偵測系統第一次偵測的結果………96 圖 4.117 第一次偵測的 error下降圖 … … … 96 圖 4.118 移除已偵測到的圖形所包含的點………97 圖 4.119 剩餘的資料………97 圖 4.120 循序式偵測系統第二次偵測的結果………98 圖 4.121 第二次偵測的 error下降圖 … … … 98 圖 4.122 移除已偵測到的圖形所包含的點………99 圖 4.123 最後偵測的結果………99 圖 4.124 傾斜地層構造圖………100 圖 4.125 模擬傾斜地層的單炸點震策訊號圖………102 圖 4.126 模擬單炸點震策訊號經過 Envelop處理………102 圖 4.127 經過 threshold後的模擬傾斜地層震測訊號………103 圖 4.128 經過 peak之後的震測訊號點………103 圖 4.129 循序式偵測系統第一次偵測的結果………104 圖 4.130 第一次偵測的 error下降圖… … … …… …… … … … …… … … … 104 圖 4.131 移除已偵測到的圖形所包含的點………105 圖 4.132 剩餘的資料………105 圖 4.133 循序式偵測系統第二次偵測的結果………106 圖 4.134 第二次偵測的 error下降圖… … … …… …… … … … …… … … … 106 圖 4.135 移除已偵測到的圖形所包含的點………107 圖 4.136 最後偵測的結果………107 圖 4.137 真實地層的單炸點震策訊號圖………109 圖 4.138 真實地層單炸點震策訊號經過 Envelop處理………110 圖 4.139 經過 threshold後的真實地層震測訊號………110 圖 4.140 取 peak之後的震測訊號點………111 圖 4.141 去除 1.7秒以下的訊號後的真實地層震測訊號(欲偵測的影像)…111 圖 4.142 循序式偵測系統第一次偵測的結果………112 圖 4.143 第一次偵測的 error下降圖… … … …… …… … … … …… … … … 112 圖 4.144 移除已偵測到的圖形所包含的點………113 圖 4.145 剩餘的資料………113

(12)

圖 4.150 把偵測出的圖形畫回真實震測圖上………116 圖 5.1 傳統還原地層原貌之流程圖………117 圖 5.2 模擬地層圖………118 圖 5.3 還原地層的範圍………119 圖 5.4 炸點在 2km的 one-shotseismogram………120 圖 5.5 炸點在 3km的 one-shotseismogram………121 圖 5.6 CDP gather訊號來源………122 圖 5.7 每張 CDP gather所含的 trace數目………122 圖 5.8 一個收集 CDP gather的簡單例子………123 圖 5.9 每張 CDP gather所含的 trace數目………124 圖 5.10 第 6 0 張 C D P g a t h e r ，共同中點為 2737.5m……..………126 圖 5.11 第 7 1 張 C D P g a t h e r ，共同中點為 3012.5m………127 圖 5.12 經過 thresholding後的訊號點………128 圖 5.13 經過 peak後的訊號點………129 圖 5.14 (a)第 1 次偵測的結果， (b)被移除的點與圖形………130 圖 5.15 (a)第一次偵測 error 的下降圖， (b)剩餘資料點………131 圖 5.16 (a)第 2次偵測的結果，(b)被移除的點與圖形(c)第二次偵測 error的下降圖 … … … … …… … … …… … … …… … … …… … … …… … … 132 圖 5.17 最後偵測出來的結果 …… … … …… … … …… … … …… … … 133 圖 5.18 每兩條 trace間的偏移量圖………135 圖 5.19 NMO stretching………137 圖 5.20 Stacked示意圖………137 圖 5.21 (a)第 71張 CDP gather做 NMO 修正後的結果，(b)stacked後的結果………138 圖 5.22 還原地層原貌………139 圖 5.23 (a) CDP 71 的 velocity spectra， (b)選取到的速度曲線

圖 … … … … …… … … …… … … …… … … …… … … …… … … 141 圖 5.24 (a)對 CDP 71 做 NM O 之後的圖， (b)stack 之後的圖 … … … … …… … … …… … … …… … … …… … … …… … … 142 圖 5.25 使用 VSPEC 做速度分析後還原的地層圖………143 圖 6.1 決定影像中圖形的數目 K。(a) K =1，(b)K =2，(c)K =3，(d)K = 4， (e) K = 5， (f) K = 1~5 與偵測誤差值的對應 圖 … … … … …… … … …… … … …… … … …… … … …… … … 147

(13)

表說明

表 3.1 a、b、f 的正負關係對應圖形的形狀………8 表 4.1 參數代表意思的對照………30 表 4.2 初始參數設定值………31 表 4.3 初始參數設定值………31 表 4.4 初始參數設定值………31 表 4.5 初始參數設定值………32 表 4.6 初始參數設定值………39 表 4.7 初始參數設定值………40 表 4.8 初始參數設定值………47 表 4.9 初始參數設定值………48 表 4.10 初始參數設定值………53 表 4.11 初始參數設定值………53 表 4.12 初始參數設定值………67 表 4.13 初始參數設定值………67 表 4.14 初始參數設定值………73 表 4.15 初始參數設定值………74 表 4.16 初始參數設定值………93 表 4.17 初始參數設定值………93 表 4.18 偵測出來的各參數值………100 表 4.19 初始參數設定值………101 表 4.20 初始參數設定值………101 表 4.21 偵測出來的各參數值………108 表 4.22 初始參數設定值………108 表 4.23 初始參數設定值………109 表 4.24 偵測出來的各參數值………116 表 5.1 每一層地層波的旅行速度………118 表 5.2 cdp#與 midpoint對照表………123 表 5.3 前 10張 one-shotseismograms中其 trace對應每張 CDPgathe與炸點和接收站位置的關係 … … … 124

表 5.4 第 60張 CDP gather中每條 trace對應炸點與接收站的位置關係 … … … 125 表 5.5 初始參數設定值 … … …… … … …… … … …… … … …… … … 129 表 5.6 CDP gather中雙曲線的參數………134 表 5.7 經換算後之每層波的旅行速度………136 表 5.8 選取到的速度………140 表 6.1 圖 6.1 中的偵測誤差值………148 表 6.2 圖 6.2 中的 CPU time………148

(14)

1 .

介紹

影像中的直線、圓、橢圓、與雙曲線常代表著一些重要的圖形，例如磁磚的接縫(直線)、硬幣(圓)、震測訊號圖形中的反射波(雙曲線)、…等其他圖形。若能精確的偵測出這些圖形，對於電腦視覺方面的應用會有很大的幫助。在圖形識別上，Houghtransform (HT)常用來偵測直線、圓、與橢圓 [1]-[2]。其方法為將影像上的點轉換到參數空間上的座標點並取參數空間上值的 peak(最大值)，然後用找到的參數還原圖形。但 HT需要龐大的計算且耗費記憶體空間。震測訊號圖形識別對於石油探勘上有很大的幫助，其中單炸點震測訊號圖包含了直線的直接波與雙曲線的反射波 [3]-[6]。在 1985年時，Huangetal.提出以 HT 偵測單炸點震測訊號圖中直線的直接波與雙曲線的反射波 [7]。但仍存在參數空間上的 peak選擇和需要大量記憶體空間的問題。 Basak與 Das在 2002年提出 Houghtransform neuralnetwork(HTNN)以類神經網路調整梯度的方法來找出圖形的參數，並用在偵測直線、圓、橢圓的圖形上[8]。然而這方法並沒有用於偵測雙曲線上。在 2006年時，Huangetal.利用 HTNN偵測在單炸點震測訊號圖中直線的直接波與雙曲線的反射波 [9]。然而在解複雜的問題時，HTNN會有落入局部最佳解的缺點。 Kirkpatricketal.在 1983年提出模擬退火演算法 simulatedannealing(SA)[10]。其內容主要重點是利用 Metropoliscriterion，有條件的接受較高能量的狀態，所以可以達到全域最佳解 [11]。SA 可用來解決最佳化的問題，如旅行銷售員問題、電路設計、與分類等問題 [10],[12]-[17]。然而 SA並沒有用於偵測雙曲線上。在 2007年，Chen與 Huang提出利用模擬退火演算法於偵測直線、圓、橢圓、與雙曲線。但其不能同時偵測橢圓和雙曲線 [18]。這裡我們提出利用 SA 的循序式圖形偵測系統，利用循序式分群的觀念，不必一次偵測出所有圖形，而是以每次偵測部分圖形，並將偵測到的圖形所包含的點去除，使得影像中的資料減少來達到正確的偵測結果。我們先將循序式圖形偵測系統應用於模擬的影像資料上，偵測模擬影像之後，再將之應用於偵測模擬的與真實的單炸點震測訊號圖形上，偵測直線的直接波與雙曲線的反射波，對於了解地層特性及後續的研究有很大的幫助。我們也將循序式偵測系統應用於偵測 CDP 圖形中的雙曲線，並利用偵測出來的雙曲線參數，推算出地層中每層波的旅行速度(intervalvelocity)，此結果可使速度分析達到自動化的效果。

(15)

2 .

模擬退火演算法

Kirkpatricketal.在 1983年提出模擬退火演算法，並用於解積體電路迴圈線路設計與旅行銷售員的問題上 [10]。模擬退火演算法是由觀察物理上玻璃結晶的技術而來。物理上的玻璃結晶必須要從很高的溫度開始，經過慢慢的降溫，使粒子達到穩定的狀態，才能得到良好的結晶結果。同樣的，利用模擬退火演算法在解最佳化的問題時，也是從一個夠高的溫度開始降溫，並使用 Metropolis criterion來決定是否接受粒子移動到較高的能量上 [11]。當系統溫度降到夠低 時，即可得到一個最佳解。圖 2.1 描述一個函式 F(x)的解 x1其能量上升與下降 方向，使用模擬退火演算法可以有條件的接受函式 F(x)的解 x1其能量往上升 的方向移動，可避免局部最小值，進而找到全域最小值的解。圖 2.2 為模擬退 火演算法的流程圖。下面為模擬退火演算法的主要執行步驟。 Algorithm:使用模擬退火演算法找到問題的最佳解。 Input:欲解問題的函式表示法。 Output:此問題的函式最佳解。 Step1:初始值設定 1.1. 設定初始溫度。 1.2. 設定函式解在問題空間中的初始位置。 1.3. 定義溫度下降的公式。 1.4. 定義能量公式。 Step2: 2.1. 給函式解在問題空間中的位置隨機的移動，獲得一個新的位置。 2.2. 使用能量公式算出函式解在新位置能量 E’與原位置能量 E，其相差為 E’ - E = △E。 2.3. 使用 Metropoliscriterion判斷是否接受新的位置: (a)如果△E ≤ 0，則接受新的位置。

(b)如果△E > 0，則計算接受機率 prob = exp[-△E/T]，T代表當時的溫度。在 [0,1]區間內產生一個 uniform random number r，若 prob > r，則接受新的位置，否則保留函式解在原位置。

2.4. 重複 step2，直到達到在此溫度下的平衡狀態為止。

Step3:當平衡狀態達到時，使用溫度下降公式降溫，並跳回 step2執行，直到溫度下降到小於設定的停止條件時，才停止運算。

(16)

(a) 圖 2.1. 描述一個函數解的能量上升與下降方向，虛線箭頭為能量上升方向，實線箭頭為能量下降方向。 x1 F(x) x

(17)

Yes 開始設定初始溫度和函式解在問題空間中的初始位置，定義溫度下降公式與能量公式。給函式解在問題空間中的位置隨機的移動，獲得一個新的位置。 使用能量公式算出函式解在新位置能量 E’ 與原位置能量 E，其相差為 E’- E=△E。

△E ≤ 0

計算 prob = exp[-△E/T]，並產生 一個 uniform random number r。

prob> r 是否達到平衡狀態? 是否達到停止條件? 停止使用溫度下降公式降溫接受新的位置。 Yes No Yes Yes No No

(18)

3 .

循序式的模擬退火圖形偵測系統

3.1. 循序式圖形偵測系統我們利用循序式的圖形偵測觀念，提出循序式圖形偵測系統。影像中包含直線、圓、橢圓、和雙曲線，循序式圖形偵測系統每次偵測出影像中一個種類的圖形，直到找出所有圖形為止。此方法可用於影像中大量的圖形偵測。 3.1.1. 移除已偵測到的圖形之觀念若影像中包含大量的圖形，很困難一次就把所有圖形都偵測出來。所以我們使用循序式的圖形偵測方法，來取代一次就把所有圖形偵測出來的方法。此觀念類似循序式分群演算法的觀念 [20]-[22]。圖 3.1 為一個使用循序式圖形偵 測的例子。使用循序式圖形偵測系統可以將影像中的直線分為一類，雙曲線分為一類，圓與橢圓分為一類。圖 3.1. 一個使用循序式圖形偵測的例子。 3.1.2. 系統架構若影像中包含直線、圓、橢圓、和雙曲線，我們可以使用循序式的模擬退火圖形偵測系統一次偵測一個種類的圖形，在此我們設定偵測圖形種類的順序為先偵測直線，接著偵測雙曲線，最後偵測橢圓(包含圓)。圖 3.2 為系統偵 測圖形種類的順序。圖 3.2. 系統偵測圖形種類的順序。輸入資料偵測直線剩餘資料偵測雙曲線剩餘資料偵測橢圓 (包含圓) Image patterns: line 1, line 2, line 3, hyperbola 4, hyperbola 5,ellipse 6, ellipse 7, ellipse 8.

Line class Hyperbola class Ellipse class

(19)

圖 3.3為我們提出的偵測系統的架構。此系統須先輸入欲偵測的資料，接 著使用模擬退火演算法針對某一個種類的圖形找出其圖形參數值，然後把偵測到的圖形資料去除，再把剩下的資料重新輸入系統，直到所有圖形種類都被偵測過為止。圖 3.3. 系統架構圖。循序式的模擬退火圖形偵測系統主要包含三個部分:1.定義系統誤差值的計算方式，2.使用模擬退火演算法找出使系統的誤差值最小的圖形參數值，3. 移除偵測到的圖形資料。為了得到系統的誤差值，首先我們定義用來表示圖形的參數，接著計算所有點到所有圖形的距離來定義系統誤差值。其中定義橢圓和雙曲線圖形參數與系統誤差值的方式參考陳凱儒於 96年發表的碩士論文 [19]。 圖 3.4 為循序式圖形偵測系統的流程圖。首先，我們依序偵測影像中的直 線、雙曲線、和橢圓(包含圓)，並設定圖形參數初始值。 接著，輸入影像中的所有點座標和影像中的圖形個數 K。若一個圖形包含 夠多靠近此圖形的點，我們則稱此圖形為已被偵測出來的圖形。所以我們必須設定點與圖形的距離必須多小(點與圖形的距離必須小於設定的門檻距離)，且圖形須包含多少靠此圖形的點(包含的點數必須大於設定的門檻點數)，其中門 檻點數的算法設定為所有點數 N，除以影像中圖形的個數 K。然而有些影像中 會有雜訊點，所以影像中的每個圖形不一定都包含 N 除以 K 個點，所以我們設 定一個偏移值，使得門檻點數降低，門檻點數的算法如下: 門檻點數=_N K/ _−偏移值 (1) 再來使用 SA偵測圖形的參數，若要偵測的圖形是直線，則使用偵測直線參數的 SA 演算法。若要偵測的圖形是雙曲線，則使用偵測雙曲線參數的 SA 演算法。若要偵測的圖形是圓與橢圓，則使用偵測圓與橢圓參數的 SA演算法。接著將已偵測到的圖形所包含的點去除，剩下的點再使用 SA偵測圖形的參數。重複此步驟直到達到我們設定的執行次數時，則換偵測下一個種類的圖形。直 Input N data

Simulated annealing for parameter detection Detected parameters Remove the detected data Remaining data

(20)

圖 3.4. 循序式圖形偵測系統的流程圖。 移除已偵測到的 m 個圖形所包含的點。圖形個數變為 K = K-m。 Yes No No No Yes 設定偵測影像中的圖形種類:直線、雙曲線、和橢圓(包含圓)。設定圖形參數初始值。 輸入影像中的所有點座標和影像中的圖形個數 K。 若設定偵測直線，則使用偵測直線圖形參數的 SA。若設定偵測雙曲線，則使用偵測雙曲線圖形參數的 SA。若設定偵測橢圓，則使用偵測橢圓圖形參數的 SA。點與圖形的距離是否小於門檻距離，且圖形所包含的點數大於門檻點數? 剩餘的點座標。是否達到設定的執行次數? 是否所有形狀的圖形都已偵測過? 停止 Yes 設定門檻距離和偏移值，計算門檻點數。

(21)

3.2. 定義圖形參數我們使用直線的公式: ax by c+ = (2) 來偵測直線。利用參數向量 lp = [a, b, c]T _{來表示一條直線圖形，例如第}_k_個直 線圖形表示為 lpk = [ak, bk, ck]T，若總共有K個直線圖形，則使用一個矩陣 LP = [lp1, lp2, …, lpK] 來表示。我們使用可表示旋轉與平移的雙曲線和橢圓的一般式: 2 2

(( x) cos ( y) sin ) ( ( x) sin ( y) cos )

a x−m θ + y−m θ + − −b x m θ + y−m θ = f (3) 表 3.1 為 a、b、f 的正負關係對應圖形的形狀。當 a、b、f 同號時，圖形為 圓或橢圓，a、b 異號時則為雙曲線，a、b 異號且 f = 0 時為雙曲線的漸進線。 在偵測雙曲線時，我們設定 f 永遠大於 0，用來避免雙曲線的漸進線與直線混 淆。我們利用式(2)中的參數 mx、my、a、b、θ 和 f 來表示任何旋轉與平移的橢 圓和雙曲線。參數 (mx,my) 表示圖形的中心座標，θ 為旋轉的角度，f 為圖形的 大小，由 a、b、f 的正負關係可得圖形的形狀，如表 3.1所示。我們使用參數 向量 p = [mx, my, a, b, θ, f]T 來表示一個圖形，例如第k個圖形表示為 pk = [mk,x, m k,y, ak, bk, θk, fk]T，若總共有K個圖形，則使用一個矩陣 P = [p1, p2, …, pK] 來表示。 表 3.1. a、b、f 的正負關係對應圖形的形狀。 a b f Graph ＋＋＋ Ellipse －－－ _Ellipse ＋－＋ Hyperbola －＋＋ Hyperbola ＋－－ Hyperbola －＋－ _Hyperbola ＋－ 0 Asymptote －＋ ₀ _Asymptote －－＋ No graph

(22)

若有 k1個直線，k2個橢圓，k3個雙曲線，總共 k1+ k2+ k3= K 個圖形在影像 中，若要表示這 K 個圖形，則使用一個矩陣 AP = [LP, P]= [lp1, lp2,…, lp k1, p k1+1, …, p k1+ k2, p k1+ k2+1, …,pK]。 3.3. 定義系統誤差值為了定義系統的誤差值，首先我們討論一個點到一個圖形的距離，接著定義一個點到K個圖形的距離為這一個點到K個圖形的距離中的最小距離。最後定系統的誤差值為N個點到K個圖形的距離的平均值。 3.3.1. 一個點到一個圖形的距離 我們定義一個座標點 xi到第 k 個直線圖形的距離為: 2 2 ( , ) k i k i k k i i k k a x b y c d x y a b + − = + (4) 而對於一個座標點 xi到第 k 個橢圓或雙曲線圖形的距離我們則定義為: 2 , , 2 , , ( , ) | (( ) cos ( ) sin ) ( ( ) sin ( ) cos ) | k i i k i k x k i k y k k i k x k i k y k k d x y a x m y m b x m y m f θ θ θ θ = − + − + − − + − − (5) 我們想找出一組參數 mx、my、a、b、θ 和 f，使得 dk為最小。當 a = 0, b = 0, 和 f = 0 時，d x y_k( ,_i _i)=0，然而這並不是我們想要的答案。為了解決這個問題， 我們對參數 a 和 b 作正規化，對式(4)中的 a 和 b 除以 ab ，改變後的 a 和 b 關係變為|ab| 1= ，即可避免 a = 0, b = 0, 和 f = 0 的情形發生。 3.3.2. 一個點到所有圖形的距離 若影像內有 K 個圖形，我們定義一個點 xi到這 K 個圖形的距離為: E( ) minx_i =

(

d₁( ), ( ), , ( ), ,x_i d₂ x_i Kd_k x_i Kd_K( )x_i

)

(6) 在偵測直線時，d x 為式(_k( )_i 4)。在偵測橢圓或雙曲線時，d x 皆為式(_k( )_i 5)。

(23)

3.3.3. 所有點到所有圖形的距離 若有 k1個直線(N1個點)，k2 個雙曲線(N2個點)，k3 個橢圓，(N3個點)總共 k1+k2+k3 = K 個圖形在影像中(N1+N2+N3 = N 為影像中的點數)。首先在偵測直線 時，設定 LP = [lp1, lp2, …, lpK]，此時系統誤差值為: 1 1 (LP) ( ) N i i i E E N = =

∑

x (7a) 移除已偵測到的 k1個直線，影像中圖形個數變為 K= k2+ k3，接著偵測雙曲線， 設定 P = [p1, p2, …, p k2+k3]，此時系統誤差值為: 2 3 1 1 (P) ( ) N i i i E E N N = = +

∑

x (7b) 移除已偵測到的 k2個雙曲線，影像中圖形個數變為 K= k3，接著偵測橢圓，設 定 P = [p1, p2, …, p k3]，此時系統誤差值為: 3 1 1 (P) ( ) N i i i E E N = =

∑

x (7c) 最後將所有偵測到的圖形參數合併為 AP =[LP, P]= [lp1, lp2,…, lpk1, p k1+1, …, p k1+k2, p k1+ k2+1, …,pK]，我們定義系統的總誤差值如下: 1 1 (AP) ( ) N i i i E E N = =

∑

x (7) 3.4. 使用模擬退火演算法尋找參數我們使用模擬退火演算法分別針對直線、雙曲線、和橢圓(包含圓)找出其參數，使的系統的誤差值可以達到全域最小值。其共有三個部分，第一部分為使用模擬退火演算法找出直線的參數，第二部分為使用模擬退火演算法找出雙曲線的參數，第三部分為使用模擬退火演算法找出橢圓(包含圓)的參數。每一部分皆使用下列的降溫公式 [10]:

(24)

3.4.1. 使用模擬退火演算法尋找直線參數 在使用模擬退火演算法尋找直線參數中，同時調整參數 a、b 和 c 並不能 使系統有效的收斂，所以我們使用 2個步驟來調整參數 [10]。首先調整圖形的參數 a 和 b，然後調整參數 c。 Algorithm: 使用SA 演算法偵測影像中直線圖形之參數向量。 Input: 所有在影像中的點 (N 個點)，設定K為圖形的數目。 Output: K 個直線參數向量。 Step1:初始值設定 1.1. 在 t= 1時，給定一夠高的溫度，即 T(1) = Tmax。 1.2. 使用式(8)為降溫公式。 1.3. 產生直線參數向量 lp1, lp2, ..., lpk, …, lpK, 其中 lpk= [ak, bk, ck]T，一個 lp 表示為一個直線圖形，設定 LP = [lp1, lp2, ..., lpk, …, lpK]。 1.4. 使用式(4)、式(6)和、式(7a)計算 E(LP)。 Step2:隨機改變參數向量並在每一個溫度中決定新的參數向量。 For m = 1 to Nt (在每個溫度中嘗試 Nt 次) For k = 1 to K (K 為影像中的圖形個數) (a) 隨機改變第 k 個圖形的參數 a 和 b: [a'_k b'_k]T =[a b_k ]_k T +α_abn (9) 其中 n = [n1 n2]T是一個 2 × 1 的隨機向量，n1 和 n2是高斯隨機變數 在 N(0, 1)範圍內，αab為一常數。改變後 lp’k = [a’k, b’k, ck]T， LP′ = [lp1, lp2, …, lp′k,…, lpK]。使用式(4)、式(6)和、式(7a)計算 E(LP′)。使用

Metropoliscriterion判斷是否接受 LP′。若∆E = E(LP′) - E(LP) ≤ 0，則 接受 LP′。若∆E = E(LP′) - E(LP) > 0，計算接受機率 prob =

exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r， 則接受 LP′，否則維持原本的 LP。 (b) 隨機改變第 k 個圖形的參數 c: c'k = ck+αcn (10) 其中 n 是一個高斯隨機變數在 N(0, 1)範圍內，αc為一常數。改變後 lp’k = [ak, bk, c’k]T， LP′ = [lp1, lp2, …, lp′k,…, lpK]。使用式(4)、式(6) 和、式(7a)計算 E(P′)。使用 Metropoliscriterion判斷是否接受 LP′。若

∆E = E(LP′) - E(LP) ≤ 0，則接受 LP′。若∆E = E(LP′) - E(LP) > 0，

(25)

機變數 r。若 prob > r，則接受 LP′，否則維持原本的 LP。

End for k End for m

Step 3: 對系統降溫

(26)

Yes 開始 設定 K 個參數向量: LP = [lp1, lp2, …, lpk, …, lpK]，其中 lpk = [ak′ bk′ ck] T_。 設定初始溫度 Tmax. 定義溫度下降公式。 隨機改變第 k 個參數向量， k = 1, 2, …, K。 使用降溫公式來降溫。所有圖形之參數向量皆被隨機改變? 完成一次循環。 循環次數達到設定的總循環次數 Nt? 停止 No Yes No 隨機改變第k個圖形之 a 和 b: lpk′ = [ak′ bk′ ck]T, LP′ = [lp1, lp2, …, lpk′, …, lpK]。限制 ak′ < 0，bk′ > 0。

計算∆E=E(LP′)-E(LP), 使用 Metropolis criterion 決定是否接受 LP′。

隨機改變第k個圖形之 c: lpk′ = [ak bk ck′]T, LP′ = [lp1, lp2, …, lpk′, …,

lpK]。

計算∆E=E(LP′)-E(LP), 使用 Metropolis criterion 決定是否接受 LP′。 開始一循環。

降溫次數達到設定的總降溫次數? No

Yes

(27)

3.4.2. 使用模擬退火演算法尋找雙曲線參數 在使用模擬退火演算法尋找雙曲線參數中，同時調整參數 mx、my、a、b、 θ 和 f 並不能使系統有效的收斂，所以我們使用 4個步驟來調整參數 [10]。首先調整圖形的中心參數 (mx, my)，接著調整參數 a 和 b，然後調整圖形旋轉角 度參數 θ，最後調整圖形大小參數 f。 雙曲線圖形中的參數 a、b 和 f 同時改變正負號，並不會改變原本的形狀， 如表 3.1 所示。所以我們在演算法中固定 f >0，以減少搜尋的範圍。其演算法詳細過程如下: Algorithm: 使用SA 演算法偵測影像中雙曲線圖形之參數向量， Input: 所有在影像中的點 (N 個點)，設定K為圖形的數目。 Output: K 個雙曲線參數向量。 Step1:初始值設定 1.5. 在 t= 1時，給定一夠高的溫度，即 T(1) = Tmax。 1.6. 使用式(8)為降溫公式。 1.7. 產生參數向量 p1, p2, ..., pk, …, pK, 其中 pk = [mk,x, m k,y, ak, bk, θk, fk]T，一 個 p 表示為一個圖形，設定 P = [p1, p2, ..., pk, …, pK]。 1.8. 使用式(5)、式(6)和、式(7b)計算 E(P)。 Step2:隨機改變參數向量並在每一個溫度中決定新的參數向量。 For m = 1 to Nt (在每個溫度中嘗試 Nt 次) For k = 1 to K (K 為影像中的圖形個數) (a) 隨機改變第 k 個圖形的中心: mn T y k x k T k,y x k m' m m m' ] =[ ] +α [ _, _, _, (11) 其中 n = [n1 n2]T是一個 2 × 1 的隨機向量，n1 和 n2是高斯隨機變數在 N(0, 1)範圍內，αm為一常數。改變後 p’k = [m’k,x, m’k,y, ak, bk, θk, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式(6)和、式(7b)計算 E(P′)。使用

Metropoliscriterion判斷是否接受 P′。若∆E = E(P′) - E(P) ≤ 0，則接受 P′。若∆E = E(P′) - E(P) > 0，計算接受機率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持

(28)

[a_k 'b_k']T =[a_k b_k]T +α_abn (12) 且使用 ' ' |a b_k _k |作正規化。其中 n = [n1 n2]T是一個 2 × 1 的隨機向量， n1 和 n2是高斯隨機變數在 N(0, 1)範圍內，αab為一常數。限制 a’k < 0 和 b’k > 0 來表示雙曲線。改變後 p’k = [mk,x, mk,y, a’k, b’k, θk, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式(6)和、式(7b)計算 E(P′)。使用 Metropolis

criterion判斷是否接受 P′。若∆E = E(P′) - E(P) ≤ 0，則接受 P′。若∆E =

E(P′) - E(P) > 0，計算接受機率 prob = exp[-∆E/T(t)]，並產生一平均分

布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持原本的 P。 (c) 隨機改變第 k 個圖形的旋轉角度: θ'_k=θ_k +α_θn (13) 其中 n 是一個高斯隨機變數在 N(0, 1)範圍內，αθ為一常數。改變後 p’k = [mk,x, mk,y, ak, bk, θ’k, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式 (6)和、式(7b)計算 E(P′)。使用 Metropoliscriterion判斷是否接受 P′。若

∆E = E(P′) - E(P) ≤ 0，則接受 P′。若∆E = E(P′) - E(P) > 0，計算接受機

率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持原本的 P。 (d) 隨機改變第 k 個圖形的中心: f'_k =| f_k +α_fn| (14) 其中 n 是一個高斯隨機變數在 N(0, 1)範圍內，αf為一常數。改變後 p’k = [mk,x, mk,y, ak, bk, θk, f’k]T， P′ = [p1, p2, …, p′k,…, pK]。限制 f >0。使用式(5)、式(6)和、式(7b)計算 E(P′)。使用 Metropoliscriterion判斷是否 接受 P′。若∆E = E(P′) - E(P) ≤ 0，則接受 P′。若∆E = E(P′) - E(P) > 0， 計算接受機率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變 數 r。若 prob > r，則接受 P′，否則維持原本的 P。 End for k End for m Step 3: 對系統降溫使用式(8)對系統降溫，並重複 Step 2 直到達到夠低的溫度為止。

(29)

Yes 開始 設定 K 個參數向量: P = [p1, p2, …, pk, …,pK]，其中 pk = [mk,x mk,y ak bk θk fk]T。 設定初始溫度 Tmax. 定義溫度下降公式。 隨機改變第 k 個參數向量， k = 1, 2, …, K。 使用降溫公式來降溫。所有圖形之參數向量皆被隨機改變? 完成一次循環。 循環次數達到設定的總循環次數 Nt? 停止 No Yes No 隨機改變第k個圖形之 mx 和 my: pk′ = [mk,x′ mk,y′ ak bk θk fk]T, P′ = [ p1, p2, …, pk′, …,pK]。

計算∆E=E(P′)-E(P), 使用 Metropolis criterion 決定是否接受 P′。

隨機改變第k個圖形之 a 和 b: pk′ = [mk,x mk,y ak′ bk′ θk fk]T, P′ = [ p1,

p2, …, pk′, …,pK]。限制 ak′ < 0，bk′ > 0。

隨機改變第k個圖形之 θ: pk′ = [mk,x mk,y ak bk θk′ fk] T

, P′ = [ p1, p2, …,

pk′, …,pK]。

隨機改變第k個圖形之 f: pk′ = [mk,x mk,y ak bk θk fk′]T, P′ = [ p1, p2, …,

pk′, …,pK]。

計算∆E=E(P′)-E(P), 使用 Metropolis criterion 決定是否接受 P′。 開始一循環。

降溫次數達到設定的總降溫次數?

No Yes

(30)

3.4.3. 使用模擬退火演算法尋找圓與橢圓參數 在使用模擬退火演算法尋找圓與橢圓參數中，同時調整參數 mx、my、a、b、 θ 和 f 並不能使系統有效的收斂，所以我們使用 4個步驟來調整參數 [10]。首先調整圖形的中心參數 (mx, my)，接著調整參數 a 和 b，然後調整圖形旋轉角 度參數 θ，最後調整圖形大小參數 f。 圓與橢圓圖形中的參數 a、b 和 f 同時改變正負號，並不會改變原本的形狀， 如表 3.1 所示。所以我們在演算法中固定 f ≥0，以減少搜尋的範圍。其演算法詳細過程如下: Algorithm: 使用SA 演算法偵測影像中圓與橢圓圖形之參數向量， Input: 所有在影像中的點 (N 個點)，設定K為圖形的數目。 Output: K 個橢圓(包含圓)參數向量。 Step1:初始值設定 1.9. 在 t= 1時，給定一夠高的溫度，即 T(1) = Tmax。 1.10. 使用式(8)為降溫公式。 1.11. 產生參數向量 p1, p2, ..., pk, …, pK, 其中 pk = [mk,x, m k,y, ak, bk, θk, fk]T，一 個 p 表示為一個圖形，設定 P = [p1, p2, ..., pk, …, pK]。 1.12. 使用式(5)、式(6)和、式(7c)計算 E(P)。 Step2:隨機改變參數向量並在每一個溫度中決定新的參數向量。 For m = 1 to Nt (在每個溫度中嘗試 Nt 次) For k = 1 to K (K 為影像中的圖形個數) (a) 隨機改變第 k 個圖形的中心: mn T y k x k T k,y x k m' m m m' ] =[ ] +α [ _, _, _, (15) 其中 n = [n1 n2]T是一個 2 × 1 的隨機向量，n1 和 n2是高斯隨機變數在 N(0, 1)範圍內，αm為一常數。改變後 p’k = [m’k,x, m’k,y, ak, bk, θk, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式(6)和、式(7c)計算 E(P′)。使用

Metropoliscriterion判斷是否接受 P′。若∆E = E(P′) - E(P) ≤ 0，則接受 P′。若∆E = E(P′) - E(P) > 0，計算接受機率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持 原本的 P。 (b) 隨機改變第 k 個圖形的形狀: abn T k k T k k b a b a ' '] =[ ] +α [ (16)

(31)

且使用 ' ' |a bk k |作正規化。其中 n = [n1 n2] T 是一個 2 × 1 的隨機向量， n1 和 n2是高斯隨機變數在 N(0, 1)範圍內，αab為一常數。限制 a’k > 0 和 b’k > 0 來表示圓與橢圓。改變後 p’k = [mk,x, mk,y, a’k, b’k, θk, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式(6)和、式(7c)計算 E(P′)。使用

Metropoliscriterion判斷是否接受 P′。若∆E = E(P′) - E(P) ≤ 0，則接受 P′。 若∆E = E(P′) - E(P) > 0，計算接受機率 prob = exp[-∆E/T(t)]，並產生一 平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持原本 的 P。 (c) 隨機改變第 k 個圖形的旋轉角度: θ'_k=θ_k +α_θn (17) 其中 n 是一個高斯隨機變數在 N(0, 1)範圍內，αθ為一常數。改變後 p’k = [mk,x, mk,y, ak, bk, θ’k, fk]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式 (6)和、式(7c)計算 E(P′)。使用 Metropoliscriterion判斷是否接受 P′。若

率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持原本的 P。 (d) 隨機改變第 k 個圖形的中心: f'_k =| f_k +α_fn| (18) 其中 n 是一個高斯隨機變數在 N(0, 1)範圍內，αf為一常數。改變後 p’k = [mk,x, mk,y, ak, bk, θk, f’k]T， P′ = [p1, p2, …, p′k,…, pK]。使用式(5)、式 (6)和、式(7c)計算 E(P′)。使用 Metropoliscriterion判斷是否接受 P′。若

率 prob = exp[-∆E/T(t)]，並產生一平均分布在(0,1)的隨機變數 r。若 prob > r，則接受 P′，否則維持原本的 P。

End for k End for m

Step 3: 對系統降溫

(32)

Yes 開始 設定 K 個參數向量: P = [p1, p2, …, pk, …,pK]，其中 pk = [mk,x mk,y ak bk θk fk]T。 設定初始溫度 Tmax. 定義溫度下降公式。 隨機改變第 k 個參數向量， k = 1, 2, …, K。 使用降溫公式來降溫。所有圖形之參數向量皆被隨機改變? 完成一次循環。 循環次數達到設定的總循環次數 Nt? 停止 No Yes No 隨機改變第k個圖形之 mx 和 my: pk′ = [mk,x′ mk,y′ ak bk θk fk]T, P′ = [ p1, p2, …, pk′, …,pK]。

隨機改變第k個圖形之 a 和 b: pk′ = [mk,x mk,y ak′ bk′ θk fk]T, P′ = [ p1,

p2, …, pk′, …,pK]。限制 ak′ > 0，bk′ > 0。

隨機改變第k個圖形之 θ: pk′ = [mk,x mk,y ak bk θk′ fk] T

, P′ = [ p1, p2, …,

pk′, …,pK]。

隨機改變第k個圖形之 f: pk′ = [mk,x mk,y ak bk θk fk′]T, P′ = [ p1, p2, …,

pk′, …,pK]。

計算∆E=E(P′)-E(P), 使用 Metropolis criterion 決定是否接受 P′。 開始一循環。

降溫次數達到設定的總降溫次數?

No Yes

(33)

3.4.4. 使用模擬退火演算法尋找圖形參數之範例我們將舉例說明如何使用模擬退火演算法來調整直線、圓、橢圓、和雙曲 線中的個別參數。在直線方面，需調整斜率參數 a 、 b 與參數 c。在橢圓(包含圓)和雙曲線方面，需調整中心參數(mx, my)、形狀參數 a 和 b、旋轉角度參 數 θ、與圖形大小參數 f。因為橢圓(包含圓)與雙曲線的參數皆相同，只在形狀 參數 a 和 b 的限制上不同，所以我們只用橢圓的例子來說明。 在調整直線參數時，首先調整參數 a 和 b。隨機調整參數 a 和 b，得到 一個新的形狀參數 a′ 和 b′。原本的參數向量 LP = [lp1] = [a1, b1, c1]T 則變為 LP′ = [lp1′] = [a1′, b1′, c 1]T,。如圖 3.8 所示，新的能量 E(LP′)大於原本的能量 E(LP)，

所以我們使用 Metropolis criterion 來決定是否接受 LP′。假設 prob <r，不接受 LP′，使得 LP←LP= [a1, b1, c1]T。

E(LP) E(LP′) ∆ =E E(LP′)- E(LP) > 0

prob=exp[-∆E/T(t)] < r Reject LP′; keep LP (a) (b) (c) 1 1 1 1 1 1 a a b b c c        ₌             隨機調整參數 a1和 b1。 0 E ∆ > ，prob r< ， 不接受 LP’。 1 1 1 a b c           ' 1 ' 1 1 a b c           圖 3.8. 調整參數 a 和 b。(a)原始參數，(b)新的參數，(c)不接受新的參數。 再來調整參數 c。隨機調整參數 c，得到一個新的參數 c′。原本的參數向量 LP = [lp1] = [a1, b1, c1]T 則變為 LP′ = [lp1′] = [a1, b1, c1′]T。如圖 3.9 所示，新的

(34)

E(LP) E(LP′) ∆ =E E(LP′)- E(LP) < 0 Accept LP′ (a) (b) (c) 1 1 1 1 ' 1 1 a a b b c c        ₌             隨機調整參數c₁。 ∆ <E 0，接受 LP’。 1 1 1 a b c           1 1 ' 1 a b c           圖 3.9. 調整參數 c。(a)原始參數，(b)新的參數，(c)接受新的參數。 在調整橢圓參數時，首先調整中心參數(mx, my)。隨機調整中心參數(mx, my)， 得到一個新的中心參數(mx′, my′)。原本的參數向量 P = [p1] = [m1,x, m1,y, a1, b1, θ1, f1]T 則變為 P′ = [p1′] = [m1,x′, m1,y′, a1, b1, θ1, f1]T。如圖 3.10 所示，新的能量 E(P′) 小於原本的能量 E(P)，所以我們接受 P′，使得 P← P′ = [p1′] = [m1,x′, m1,y′, a1, b1, θ1, f1]T 。 E(P) E(P′) ∆E =E(P')−E(P)≤0 Accept P’ (a) (b) (c) , , , , ' 1 x 1 x ' 1 y 1 y 1 1 1 1 1 1 1 1 m m m m a a b b f f θ θ                 =                         隨機調整參數 m1,x和 m1,y。 0 E ∆ ≤ ，接受 P’。 , , 1 x 1 y 1 1 1 1 m m a b f θ                     , , ' 1 x ' 1 y 1 1 1 1 m m a b f θ                     圖 3.10. 調整中心參數(mx, my)。(a)原始參數，(b)新的參數，(c)接受新的參數。

(35)

接著調整形狀參數 a 和 b。隨機調整形狀參數 a 和 b，得到一個新的形 狀參數 a′ 和 b′。原本的參數向量 P = [p1] = [m1,x, m1,y, a1, b1, θ1, f1]T 則變為 P′ =

[p1′] = [m1,x, m1,y, a1′, b1′, θ1, f1]T,。如圖 3.11 所示，新的能量 E(P′)大於原本的能

量 E(P)，所以我們使用 Metropolis criterion 來決定是否接受 P′。假設 prob ≥ r， 則接受 P′，使得 P← P′ = [p1′] = [m1,x, m1,y, a1′, b1′, θ1, f1]T。 E(P) E(P′) ∆E =E(P')−E(P)>0 prob=exp[-∆E/T(t)] ≥ r Accept P′ (a) (b) (c) , , , , 1 x 1 x 1 y 1 y ' 1 1 ' 1 1 1 1 1 1 m m m m a a b b f f θ θ                 =                         隨機調整參數 a1和 b1。 0 E ∆ > ， prob r≥ ， 接受 P’。 , , 1 x 1 y 1 1 1 1 m m a b f θ                     , , 1 x 1 y ' 1 ' 1 1 1 m m a b f θ                     圖 3.11. 調整形狀參數 a 和 b。(a)原始參數，(b)新的參數，(c)接受新的參數。 再來調整旋轉角度參數 θ。隨機調整旋轉角度參數 θ，得到一個新的中心 參旋轉角度參數 θ′。原本的參數向量 P = [p1] = [m1,x, m1,y, a1, b1, θ1, f1]T 則變為 P′ = [p1′] = [m1,x, m1,y, a1, b1, θ1′, f1]T。如圖 3.12 所示，新的能量 E(P′)小於原本的能 量 E(P)，所以我們接受 P′，使得 P← P′ = [p1′] = [m1,x, m1,y, a1, b1, θ1′, f1]T。

(36)

E(P) E(P′) ∆E =E(P')−E(P)<0 Accept P′ (a) (b) (c) , , , , 1 x 1 x 1 y 1 y 1 1 1 1 ' 1 1 1 1 m m m m a a b b f f θ θ                 =                         隨機調整參數 1 θ 。 ∆ <E 0，接受 P’。 , , 1 x 1 y 1 1 1 1 m m a b f θ                     , , 1 x 1 y 1 1 ' 1 1 m m a b f θ                     圖 3.12. 調整旋轉角度參數 θ。(a)原始參數，(b)新的參數，(c)接受新的參數。 最後調整圖形大小參數 f。隨機調整圖形大小參數 f，得到一個新的圖形大 小參數 f′。原本的參數向量 P = [p1] = [m1,x, m1,y, a1, b1, θ1, f1]T 則變為 P′ = [p1′] =

[m1,x, m1,y, a1, b1, θ1, f′1]T,。如圖 3.13 所示，新的能量 E(P′)大於原本的能量 E(P)，

所以我們使用 Metropolis criterion 來決定是否接受 P′。假設 prob < r，不接受 P′， 使得 P← P，仍維持原本的參數向量。

(37)

E(P) E(P′) ∆E =E(P')−E(P)>0 prob=exp[-∆E/T(t)] < r Reject P′; keep P (a) (b) (c) 隨機調整參數 1 f 。 0 E ∆ > ，prob< r 不接受 P’，保留 P。 , , , , 1 x 1 x 1 y 1 y 1 1 1 1 1 1 1 1 m m m m a a b b f f θ θ                 =                         , , 1 x 1 y 1 1 1 1 m m a b f θ                     , , 1 x 1 y 1 1 1 ' 1 m m a b f θ                     圖 3.13. 調整圖形大小參數 f。(a)原始參數，(b)新的參數，(c)不接受新的參數。 3.5. 循序式震測圖形偵測系統 3.5.1. 系統架構我們提出循序式的震測圖形偵測系統來偵測單炸點的震測圖形。因為單炸點的震測圖形中的直接波為兩條相交的直線，所以我們使用雙曲線的漸進線來偵測。 圖 3.14 為偵測震測圖形種類的順序。因為直接波為漸進線的雙曲線且反射波為 南北開口的雙曲線，所以我們在偵測圖形時，並不會調整到 θ 參數，θ 參數的 直接固定為零。 圖 3.14. 偵測震測圖形種類的順序。 輸入資料使用雙曲線的漸進線來偵測直接波剩餘資料偵測雙曲線的反射波

(38)

首先，我們依序偵測影像中的直線的直接波和雙曲線的反射波，並設定圖 形參數初始值。接著，輸入影像中的所有點座標和影像中的圖形個數 K。設定 門檻距離和門檻點數。再來使用 SA偵測圖形的參數，若要偵測的圖形是直線的直接波，則使用偵測雙曲線的漸進線參數的 SA演算法。若要偵測的圖形是雙曲線，則使用偵測雙曲線參數的 SA演算法。接著將已偵測到的圖形所包含的點去除，剩下的點再使用 SA偵測圖形的參數。重複此步驟直到達到我們設定的執行次數時，則換偵測下一個種類的圖形。直到所有圖形種類皆被偵測過後，則停止系統。

(39)

移除已偵測到的 m 個圖形所包含的點。圖形個數變為 K = K-m。 Yes No No No Yes 設定偵測影像中的圖形種類:直線、雙曲線、和橢圓(包含圓)。設定圖形參數初始值。 輸入影像中的所有點座標和影像中的圖形個數 K。 若設定偵測直接波，則使用偵測雙曲線的漸進線圖形參數的 SA。若設定偵測反射波，則使用偵測雙曲線圖形參數的 SA。點與圖形的距離是否小於門檻距離，且圖形所包含的點數大於門檻點數? 剩餘的點座標。是否達到設定的執行次數? 是否所有形狀的圖形都已偵測過? 停止 Yes 設定門檻距離和偏移值，計算門檻點數。