核估計的最佳閾值SAR圖像海岸線提取

(遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院, 遼寧阜新 123000)

0 引言

海岸線的實時監(jiān)測對國民經(jīng)濟(jì)發(fā)展來說是非常重要的，而當(dāng)今遙感技術(shù)是完成此項任務(wù)的最佳手段[1]。在眾多傳感器中，合成孔徑雷達(dá)(Synthetic Aperture Radar, SAR)相比于光學(xué)傳感器可以在不受時間和天氣的限制下連續(xù)工作，因此其更加適合作為海岸線監(jiān)測的原始數(shù)據(jù)。但SAR圖像固有的斑點(diǎn)噪聲和海面的不均勻散射特性極大地增加了海岸線提取的難度，因此準(zhǔn)確、迅速、自動的SAR圖像海岸線提取仍然是一個具有挑戰(zhàn)性的研究領(lǐng)域。

到目前為止，SAR圖像海岸線提取方法主要分為3大類：基于邊緣、基于區(qū)域和基于閾值的方法[2]。基于邊緣的方法著重研究具有灰度跳變的邊緣，故其具有高效率和高定位精度的優(yōu)點(diǎn)[3-5]。但它們最致命的缺陷是它們檢測出的海岸線是不連續(xù)的。更重要的是，有許多位于陸地區(qū)域的虛假海岸線也被檢測了出來，這無疑會對后續(xù)的真實海岸線提取和判讀造成不小的困難。

為了克服這些缺陷，基于區(qū)域的方法就應(yīng)運(yùn)而生了。它們通常以圖像全局或某些同質(zhì)區(qū)域為研究對象，故它們具有抗噪性好以及檢測邊緣連續(xù)的優(yōu)點(diǎn)。其中它們又可再細(xì)分為：基于統(tǒng)計模型的方法、水平集法和其他法[6]。其中，其他法具有耗時長、難以適用于工程實踐的缺點(diǎn)[7]，此外，其還需要預(yù)處理操作(如濾波)來進(jìn)行輔助[8-9]，但這就使得邊緣的細(xì)節(jié)信息不可避免地遭到了一定程度的破壞，進(jìn)而使提取結(jié)果的精度下降。而由于基于統(tǒng)計模型的方法中參數(shù)的數(shù)值求解較為繁瑣且耗時較長，故需要進(jìn)一步研究與優(yōu)化才能被廣泛應(yīng)用[10]。而趙泉華等[11]所用G0分布的參數(shù)估計雖然簡單，但他們僅僅根據(jù)參數(shù)的變化而確定水陸邊界點(diǎn)使得誤差很大，提取的精度較低。水平集法的參數(shù)變化對SAR圖像中較強(qiáng)斑點(diǎn)噪聲的敏感性高，因此仍有待改進(jìn)[12-14]。

由于圖像中只有海陸兩類地物，因此將圖像分類是海岸線提取的最直接手段，而閾值法則是將圖像分為兩類的最快方法，且其可以不需要使用濾波等降低精度的預(yù)處理，故可以保留精細(xì)的海岸線細(xì)節(jié)。Al Fugura等[15]將Lee濾波后的圖像的所有像素的灰度均值作為閾值，但這種取值原則并不具有普適性。而對于著名二值化算法OTSU法，很多學(xué)者都應(yīng)用它進(jìn)行了海陸分割[16-18]。但其并非是針對海陸分割而設(shè)計的，故其獲得的閾值并非最佳閾值。為了尋找最佳閾值，陳祥等[19]在OTSU分割的基礎(chǔ)上獲得了最終閾值，但因OTSU獲得的閾值本身距最佳閾值的差距就較大，且其只考慮了海水的分布而沒有考慮陸地，故分割效果仍然不理想。為了兼顧海水與陸地的分布特性，Liu和Jezek[20]假設(shè)整幅圖像服從雙峰高斯分布，并認(rèn)為雙峰中間的“谷”所對應(yīng)的灰度值就是最佳閾值。但在很多情況下，由于陸地上復(fù)雜的地物分布與不同的散射機(jī)制，使得陸地區(qū)域是不服從高斯分布的，進(jìn)而使得到的閾值也是錯誤的。為了顧及不同的散射機(jī)制與統(tǒng)計分布，Nunziata等[21-22]組成的研究小組根據(jù)極化SAR數(shù)據(jù)的物理模型確定最佳閾值，對閾值分割法作出了較大貢獻(xiàn)[23-24]。但該方法仍然受不同統(tǒng)計分布的限制，無法獲得在任意情況下的最佳閾值。

為了擺脫分布的限制，本文借助非參數(shù)核估計為基本工具[25]，在分析了海水和陸地的分布規(guī)律后利用跳變檢測的基本思想估計出了海水和陸地分布的臨界點(diǎn)[26]，即最佳閾值?？蔀殚撝捣ㄌ崛『０毒€提供一種新思路。

1 差值核估計

對于n個數(shù)據(jù)點(diǎn)組成的點(diǎn)集{(xi,yi),i=1, 2,…,n}來說，每一個xi和yi都具有一一對應(yīng)的關(guān)系，且有不確定因素所造成的微小偏差。那么可以把它們的關(guān)系描述為如下的模型：

yi=f(xi)+εi,i=1,2,…,n

(1)

式中，f表示xi和yi對應(yīng)的關(guān)系，i為數(shù)據(jù)點(diǎn)的索引，ε表示不確定因素對第i個點(diǎn)的值所造成的偏差。假設(shè)某一點(diǎn)的y值相對于其他點(diǎn)發(fā)生了明顯的跳變，那么可以利用差值核估計來進(jìn)行跳變檢測，分別估計出跳變點(diǎn)的所在位置與相應(yīng)的跳變幅度。

1.1 核函數(shù)

(2)

式中，k(>0)用以調(diào)整xi+k所獲得的權(quán)重幅度，這里取k=0.618。采用式(2)可以將點(diǎn)i的左右兩個鄰域內(nèi)的數(shù)據(jù)點(diǎn)的y值進(jìn)行加權(quán)，使越接近xi的xi+k獲得的核函數(shù)值(權(quán)重)越大，反之越小。隨著xi逐漸接近跳變點(diǎn)的所在位置，會使得兩個鄰域內(nèi)數(shù)據(jù)點(diǎn)的加權(quán)均值之差的絕對值增大的速度呈現(xiàn)越來越快的非線性增長。最終可使估計出的跳變幅度在靠近真實跳變位置時更大，遠(yuǎn)離時更小，從而估計出的跳變位置更加接近真實跳變位置。圖1為鄰域內(nèi)不同數(shù)據(jù)點(diǎn)的所在位置獲得的權(quán)重示意圖，其中黑點(diǎn)代表xi，紅點(diǎn)代表xi+k，而綠線代表權(quán)重大小，越長的權(quán)重越大，反之越小。

圖1 鄰域內(nèi)不同數(shù)據(jù)點(diǎn)的權(quán)重

1.2 跳變檢測

以點(diǎn)i為中心點(diǎn)，其差值核估計為

(3)

(4)

M(xi)=|M1(xi)-M2(xi)|

(5)

式(3)和式(4)分別表示點(diǎn)i左鄰域內(nèi)和右鄰域內(nèi)數(shù)據(jù)點(diǎn)y值的加權(quán)平均值；式(5)表示兩個鄰域內(nèi)加權(quán)均值之差的絕對值。

假設(shè)數(shù)據(jù)點(diǎn)集中只有一處跳變，那么相對于給定尺寸的鄰域，對遠(yuǎn)離跳變位置的數(shù)據(jù)點(diǎn)，其左右鄰域內(nèi)的數(shù)據(jù)點(diǎn)值一般都是平穩(wěn)變化的，因此其M值較小且在誤差范圍內(nèi)恒定不變；對靠近跳變位置的數(shù)據(jù)點(diǎn)，存在某一鄰域內(nèi)包含跳變前與跳變后的數(shù)據(jù)點(diǎn)，使得M值變大；而對于跳變點(diǎn)，會使得兩個鄰域內(nèi)分別全部包含的是跳變前的數(shù)據(jù)點(diǎn)和跳變后的數(shù)據(jù)點(diǎn)，這時M取得最大值。由此可知最大M值所對的位置(x值)即為跳變位置i*，而跳變位置對應(yīng)的最大M值即為跳變幅度s：

(6)

s=M(xi*)

(7)

以100個數(shù)據(jù)點(diǎn)為例，它們的橫坐標(biāo)x分別是它們的序數(shù)。而x與縱坐標(biāo)y的關(guān)系滿足：當(dāng)x≤50時，y=x；而當(dāng)x>50時，y=150。因此該例中的第50個點(diǎn)就是跳變點(diǎn)，而跳變的幅度為150-50=100，如圖2(a)所示。對每點(diǎn)進(jìn)行差值核估計后的結(jié)果如圖2(b)所示，這里帶寬取h=10，可以看出在跳變點(diǎn)附近時M值開始增大，而跳變點(diǎn)處M值則達(dá)到了最大(接近100)，因此較好地估計出了跳變位置與跳變幅度。

(a) 原數(shù)據(jù)點(diǎn)集

(b) 經(jīng)跳變檢測后的數(shù)據(jù)點(diǎn)集圖2 差值核估計的跳變檢測

2 海岸線提取

2.1 最佳閾值分割

在高分辨率雷達(dá)成像后的SAR強(qiáng)度圖像中，海面雜波的分布呈現(xiàn)出一種有不同程度“拖尾”特性的近似高斯分布的現(xiàn)象[19]。由于強(qiáng)散射與斑點(diǎn)噪聲，陸地上灰度值的統(tǒng)計特性則是均值和方差都較大。這就使得陸地中某些相對較暗的像素的灰度值與海面中某些相對較亮的像素的灰度值很接近或相同，表現(xiàn)在分布直方圖上則為：隸屬于陸地分布的線段與隸屬于海面分布的線段相交在了陸地與海面中互相接近的像素灰度值附近。當(dāng)陸地整體相對較亮?xí)r，該相交處通常位于近似高斯分布中的“尾”處；而當(dāng)陸地整體相對較暗時，較暗像素的數(shù)目更多，因此通常在“尾”之前就相交了。以圖3(a)中256×256像素的SAR圖像為例，該現(xiàn)象可見于圖3(b)中的直方圖，且紅色箭頭所指處即為相交處。

(a) SAR圖像

(b) SAR圖像的直方圖圖3 SAR圖像與其直方圖

而直方圖的分布特性為：絕大多數(shù)位于海面的像素的灰度值都小于該相交處所對應(yīng)的灰度值(橫坐標(biāo))，而絕大多數(shù)位于陸地的像素的灰度值都大于該相交處所對應(yīng)的灰度值。因此可以認(rèn)為該相交處是海面和陸地分布的臨界點(diǎn)，那么該臨界點(diǎn)所對應(yīng)的灰度值就可以作為海陸分割的閾值。為了自適應(yīng)地確定該閾值，需要分析該臨界點(diǎn)周圍的幾何特性。通過觀察直方圖可知，在該臨界點(diǎn)左側(cè)一定范圍的鄰域內(nèi)由線段組成的“坡”通常較陡，而右側(cè)的“坡”通常較緩。即該點(diǎn)是“坡”從陡到緩的轉(zhuǎn)折點(diǎn)，而縱觀整個直方圖，以一定范圍的鄰域來衡量每個點(diǎn)的幾何特性時，沒有其他點(diǎn)有該臨界點(diǎn)的幾何特性。所以可認(rèn)為該直方圖中只有一處坡度從大到小的跳變點(diǎn)。為了衡量坡度，可以以直方圖中的某點(diǎn)為中心點(diǎn)，計算某側(cè)鄰域內(nèi)所有點(diǎn)的縱坐標(biāo)與中心點(diǎn)縱坐標(biāo)之差的加權(quán)平均值。若該加權(quán)平均值的絕對值較大，則某側(cè)鄰域內(nèi)坡度較陡，反之坡度較緩。設(shè){(xi,Yi),i=1, 2, …, 256}為組成直方圖的點(diǎn)集，x和Y分別為橫坐標(biāo)與縱坐標(biāo)，那么利用差值核估計來檢測從大到小的坡度跳變點(diǎn)：

(8)

(9)

M(xi)=M2(xi)-M1(xi)

(10)

其中，式(8)表示直方圖中點(diǎn)i的縱坐標(biāo)與點(diǎn)i左鄰域內(nèi)所有點(diǎn)的縱坐標(biāo)之差的加權(quán)平均值；式(9)表示點(diǎn)i右鄰域內(nèi)所有點(diǎn)的縱坐標(biāo)與點(diǎn)i的縱坐標(biāo)之差的加權(quán)平均值；式(10)表示前兩個加權(quán)平均值之差。

那么相對于給定尺寸的鄰域，當(dāng)遠(yuǎn)離坡度從大到小的跳變點(diǎn)以及遠(yuǎn)離近似高斯分布的峰時，其M1與M2相差不大，因此其M較小且波動不大；當(dāng)靠近近似高斯分布的峰時，存在某一鄰域內(nèi)包含該峰兩側(cè)的“坡”的點(diǎn)，即某一鄰域內(nèi)會同時包含縱坐標(biāo)比中心點(diǎn)大的與縱坐標(biāo)比中心點(diǎn)小的點(diǎn)，使得M減小；而對于近似高斯分布的峰，會使得兩個鄰域內(nèi)分別包含的是峰左側(cè)的“上坡”的點(diǎn)與峰右側(cè)的“下坡”的點(diǎn)，即右鄰域內(nèi)所有點(diǎn)的平均縱坐標(biāo)都小于中心點(diǎn)很多，而中心點(diǎn)的縱坐標(biāo)都遠(yuǎn)大于左鄰域內(nèi)所有的點(diǎn)平均縱坐標(biāo)，故這時M取得最小值；當(dāng)靠近坡度從大到小的跳變點(diǎn)時，會使得左、右鄰域內(nèi)分別包含的是“陡坡”的點(diǎn)與“緩坡”的點(diǎn)，其右鄰域內(nèi)絕大多數(shù)點(diǎn)的平均縱坐標(biāo)通常會與中心點(diǎn)相差不大，而中心點(diǎn)的縱坐標(biāo)都遠(yuǎn)小于左鄰域內(nèi)點(diǎn)的平均縱坐標(biāo)，這時M取得最大值。由此可知最大M所對應(yīng)的灰度值即為坡度從大到小的跳變位置i*，也即最佳閾值T：

(11)

以圖3(b)中的直方圖為例，對其進(jìn)行從大到小的坡度跳變點(diǎn)檢測后的結(jié)果如圖4(a)所示。這里帶寬取h=27，且由于SAR圖像沒有絕對精確的最佳閾值，位于坡度跳變點(diǎn)附近的位置都可以認(rèn)為是理論上的最佳閾值，因此可以通過改變帶寬來調(diào)整最佳閾值。從圖4(a)可知，最大M值所對應(yīng)的橫坐標(biāo)(估計的閾值)與理論上的最佳閾值還是很接近的。用估計的閾值對SAR圖像進(jìn)行海陸分割的結(jié)果如圖4(b)所示，可以看出陸地和海面都存在少量未被正確分割的像素(即雜散點(diǎn))，這些未被正確分割的像素就對應(yīng)于前文所述的陸地中某些相對較暗的像素與海面中某些相對較亮的像素，這個現(xiàn)象同樣證明了估計的閾值與理論閾值相契合。

(a) 坡度跳變檢測后的M值

(b) 經(jīng)閾值分割后的結(jié)果圖4 最佳閾值估計與海陸分割結(jié)果

2.2 分割結(jié)果后處理

為了去除分割后產(chǎn)生的雜散點(diǎn)，需要借助形態(tài)學(xué)中尋找最大連通域的方法。首先，將連通標(biāo)準(zhǔn)設(shè)為4鄰域連通，然后尋找到值為1的最大連通域后將其保留，那么值為1的非最大連通域則被自動去除。對應(yīng)到二值圖像中：大片的白色陸地區(qū)域被保留了，而海面中白色的雜散點(diǎn)則被去除了。接下來為了去除陸地區(qū)域中的黑色雜散點(diǎn)，需要將圖像取反，然后再次保留最大連通域，然后再次將圖像取反以還原視覺效果。如果圖像中存在多片陸地區(qū)域，那么則需標(biāo)定所有的連通域并設(shè)置閾值，大于一定閾值的連通域才能被保留，這樣則可以保留多片陸地區(qū)域且去除雜散點(diǎn)。

至此，雖然已經(jīng)獲得了理想的分割結(jié)果，但理想的海岸線目前還不能在該結(jié)果中直接提取。這主要是因為經(jīng)初始閾值分割后海陸交界處有許多的“毛刺”以及不規(guī)則形狀出現(xiàn)。為此本文設(shè)計了一種用于消除“毛刺”及不規(guī)則形狀的濾波手段。設(shè)H為去除雜散點(diǎn)后的二值分割結(jié)果，則消除過程可表示為

(12)

Hf={Hf(x,y):(x,y)∈H}

(13)

式中，N(xi,yi)表示以H或Hf中值為1的像素(xi,yi)為中心的3×3鄰域，s表示N(xi,yi)中值為1的像素數(shù)目，N0(xi,yi)表示N(xi,yi)中所有的像素的值都為0。將式(12)的操作遍歷完H后得到式(13)中的Hf，再將式(12)遍歷Hf后得到新的Hf，以此重復(fù)迭代，直到Hf不再發(fā)生變化為止。

當(dāng)N(xi,yi)經(jīng)過海陸交界的規(guī)則形狀處時，會使其中超過半數(shù)的像素都是值為1的像素，故此時需保留；反之，當(dāng)N(xi,yi)經(jīng)過海陸交界的不規(guī)則形狀處時則會將少于半數(shù)的值為1的像素的值都置為0，起到填充孔洞從而將其改變?yōu)橐?guī)則形狀的作用。

圖5為消除不規(guī)則形狀示意圖，圖5(a)為圖4(b)去除雜散點(diǎn)后的二值分割結(jié)果H；圖5(b)為圖5(a)中藍(lán)色方框所在區(qū)域的局部放大，且其中的紅框表示N(xi,yi)正位于不規(guī)則形狀處，而綠框表示N(xi,yi)正位于規(guī)則形狀處；圖5(c)為圖5(b)去除了不規(guī)則形狀后的結(jié)果，且為圖5(d)中橙色方框的局部放大；圖5(d)為圖5(a)去除了不規(guī)則形狀后的結(jié)果。至此已經(jīng)得到了最終的分割結(jié)果，然后再利用邊緣檢測算子就可以得到最終的海岸線。

(a) 去除雜散點(diǎn)后的結(jié)果 (b) 去除不規(guī)則形狀

(c) 去除不規(guī)則形狀后的結(jié)果 (d) 最終分割結(jié)果圖5 去除不規(guī)則形狀

3 實驗結(jié)果與討論

為了證明提出方法的有效性，將其分別應(yīng)用到了由RADARSAT-2和Sentinel-1A衛(wèi)星獲得的4幅不同特點(diǎn)的大尺度SAR強(qiáng)度圖像中。第1個實驗圖像尺寸為2 066×2 066像素，其特點(diǎn)是陸地區(qū)域有較為均勻的紋理結(jié)構(gòu)；第2個圖像的尺寸為1 872×2 419像素，其特點(diǎn)為有許多山位于陸地上以致陸地區(qū)域忽明忽暗；第3個圖像的尺寸為1 617×2 005像素，其陸地區(qū)域中有很多建筑物與山；第4個圖像的尺寸為2 331×2 224像素，其陸地區(qū)域主要由山與城區(qū)所組成。圖6為該4幅實驗圖像，且在4個實驗中，帶寬h皆為17。

(a) 實驗圖像1 (b) 實驗圖像2

(c) 實驗圖像3 (d) 實驗圖像4圖6 4幅大尺度實驗SAR圖像

圖7(a)～(e)分別為本文方法應(yīng)用到第1個至第4個圖像后的直方圖、最佳閾值估計結(jié)果、初始閾值分割結(jié)果、最終二值分割結(jié)果與最終提取的海岸線。為了目視效果更好，將提取的海岸線(紅線)進(jìn)行了加粗并疊加到了原圖上。

為測試提出方法的準(zhǔn)確性，采用Modava和Akbarizadeh[17]提出的一種基于鄰域像素的評價標(biāo)準(zhǔn)。首先將手畫的海岸線作為用于評價的正確海岸線，然后以正確海岸線為中心，9個像素為半徑的范圍作為評價區(qū)，計算提出方法所提取的海岸線像素分別落入不同半徑評價區(qū)的累加百分比，精度評價結(jié)果如表1所示。

表1 提取的4幅SAR圖像的海岸線像素分別落入不同半徑評價區(qū)的累加百分比 %

(a) 直方圖

(b) 最佳閾值估計

(c) 初始分割結(jié)果

(d) 最終分割結(jié)果

(e) 海岸線圖7 4個實驗的相關(guān)結(jié)果

由表1可知，本文方法對4幅大尺度SAR圖像的提取精度均達(dá)到了90%以上，可以滿足工程上的應(yīng)用需求。

為了證明提出算法相對其他閾值分割法的有效性，將OTSU算法、陳祥等[19]提出的方法、Al Fugura等[15]提出的方法與Liu和Jezek[20]提出的方法作為對比方法來分割上述4幅大尺度SAR圖像。其中依次用上述幾個對比方法進(jìn)行初始分割的結(jié)果與保留最大連通域后的最終分割結(jié)果分別如圖8～圖11所示。

從對比方法的結(jié)果可以看出，OTSU算法因為并非是針對海陸分割而設(shè)計的，因此陸地上的大量暗區(qū)也被誤分為海域，導(dǎo)致后面的保留最大連通域的步驟中將大量的陸地區(qū)域都去除了，根本無法直接用來實現(xiàn)海陸分割；而陳祥等提出的方法由于是借助OTSU來確定最佳閾值的，故也無法較好地實現(xiàn)海陸分割；Al Fugura等提出的方法由于事先使用了Lee濾波操作來去噪，所以使得利用圖像的灰度均值作閾值的分割結(jié)果好于前兩種對比方法，但這無疑會降低原始圖像的精度與最終提取的精度，且仍然無法避免陸地上暗區(qū)的錯分現(xiàn)象；Liu和Jezek使用雙峰高斯分布無法較好地擬合SAR強(qiáng)度圖像中的分布，所以初始的與最終的分割結(jié)果也很不理想。因此直接利用這些閾值分割法無法獲得理想的海岸線。

相比之下，本文提出的方法無須事先對圖像進(jìn)行濾波等預(yù)處理，也并不依賴于一種或幾種特定的統(tǒng)計分布。而只是根據(jù)SAR強(qiáng)度圖像的直方圖中海域與陸地分布交界處有轉(zhuǎn)折的特點(diǎn)，尋找到該轉(zhuǎn)折處所對應(yīng)的灰度值，因此更接近理論上的最佳閾值。最佳閾值分割結(jié)果中有少量的雜散點(diǎn)，結(jié)合形態(tài)學(xué)處理與本文設(shè)計的去除不規(guī)則形狀方法就可以獲得連續(xù)、光滑的理想海岸線。

(a) 初始分割

(a) 初始分割

(b) 最終分割圖9 陳祥等提出的方法的初始分割結(jié)果與最終分割結(jié)果

(a) 初始分割

(b) 最終分割圖10 Al Fugura等提出的方法的初始分割結(jié)果與最終分割結(jié)果

(a) 初始分割

(b) 最終分割圖11 Liu和Jezek提出的方法的初始分割結(jié)果與最終分割結(jié)果

4 結(jié)束語

針對現(xiàn)有的基于閾值的SAR圖像海岸線提取方法具有的難以確定最佳閾值的問題，本文首先介紹了基于差值核估計的跳變檢測方法，然后分析了包含有海域和陸地的SAR強(qiáng)度圖像的直方圖分布特點(diǎn)。結(jié)合這種一維直方圖獨(dú)特的分布特點(diǎn)，再借助于跳變檢測方法從而估計出了SAR強(qiáng)度圖像海陸分割最佳閾值，進(jìn)而利用該閾值完成了海陸初始分割并利用設(shè)計的后處理方法實現(xiàn)了較為理想的最終分割，從而得到了海岸線。最后以RADARSAT-2和Sentinel-1A衛(wèi)星獲得的SAR強(qiáng)度圖像作為實驗數(shù)據(jù)并進(jìn)行了相關(guān)實驗，實驗結(jié)果表明，本文提出的方法可以快速有效地提取大尺度復(fù)雜SAR強(qiáng)度圖像中的海岸線，且對實驗結(jié)果進(jìn)行的精度評價也顯示提出方法的提取精度基本都達(dá)到了90%以上，從而驗證了該算法在工程上的可行性與有效性。