傅興玉 尤紅建 付 琨
①(中國科學院空間信息處理與應用系統技術重點實驗室 北京 100190)
②(中國科學院電子學研究所 北京 100190)
③(中國科學院研究生院 北京 100049)
由于SAR圖像乘性斑點噪聲的影響,傳統基于梯度的邊緣檢測算子不具有恒定虛警率的特性,難以確定邊緣的準確位置,檢測效果不理想[1]。為此,文獻[2]提出了基于ROA (Ratio Of Average)的SAR圖像邊緣算子以克服基于梯度的邊緣算子對乘性噪聲敏感的缺點。之后陸續出現了許多以ROA為基礎的SAR圖像邊緣檢測算子,比如MROA(Modified Ratio Of Averages),RGOA(Ratio and Gradient Of Averages),ROEWA(Ratio Of Exponentially Weighted Averages)[3]。文獻[4]提出了一種基于最佳似然函數的SAR圖像邊緣檢測算法,并提出基于FWSE(Fixed Window Scanning Edge)或SWCE (Scanning Window Central Edge)的窗口分析提高邊緣定位精確度。文獻[5]利用改進的ROEWA算法,在計算邊緣強度的同時,利用二次曲線進行估計邊緣方向。
基于ROA的SAR圖像邊緣算子,利用4個方向兩側均值最小比值作為邊緣強度,對SAR圖像斑點噪聲具有一定的抑制能力,但沒有充分考慮到像素周圍均值分布和邊緣走向,對噪聲的抑制能力有限。本文提出了一種基于鄰域均值連續差分平方和的邊緣算子,以提高SAR圖像邊緣檢測性能。首先,將滑動窗口分成多個互不重疊的子區域,采用鄰域均值連續差分平方和作為像素的邊緣強度,并采用分區均值差別最大的方法估計邊緣走向。然后,根據邊緣走向對邊緣強度圖像進行細化,消除真實邊緣附近的虛假邊緣,并采用基于平均強度變化率的自適應雙閾值連接算法提取邊緣。仿真和實測SAR圖像實驗結果表明,本文邊緣算子在不同亮度的同質區域表現出恒定虛警率的特性,具有較高的檢測率和邊緣定位精度,提取得到的邊緣線段的連續性保持也較好。
基于鄰域均值連續差分平方和(均方連續差分)的邊緣算子分析如下:對大小為N×N的滑動窗口,首先將其分割成多個互不重疊的等面積子區域,然后計算各子區域像素均值,并對均值連續求差,以連續差值平方和作為中心像素的邊緣強度,如式(1)所示。

其中n表示鄰域分塊的個數,μi,μj分別表示標號為i和j分塊的均值,mod表示求余操作。圖1給出了將9×9的滑動窗口劃分成8個鄰域分塊的示意圖。

圖1 鄰域分區示意圖 (分塊數量n=8)
本文邊緣算子統計像素點周圍鄰域均值的分布特性來計算邊緣強度,通過檢驗鄰域均值是否一致來判斷窗口中心像素是處在同質區域內部還是區域邊界。對于每個鄰域子區域,本文假定其滿足同質性要求,根據中心極限定理,若干同分布的獨立隨機變量之和漸進正態分布,則有

其中X表示SAR圖像像素值對應的隨機變量,E表示X的數學期望,D表示其方差,m表示樣本數量,對應子分區內的像素數目,μi表示標號i分塊的均值。

設ξ=E(X),σ2=D(X),則有

如果像素點處在同質區域內部,則鄰域分塊均值具有相同正態分布,同時鄰域分塊互不重疊,均值互相獨立。由正態分布的可加性可知

則有


因此,通過正交變換可以將式(1)中的邊緣強度轉換成n個互相獨立的隨機變量之和[6]。由 Fisher定理可知,獨立同分布的N(0,σ2)變量的線性組合仍服從零均值的正態分布,并且互相統計獨立。因此,式(1)中邊緣強度滿足自由度為n的卡方分布,其分布密度函數如式(8)所示

其中Γ表示 Gamma函數,σ表示區域方差,n表示鄰域分塊數量,NL表示滑動窗口的像素數量。
如果像素點處在區域邊界上,則鄰域分塊可以分成兩類,分別記為類別I和類別II,每個類別內部的鄰域分塊具有相同的分布。因此,式(1)中的邊緣強度由以下3部分組成,類別I內部的均方差,類別II內部的均方差,以及類別I和類別II之間的均方差。對于類別I和類別II內部的均方差,其形式與同質區域類似,下面著重討論類別之間的均方差。
設類別I和類別II分布的均值分別為ξI,ξII,方差分別為,則類別I和類別II之間的均值差服從下述分布,并且互相統計獨立。


由

綜上,在同質區域內部,本文提出的邊緣算子與區域亮度無關,表現出恒定虛警率的特性,而在邊界處強度值與均值差平方有關,具有較大響應。
邊緣走向是邊緣一個重要量度,對后續邊緣細化和邊緣連接起著重要的作用。與傳統使用梯度變化最小的方向作為邊緣走向[7,8]不同的是,本文提出一種基于類間均值差別最大的邊緣方向估計方法,如圖2所示。給定兩個方向,將鄰域子塊分成兩組,分別統計兩組均值的算術平均值,使用均值差別最大的方向組合作為像素點的邊緣方向。對于圖1的窗口劃分方式,則共有28種可能的方向組合,圖2給出部分邊緣走向示意圖。
由于滑動窗口的影響,單邊緣像素點會產生響應多次,產生邊緣展寬效應,一般借助邊緣細化算法對邊緣強度圖進行細化。ROA算子僅能估計4種不同邊緣走向,當實際走向與窗口劃分方向不一致時,會導致真實邊緣點被消除,造成邊緣線段斷裂。與傳統ROA算子不同的是,本文方法可以更好地估計出真實邊緣點走向,較好地在保留真實邊緣的同時,消除真實邊緣附近的虛假邊緣。
對于圖像上的每個像素點(x,y),記其邊緣強度值為s(x,y),邊緣方向記為(dir1,dir2),然后沿著邊緣走向中心線方向比較邊緣強度值,搜索距離為d,如圖3所示。記搜索路徑上的像素點(x',y')的強度值為s(x',y'),如果s(x,y)不小于搜索路徑上的任一點的強度s(x',y'),則認為該像素點是局部極大值點,予以保留。否則認為點(x,y)不是局部極值點,予以消除。搜索距離d是一個尺度參數,一般取窗口大小N的1/2。
如果點(x',y')正好落在像素點上,s(x',y')等于相應像素點的邊緣強度值。否則,采用雙線性插值計算點s(x',y')的邊緣強度值,如圖3(d)的情況所示。
由于SAR圖像受斑點噪聲影響嚴重,在灰度變化緩和的邊界處,真實邊緣點在細化過程中可能被錯誤地消除,導致邊緣不連續。通過實驗發現,邊緣斷續主要是由于像素點的估計方向與實際方向不一致產生的,本文采用下述方法對抑制路徑方向進行修正。

圖2 像素點不同邊緣方向示意圖

圖3 不同方向組合下邊緣細化示意圖 (箭頭表示估計的邊緣方向,虛線表示邊緣細化路徑)
對于每個像素點,在非極大值抑制過程中,如果發現抑制路徑上已經存在多個被標記的邊緣像素點,則認為當前像素點邊緣方向與實際不一致,應予以修正。將中心線方向作為像素點的修正走向,并沿著中心線的垂直方向進行非極大值抑制,修正示意圖如圖4所示。
通過邊緣細化處理,消除了真實邊緣附近的虛假邊緣,但細化結果中存在部分漏檢邊緣像素,導致邊緣線段不連續,本文采用一種自適應雙閾值連接算法連接斷裂的邊緣線段,提高邊緣的完整性。傳統Canny算子進行邊緣連接時,需要確定高,低閾值參數,不同的閾值對邊緣檢測的結果影響很大[9]。同一圖像尤其是SAR圖像中,不同邊緣取得最佳邊緣檢測效果的閾值各不相同,如果采用固定閾值則會檢測出虛假邊緣或丟失局部邊緣。本文采用基于平均強度的方法對上述方法進行改進,算法的主要步驟如下:
首先,確定高閾值系數Rhigh,其定義為邊緣強度值不小于閾值的像素點所占總像素的比例

其中f(x)表示邊緣強度的概率分布密度函數,Thigh表示高閾值,通過邊緣強度的累計分布函數和高閾值系數可求出該值。將邊緣強度不小于Thigh的像素定義為強邊緣點,并標記為種子點。
其次,設定邊緣線段平均強度最大變化率Δmax,其值為0到1之間的小數。然后從種子點出發進行邊緣連接。以連接點八鄰域內未連接像素中強度變化最小的像素作為下一個連接點。如果下一連接點為種子點或者平均強度的變化率不大于Δmax,連接該點,并更新邊緣線段的平均強度,繼續跟蹤。如果在八鄰域內找不到滿足條件的邊緣像素點,結束跟蹤。邊緣線段平均強度的變化率計算公式如式(13)所示

其中k和μk分別表示連接前邊緣線段的像素數目和平均強度,s表示下一連接像素點的邊緣強度值。
最后,按照線段長度,平均強度等條件對提取結果進行過濾,消除長度較小或者孤立的邊緣像素。
上述自適應雙閾值連接算法中使用了兩個比例參數,高閾值系數Rhigh和平均強度最大變化率Δmax。Rhigh用來控制邊緣連接的種子點數目,在進行邊緣連接之前已進行了邊緣細化處理,真實邊緣附近的虛假邊緣像素已被較好地消除,在細化后的圖像中,真實邊緣所占比例較大,因此高閾值系數一般取值較大,例如0.7~0.9即可達到較好的效果。Δmax用來動態確定邊緣連接過程中弱邊緣連接閾值,對于不同的邊緣線段或者同一線段連接不同階段,弱邊緣連接與否的閾值不盡相同,實現了根據邊緣線段的平均強度進行自適應的閾值處理。其值一般根據弱邊緣的保留程度進行設定,在實際中可設定為低閾值與高閾值的比例因子,一般取0.5即可。
優秀SAR圖像邊緣算子,一方面在真實邊緣處取得極大響應,并盡可能地消除虛警邊緣,具有較高的檢測率,另一方面檢測邊緣距離真實邊緣盡可能近,具有較高的定位精度。文獻[10]給出了邊緣算子檢測率的評估準則,

其中E為檢測結果和真實邊緣集合的交集,FN為漏檢集合,定義為未檢測到的真實邊緣像素集合,FP為誤檢集合,定義為實際為背景像素點,卻誤檢為邊緣像素點的集合。card表示集合的元素數目。Pdet是處在0到1之間的小數,如果真實邊緣都被正確檢測,同時沒有背景像素點被誤檢為邊緣點,則有Pdet=1。Pdet越接近0,說明有更多的真實邊緣點被漏檢或者背景像素點被誤檢為邊緣點。
對于每個真實邊緣點,標記檢測集合中與其距離最近的邊緣點作為該點的檢測結果,并計算兩者的位置偏差d。設定距離閾值D。集合E為真實邊緣集合中距離偏差較小(d≤D)的邊緣集合;FN為真實邊緣中實際位置與檢測位置偏差較大(d>D)的邊緣像素集合;FP為不存在真實邊緣與其對應的檢測邊緣集合。
邊緣算子得到的邊緣應該盡可能接近實際邊緣點的位置。為了評價算子的邊緣定位精度,本文對于真實邊緣像素點,計算其實際位置(xtrue,yt
rue)和檢測位置(xdet,ydet)的偏差,以正確檢測集合E中的所有邊緣點的平均偏差作為邊緣定位精度的評價準則。

其中n表示集合E中的像素數目。PLoc是大于零的實數,其值越接近零,說明算子的邊緣定位精度越好;反之定位精度越差。
選用仿真測試SAR圖像進行實驗,仿真圖像如圖5(a)所示。分別采用Canny算子,4方向ROA算子和本文算子對仿真SAR圖像進行處理,得到圖像的梯度圖或邊緣強度圖像,分別如圖5(b),5(c)和5(d)所示。ROA算子和本文算子的滑動窗口大小均為9×9。
通過對比可以發現,與Canny算子相比,ROA算子和本文算子得到的邊緣強度與區域亮度無關,在高低亮度同質紋理區域內部的響應都比較小,而在邊界處具有局部極大響應。圖6給出了 4方向ROA算子和本文算子在低亮度同質區域,高亮度同質區域以及邊界上的邊緣強度分布直方圖。不難看出,本文算子的邊緣強度要明顯優于 4方向 ROA算子,前者得到的邊緣強度分布曲線較為平滑,波動較小,表明對斑點噪聲的抑制能力較好。由于滑動窗口的影響,本文算子和 ROA算子得到的邊緣響應較寬,借助細化算法可消除真實邊緣附近的虛假邊緣,確定邊緣的準確位置。
圖8(a)和 8(b)分別給出了使用 Canny算子和ROA算子進行邊緣檢測的結果圖像。其中,Canny算子使用的參數如下:高斯平滑窗口噪聲標準差為0.55, 高閾值為380,低閾值為190。對于ROA算子,以取得最小比值的窗口劃分方向作為邊緣走向(如圖7所示),并沿著與走向垂直的方向進行細化,最后采用本文自適應雙閾值連接算法對細化結果進行連接,高閾值系數Rhigh=0.82,最大平均強度變化率Δmax=0.5,最后得到ROA邊緣檢測結果如圖8(b)所示。

圖5 不同梯度算子實驗結果對比圖

圖6 ROA算子和本文算子在不同類型區域的分布直方圖

圖7 ROA算子邊緣走向示意圖 (中間白色線框表示邊緣走向,虛線表示邊緣細化方向)

圖8 仿真SAR圖像邊緣提取實驗
圖8(c)是采用本文邊緣細化算子對得到的邊緣強度圖像進行細化的結果,抑制距離為5像素。圖8(d)進行自適應雙閾值連接處理后的結果圖像,使用參數與上述 ROA后處理中相同。對比提取結果不難看出,Canny算子得到邊緣結果中在高灰度同質區域邊界附近存在較多的誤檢情況,而 ROA算子和本文算子的檢測結果較為理想,誤檢和漏檢較少。相比之下,本文方法得到區域邊緣比 ROA算子更平滑,邊緣完整性保持較好,說明本文算子對SAR圖像斑點噪聲抑制能力較好。同時本文算子對拐點處的邊緣提取效果較好。
實驗1 選用一幅實測機場跑道SAR圖像進行邊緣提取實驗,原始SAR圖像如圖9(a)所示,提取的主要目的是得到SAR圖像中跑道的輪廓。圖9(b)是手工勾繪的邊緣結果,用來對邊緣提取結果進行性能評估。分別利用 Canny,ROA和本文算子對SAR圖像中的邊緣進行提取,結果分別如圖9(c),9(d)和 9(e)所示。Canny算子高斯平滑窗口同 4.2節中仿真圖像實驗相同,高低梯度閾值分別為 240和120。ROA算子和本文算子使用的滑動窗口大小均為9×9,后處理過程中的邊緣細化抑制長度為5像素,自適應雙閾值邊緣連接中高閾值系數Rhigh=0.80,平均強度最大變化率Δmax=0.5。從結果圖像不難看出,Canny算子得到邊緣存在較為嚴重的斷裂和毛刺現象,跑道邊界提取不完整。ROA算子和本文算子的提取效果較為理想,跑道輪廓提取相對完整。相比之下,ROA算子在跑道輪廓拐角處的提取效果不理想,部分邊界像素點漏檢,導致邊緣斷裂,而在這些地方本文算子提取效果較好,提取的邊緣具有較好的連續性。由于滑動窗口的影響,ROA算子和本文算子對SAR圖像中的細線結構檢測效果均不理想。

圖9 實測機場跑道SAR圖像邊緣提取實驗
實驗2 選用一幅實測海域高分辨率SAR圖像進行邊緣提取。原始SAR圖像如圖10(a)所示,圖像左下角為海水區域,右上角為陸地部分,包括機場跑道和建筑物等。圖10(b)是手工勾繪的邊緣結果圖像,圖10(c),10(d)和10(e)分別是Canny,ROA和本文算子的邊緣提取結果圖像,3種算子的實驗參數與實驗1相同。不難看出,Canny算子邊緣提取結果中,海岸線部分邊緣和其上側的整條邊界漏檢,同時高亮度建筑區域存在嚴重的誤檢情況。而ROA算子提取結果中,同質區域內部存在噪聲形成的虛警邊緣,說明對斑點噪聲的抑制能力不強,同時邊緣斷裂嚴重,提取結果不理想。雖然本文算子在緩變的邊緣存在漏檢情況,但總體的檢測效果要明顯優于Canny算子和ROA算子,海岸線,跑道等強邊緣以及位于海岸線上側的弱邊緣保持完整,同時在亮度較低的海水紋理區域和亮度較高的建筑物紋理區域,較好地抑制了斑點噪聲形成的虛警邊緣,誤檢率較低。
表1給出了上述3種邊緣檢測算子性能評估對比結果,其中用來評估邊緣檢測性能和定位精度的最近真實邊緣距離閾值D=5像素。從表中結果可以看出,本文算子和ROA算子對SAR圖像邊緣的檢測性能明顯優于Canny算子,說明Canny算子在SAR圖像上具有較多誤檢或者漏檢,在邊緣定位精度上兩者也略優于Canny算子。相比之下,本文算子比 ROA算子具有更好的檢測率和更高的邊緣定位精度。從得到的邊緣線段的數目和平均長度的結果可以看出,本文算子得到的數目最少,平常長度最長,說明本文算子提取得到的邊緣線段的連續性較好,而ROA算子和Canny算子邊緣斷裂較多,連續性保持較差。在運算時間上,由于本文算子和ROA算子都需要統計滑動窗口內的像素值,運算時間主要花費在邊緣強度計算和邊緣方向估計,運算時間要明顯長于Canny算子,在后處理階段,運算時間主要與提取邊緣線段數目和長度有關。本文算子提取的邊緣線段連續性保持最好,線段數目較少,運算較快。相比之下,本文算子邊緣強度估計運算時間明顯少于ROA算子。

圖10 實測海域SAR圖像邊緣提取實驗

表1 3種邊緣算子評估對比表
通過以上 3組實驗對比分析不難發現,對于SAR圖像邊緣檢測來說,傳統基于梯度的 Canny算子,在梯度較大的強邊緣處具有較好的檢測效果,但對于兩側亮度值差別較小的邊緣點漏檢率較高,同時不具有恒虛警性,在高亮度同質區域內部誤檢率較高,對斑點噪聲抑制能力差。基于ROA的邊緣檢測算子,沿著4個不同的方向將窗口分成互不相交的兩部分,計算兩部分均值比,以4個方向最小比值作為邊緣強度,該算子與區域均值無關,具有恒定虛警率的特性,但沒有充分考慮像素點鄰域像素值分布,對斑點噪聲的抑制能力不強,同時對邊緣走向估計有限,造成較多的邊緣斷裂。本文提出的基于鄰域均方連續差分的邊緣算子從像素鄰域均值分布特性出發,以鄰域均值連續差分平方和作為邊緣強度,邊緣強度衡量因子與區域亮度無關,表現出恒定虛警率的特性,同時可以較好地確定邊緣方向,并結合邊緣細化和自適應雙閾值連接處理,較好地提取得到SAR圖像中邊緣。通過仿真和實測SAR圖像實驗對比發現,本文提出的邊緣算子可以較好地提取SAR圖像中的邊緣,檢測率和邊緣定位精度較好,同時邊緣線段的連續性保持較好。下一步將研究多邊緣模型和多尺度方法來對本文方法進行改進,達到更好的SAR圖像邊緣檢測性能[11,12]。
[1]安成錦,杜琳琳,王衛華,等.基于融合邊緣檢測的SAR圖像線性特征提取算法[J].電子與信息學報,2009,31(6):1279-1282.An Cheng-jin,Du Lin-lin,Wang Wei-hua,et al..Linear feature extraction for SAR image based on fused edge detector[J].Journal of Electronics&Information Technology,2009,31(6): 1279-1282.
[2]Touzi R,Lopbs A,and Bousquet P.A statistical and geometrical edge detector for SAR images[J].IEEE Transactions on Geosciences and Remote Sensing,1988,26(6): 764-773.
[3]Fjortoft R,Marthon P,Lopes A,et al..Multiedge detection in SAR images[C].IEEE International Conference on Acoustics,Speech and Signal Processing,Munich,Germany,April 1997,Vol.4: 2761-2764.
[4]Oliver C J,Blacknell D,and White R G.Optimum edge detection in SAR[J].IEE Proceedings Radar,Sonar&Navigation,1996,143(1): 31-40.
[5]陳天澤,吳俞昊,李燕.基于DS理論的高分辨率SAR圖像復制背景直線邊緣提取方法[J].信號處理,2011,27(1): 94-101.Chen Tian-ze,Wu Yu-hao,and Li Yan.The linear edge extraction with complicated background in high resolution SAR images based on the DS evidence theory[J].Signal Processing,2011,27(1): 94-101.
[6]李慶海,陶本藻.概率統計原理和在測量中的應用(修訂本)[M].北京: 測繪出版社,1990: 185-186.Li Qing-hai and Tao Ben-zao.Principles of Probability Statistics and Its Application to Surveying (2nd Ed.)[M].Beijing: surveying and Cartography Press,1990: 185-186.
[7]Papari G and Petkov N.Edge and line oriented contour detection: State of the art[J].Image and Vision Computing,2011,29(2): 79-103.
[8]Zhou G Y,Cui Y,Chen Y L,et al..SAR image edge detection using curvelet transform and Duda operator[J].Electronics Letters,2010,46(2): 167-169.
[9]Sen D and Pal S K.Thresholding for edge detection in SAR images[C].International Conference on Signal Processing,Communications and Networking,Chennai,Jan.2008:311-316.
[10]楊朝輝,陳鷹.基于ROC融合準則的SAR邊緣檢測算法[J].光電子·激光,2010,21(7): 1053-1057.Yang Chao-hui and Chen Ying.Edge detection algorithm of SAR images based on ROC fusion rule[J].Journal of Optoelectronics·Laser,2010,21(7): 1053-1057.
[11]Dachasilaruk S.Multiscale edge detection for SAR image de-speckling[C].5th International Conference on Visual Information Engineering,Xi,an,China,July 2008: 582-587.
[12]鐘釩,周激流,郎方年,等.邊緣檢測濾波尺度自適應選擇算法[J].自動化學報,2007,33(8): 867-870.Zhong Fan,Zhou Ji-liu,Lang Fang-nian,et al..Adaptive scale filtering for edge detection[J].Acta Automatica Sinica,2007,33(8): 867-870.