李軍軍， 曹建農(nóng)， 廖娟， 程貝貝

(長安大學(xué)地球科學(xué)與資源學(xué)院，西安 710054)

0 引言

隨著遙感影像的空間分辨率提高，地物內(nèi)部細節(jié)得到了更好表征，空間結(jié)構(gòu)信息更加豐富，地物形狀以及與相鄰地物之間的關(guān)系得到更加詳盡的反映。因此，如何利用計算機更有效地處理這些高分影像成為當(dāng)前一項熱門研究課題。高空間分辨率遙感影像的邊緣提取及其相關(guān)的應(yīng)用，在地區(qū)規(guī)劃和防災(zāi)救災(zāi)等過程中占據(jù)著重要位置。傳統(tǒng)的邊緣檢測方法如Sobel算子、Robert算子和Canny算子在許多應(yīng)用中都具有較好的效果[1-2]。其中Canny算子在邊緣檢測中能得到較為清晰連續(xù)的邊界而得到許多學(xué)者的研究與改進[3-6]。但上述算子的本質(zhì)仍是一種基于局部梯度的算法，容易受到噪聲的影響，遙感影像的獲取過程中由于受到自然天氣條件等處界因素的影響通常帶有噪聲。而使用濾波方法對圖像去噪會平滑掉部分弱邊緣,使得一部分邊緣丟失或不連續(xù)。因此，這些方法存在著抑制噪聲與邊緣提取結(jié)果精確性間的矛盾[7]。

小波變換作為非平穩(wěn)信號分析的有力工具，其基函數(shù)具有局部性和衰減性，不僅能夠得到信號分解后的頻率，而且能夠知道頻率發(fā)生的位置[8]。圖像邊緣通常對應(yīng)高頻信息，使用小波變換能夠有效檢測圖像邊緣，并且具有較好的邊緣定位性質(zhì)，所以，采用小波變換提取圖像邊緣成為當(dāng)前一個研究方向。文獻[9-11]中提出以小波變換為基礎(chǔ)的不同影像邊緣提取方法，提取結(jié)果相較于傳統(tǒng)方法有一定提高。但是，標(biāo)準(zhǔn)小波變換方向數(shù)較少，小波系數(shù)僅被分解到水平、垂直和對角3個方向上，且對角子帶系數(shù)一般由于含有較多噪聲而舍棄。而高分遙感影像中地物內(nèi)部空間結(jié)構(gòu)信息豐富，地物邊緣存在于各個方向，使得基于小波變換的邊緣提取存在結(jié)果不完整問題。在此基礎(chǔ)上，研究人員發(fā)展了一些二維圖像的表示方法，如Curvelet，Ridgelet，Contourlet和Brushlet等[12-15]，相比于小波變換能更加精確地表示二維圖像，因此被廣泛應(yīng)用于圖像壓縮、分割及去噪等領(lǐng)域。文獻[16-18]等使用上述方法提取不同目標(biāo)的邊緣信息，但這些方法也存在著明顯的缺點，濾波器的結(jié)構(gòu)設(shè)計較為復(fù)雜，計算復(fù)雜度較高等。例如Contourlet變換是過完備的，即圖像變換后的大小約為原始圖像大小的4/3； Curvelet變換冗余度高，存在加窗效應(yīng)，實現(xiàn)效率差； Ridgelet變換是單尺度的，因此這些方法在圖像邊緣分析中都存在一定的不足[19]，尤其是對于數(shù)據(jù)量較大的遙感影像。

針對上述方法中的部分缺陷，本文設(shè)想提出一種理想的方法，通過采用具有多方向性基函數(shù)對圖像進行不同方向的一維小波變換，能夠從多個方向檢測圖像的奇異點，獲取影像邊緣信息。這樣既具有一維小波變換的計算簡單性，又具有二維不可分離小波變換的多方向性。相較于文獻[16-18]等方法具有較低的計算復(fù)雜度，同時也具有其多方向性的優(yōu)勢，在獲得更加完整的邊緣提取結(jié)果的同時減少時間消耗。

1 Directionlet變換

圖1(a)顯示了使用某一方向性基函數(shù)檢測與其方法正交或近似正交的奇異曲線邊緣，以獲取較大的邊緣系數(shù)。這種方法處理的好處在于，使用多個方向處理后選取較大邊緣系數(shù)能夠得到較為連續(xù)的邊緣，不會因為檢測方向數(shù)較少而產(chǎn)生一些曲線部分系數(shù)與兩側(cè)差異過大從而導(dǎo)致最后的邊緣不連續(xù)； 而且不具備線狀結(jié)構(gòu)的噪聲則仍會在后續(xù)形態(tài)學(xué)處理中被消除。本文從Directionlet變換[20]出發(fā)，根據(jù)其基本理論與方法完成圖1(b)所示的過程。

(a) 方向性基函數(shù)原始高分影像(b) 垂直基函數(shù)與經(jīng)過調(diào)制的高分影像

圖1影像方向性小波變換

Fig.1DirectionalwavelettransformofHSRimage

利用具有不同方向特性的小波基函數(shù)對高分影像進行變換，等價于對影像預(yù)先調(diào)制再利用水平或垂直方向的小波變換。

Directionlet變換是一種多尺度幾何變換工具，相比于其他常用的多尺度幾何變換工具除了具有多尺度、多方向和局部化的優(yōu)點以外還具有標(biāo)準(zhǔn)二維小波變換的計算簡單性。其基函數(shù)被稱為斜各向異性小波基S-AWT(n1，n2)，n1和n2分別表示沿不同方向的變換次數(shù)，當(dāng)n1=n2時表明變換是各向同性的，否則為各向異性變換。

定義由2個線性無關(guān)的方向向量d1=[a1，b1]，d2=[a2，b2](斜率分別為r1=b1/a1和r2=b2/a2，且r1和r2都為有理數(shù))，d1和d2線性組合構(gòu)成滿秩整數(shù)格Λ，其基本形式用MΛ表示為

(1)

任一整數(shù)格Λ均為立方整數(shù)格(cubic integer lattice)Ζ2的子格，即Λ?Ζ2。

由格理論[21]可知，完整的格空間被MΛ劃分為|det(MΛ)|個關(guān)于Λ的陪集(co-set，CS)。每一個陪集CS[Sk,1 ，Sk,2]可由格Λ通過偏移矢量Sk=[s1，s2](k=1，2，...，|det(MΛ)| ,s1，s2∈Ζ)得到。位于格Λ及其CS上的像素點的坐標(biāo)可以表示為(i，j)，即

i=c1a1+c2a2+Sk,1，

(2)

j=c1b1+c2b2+Sk,2，

(3)

式中：c1和c2為任意整數(shù)，c1，c2∈Ζ,具體取值視影像大小而定；Sk,1=s1；Sk,2=s2。當(dāng)固定任意c2，遍歷c1時可以得到嚴(yán)格位于d1方向上的一條陪線(co-line，CL)。同樣，由偏移矢量可以得到其他CS上的CL。反之可以得到d2方向上CL。分別沿d1和d2方向進行n1和n2次一維小波變換則完成了單尺度下的Directionlet變換，詳細內(nèi)容請參考文獻[20]。在此不再詳細討論Directionlet變換的整個過程，而是根據(jù)其基本理論方法獲取多方向一維CL用于進行多方向小波變換。

以d1=[-1，2]和d2=[1，1]所構(gòu)成的一個基本單元格為例。圖2(a)中黑色點表示了由MΛ生成的一個基本單元格，|det(MΛ)|=3表示平面可以被MΛ分解為3個子集，偏移矢量分別為S1=[0，0]，S2=[0，1]，S3=[0，2]。變換c1和c2的值完成整個分解過程得到CS[0，0]，其余2個CS[0，1]和CS[0，2]上的點可以由CS[0，0]上的點經(jīng)過偏移矢量得到。圖2(b)表示3條CL，CL上的點都嚴(yán)格位于d2方向上，同樣其余的CLSk也可以通過偏移矢量獲得。將上述方法應(yīng)用于圖像，即圖2中圓圈視為像素點，完成整個分解過程獲取任意有理斜率方向的像素點集合。

(a) 一個基本單元格(b) 陪線

圖2格與陪線

Fig.2Latticeandco-line

2 多方向小波變換

小波變換是對傅里葉變換的發(fā)展，是一種在空間域和頻率域同時擁有分辨性的多尺度分析方法。小波變換使用小波函數(shù)族及相應(yīng)的尺度函數(shù)將原信號分解成不同頻率子帶。對于連續(xù)函數(shù)f(x)，定義其小波變換為

(4)

根據(jù)Directionlet變換的基本理論，對于給定的圖像可以獲取任意有理斜率方向上的CL集合 。例如分別采用[1，0]，[2，1]，[1，1]，[1，2]4個方向(分別對應(yīng)0°，30°，45°，60°)對圖像進行多方向性小波變換，并給出變換后的子帶系數(shù)，其結(jié)果如圖3所示。

(a) 變換方向(b) 高分影像(c) [1，0]方向

(d) [2，1]方向(e) [1，1]方向(f) [1，2]方向

圖3高分影像不同方向小波變換結(jié)果

Fig.3WavelettransformresultsofHSRimageindifferentdirections

從圖3結(jié)果可以看出，在不同方向小波的作用下其變換結(jié)果不同。由于變換方向與紅色圈中道路的方向一致，在Directionlet變換方向消失矩(directional vanishing moments ，DVM)的作用下，邊緣在沿30°方向的子帶中基本消失，即小波系數(shù)較小，在與30°方向相鄰的水平和45°方向上的子帶中重現(xiàn)但數(shù)值仍不大，而在60°方向子帶上系數(shù)的響應(yīng)相對較大，且呈高亮狀態(tài)。

2.1 重定義小波模值計算方法

圖像經(jīng)變換后最重要的是需要選合適的系數(shù)進行后處理。以標(biāo)準(zhǔn)小波變換為例，fd1和fd2分別是圖像f經(jīng)過水平與垂直方向小波變換后的系數(shù)，對于圖像中的任意一點(x,y)，fd1(x,y)和fd2(x,y)正好構(gòu)成了一組梯度矢量，由此可得f(x,y)小波變換后的梯度模值Mf(x，y)和幅角Af(x，y)分別為

(5)

(6)

由于圖像中突變點主要位于圖像邊緣，這些突變點對應(yīng)于模極大值，梯度方向指向梯度模極大值，這就意味著可以通過梯度方向?qū)ふ夷O大值來檢測圖像邊緣[22]得到初始結(jié)果。

傳統(tǒng)的方法中，梯度分量只有2個，對其求模和幅角就可以確定梯度大小和方向。而在本方法中變換方向一般為8個方向或4個方向。若以4個方向為例，則圖像變換后得到4個梯度分量分別為fd1,fd2,fd3,fd4。文獻[23]對此處理方法為

(7)

但這種做法太過籠統(tǒng)，且當(dāng)變換方向較多時計算量會很大。在實際情況中，多數(shù)邊緣提取實驗中發(fā)現(xiàn)，只有少數(shù)幾個梯度分量的模值較大。這符合人們的認(rèn)知，由于邊緣一般出現(xiàn)在像素灰度值差異較大的地方，通常情況下，f(x,y)點處的梯度分量最大值出現(xiàn)在該點法線方向，而其他方向梯度分量的模沿法線兩側(cè)遞減，最小為0，如圖4所示。

圖4 多方向梯度分量Fig.4 Multi-direction gradient components

(8)

2.2 非極大值抑制

非極大值抑制是進行邊緣檢測的重要步驟，用于尋找局部梯度最大值，細化邊緣。其基本思想是： 如果圖像(x,y)位置上的系數(shù)模值大于沿梯度方向兩端的值，則認(rèn)為該像素位置為候選圖像邊緣，否則為非邊緣點。而通常梯度方向介于2個梯度分量方向之間的亞像素點位置上(如圖5(a)所示)，因此需要對梯度方向兩邊的梯度值進行插值，得到沿梯度方向上的梯度值。在本文中，選取4個變換方向{[1，0]，[1，1]，[0，1]，[-1，1]}，對圖像進行實驗并加以說明。

(a) 2個最大梯度分量方向相鄰(b) 2個最大梯度分量方向不相鄰

圖52個最大梯度分量方向的位置關(guān)系

Fig.5Positionrelationofdirectionofthetwomaximumgradientcomponents

多數(shù)情況下，梯度值最大的方向與次大梯度方向是相鄰的，如圖5(a)所示。但是當(dāng)采用多方向小波變換時，由于變換數(shù)量較多及影像本身的幾何結(jié)構(gòu)等因素，存在部分2個最大梯度方向不相鄰的情況，如圖5(b)所示。對于第一種情形，采用一般方法，對梯度分量f1和f2插值得到梯度方向上的2個梯度值M1和M2，即

(9)

(10)

比較Mf(x，y)與M1，M2值的大小確定當(dāng)前位置是否為候選邊緣點。對于第二種情況，由于系數(shù)最大方向與系數(shù)第二大方向的位置不相鄰，相關(guān)性弱，無法直接對其插值求梯度方向模值。因此，將f(x，y)處的值分別與f1和f2方向前后值比較，若為局部最大值則保留f(x，y)為候選邊緣。圖6顯示了一幅高分影像經(jīng)非極大值抑制后的結(jié)果。最后采用雙閾值法[24-25]對非極大值抑制后的邊緣圖像進行弱邊緣去除，得到主要邊緣。由于圖像中通常含有噪聲，雖然噪聲形成的邊緣像素個數(shù)較少，但會在部分邊緣產(chǎn)生毛刺現(xiàn)象，因此需要采用數(shù)學(xué)形態(tài)學(xué)的方法對其進行一定的后處理[26]，得到最終的圖像邊緣。

(a) 高分影像(b) 非極大值抑制結(jié)果

圖6高分影像與非極大值抑制后的結(jié)果

Fig.6HSRimageandresultofnon-maximumsuppression

2.3 實驗流程

本文方法具體步驟如下： ①對每個給定的方向向量d=[a，b]，將輸入圖像進行基于格的分解，獲取位于給定方向上的一維線集合； ②對第①步中每一條線進行小波變換，使用Haar小波高通濾波器，將變換后的系數(shù)恢復(fù)為圖像矩陣形式，便于計算，從而完成方向性小波變換，得到不同方向小波變換的高頻子帶； ③重新定義小波系數(shù)模值的方法計算每個像素位置模和梯度方向，用非極大值抑制確定候選邊緣點，分別采用不同方法處理2個梯度方向是否相鄰時的情況； ④采用雙閾值法得到邊緣結(jié)果； ⑤使用數(shù)學(xué)形態(tài)學(xué)方法對細節(jié)進行適當(dāng)后處理。多方向小波變換的高分影像邊緣提取具體流程如圖7所示。

圖7 多方向小波變換的高分影像邊緣提取流程Fig.7 HSR image edge extraction flow chart based on multi-direction wavelet transform

3 實驗結(jié)果與分析

本實驗在Intel Core i5 2.7 GHz，4 GB內(nèi)存的計算機上進行，編程環(huán)境為Matlab2014a。數(shù)據(jù)選擇3個地區(qū)空間分辨率為0.5 m的彩色航空影像，由于光譜特征對圖像邊緣提取影響不大，因此將影像預(yù)先轉(zhuǎn)為灰度影像，選擇地物主體類型為建筑物和道路。大小均為256像素×256像素。選擇{[1，0]，[1，1]，[0，1]，[-1，1]}4個方向進行方向性小波變換。將與基于Canny算子的邊緣提取結(jié)果和基于標(biāo)準(zhǔn)小波變換的邊緣提取結(jié)果進行對比分析。

3幅實驗數(shù)據(jù)的Canny算子邊緣提取結(jié)果、小波變換邊緣提取結(jié)果及本文方法邊緣提取結(jié)果如圖8所示。3幅實驗數(shù)據(jù)Canny算子邊緣提取的閾值參數(shù)分別選擇0.2，0.23和0.23。

(a) 實驗1高分影像(b) 實驗1Canny算子結(jié)果(c) 實驗1小波變換結(jié)果(d) 實驗1本文方法結(jié)果

(e) 實驗2高分影像(f) 實驗2Canny算子結(jié)果(g) 實驗2小波變換結(jié)果(h) 實驗2本文方法結(jié)果

(i) 實驗3高分影像(j) 實驗3Canny算子結(jié)果(k) 實驗3小波變換結(jié)果(l) 實驗3本文方法結(jié)果

圖8實驗結(jié)果

Fig.8Resultsofexperiment

由圖8(a)—(h)可以看出，由于這2類影像的幾何結(jié)構(gòu)較為簡單，主要為規(guī)則線狀結(jié)構(gòu)，邊緣與內(nèi)容部分灰度值差異明顯，3種方法均能識別主體邊緣，但前2種方法結(jié)果在某些較弱邊緣線處都發(fā)生斷裂現(xiàn)象，而本文方法則有效地避免了這種情況。而對于結(jié)構(gòu)較為復(fù)雜的高分圖像，如圖8(i)，本文方法明顯要優(yōu)于Canny算子，從圖8(l)中可以看出，本文方法獲得的結(jié)果包含更多更完整的線狀邊緣。雖然傳統(tǒng)方法中選擇較小的閾值時可以獲得更多更完整的邊緣，但是由草地和樹冠等形成的紋理區(qū)域也包含大量的高頻信息，容易形成偽邊緣，因此結(jié)果很難同時兼顧邊緣完整性與去除偽邊緣。本文方法在控制由于光照和紋理區(qū)的噪聲和偽邊緣的同時，能夠獲得更加完整的邊緣，這是因為本文方法更好地利用了圖像本身的幾何結(jié)構(gòu)，有針對性地檢測圖像的奇異信息，而噪聲和光照不均勻造成的偽邊緣通常不具備邊緣的幾何結(jié)構(gòu)。

對比采用小波變換的邊緣提取結(jié)果可以明顯看出，本文方法優(yōu)于標(biāo)準(zhǔn)小波方法，因為標(biāo)準(zhǔn)小波變換只能檢測水平和垂直方向奇異點，對于非垂直與水平的弱邊緣小波系數(shù)被分解到垂直與水平子帶，沒有充分利用圖像的幾何結(jié)構(gòu)信息。而基于多方向性小波變換可以檢測任意方向的奇異點，從中選取較大系數(shù)，有利于獲取豐富連續(xù)的圖像邊緣。

但是在多次的實驗結(jié)果對比分析中，本文方法也存在一些缺陷。例如對于某些較弱的邊緣或者較小地物的邊緣提取結(jié)果中發(fā)現(xiàn)仍有部分像素點斷裂的情況?？赡艽嬖诘脑蛟谟诓糠珠撝翟O(shè)置不恰當(dāng)，或未考慮多尺度的方法。某種程度上地物的邊緣也可以視為是多尺度的統(tǒng)一，使用單一或區(qū)間閾值時能夠得到某個尺度下較好的邊緣結(jié)果，但在多尺度下提取結(jié)果會不可避免地出現(xiàn)問題。

4 結(jié)論

本文提出一種圖像邊緣提取方法，對高空間分辨率遙感影像進行邊緣信息提取。首先分析了傳統(tǒng)邊緣提取方法和小波邊緣提取方法存在的不足之處，再根據(jù)高分遙感影像的自身幾何結(jié)構(gòu)特性提出多方向性小波變換邊緣提取方法，兼顧了一維小波變換計算簡單性和二維多尺度幾何分析多方向性。對比實驗結(jié)果表明，本文方法相對于傳統(tǒng)方法能夠更好地提取高分影像邊緣，具有應(yīng)用有效性。

小波變換是一種多尺度的表達工具，高分影像空間分辨率較高，地物內(nèi)部細節(jié)更加詳細，使用單一尺度的邊緣提取過于籠統(tǒng)。在應(yīng)用中，有時需要忽略較小地物目標(biāo)而更多關(guān)注整體的輪廓信息，有時則不然。因此，充分利用高分影像的優(yōu)勢，在不同尺度空間下提取圖像邊緣信息的更有效方法是今后仍需繼續(xù)研究的問題。