王燕紅，程博，尤淑撐，武盟盟

（1.中國(guó)科學(xué)院遙感與數(shù)字地球研究所，北京 100094；2.中國(guó)科學(xué)院大學(xué)，北京 100049；3.中國(guó)土地勘測(cè)規(guī)劃院，北京 100035）

1 引 言

近些年，城市擴(kuò)張和城市建設(shè)的迅猛發(fā)展，運(yùn)用遙感技術(shù)進(jìn)行城市環(huán)境研究具有重要應(yīng)用價(jià)值。星載合成孔徑雷達(dá)（Synthetic Aperture Radar，SAR）具有全天時(shí)、全天候、穿云透霧的優(yōu)勢(shì)，對(duì)我國(guó)幅員遼闊、氣象氣候和下墊面條件復(fù)雜、且急需快速獲取監(jiān)測(cè)數(shù)據(jù)的前提下，具有明顯的優(yōu)勢(shì)。2007年以來(lái)，加拿大 RADARSAT－2、德國(guó) TerraSAR－X及意大利COSMO－SkyMed等分辨率達(dá)1m雷達(dá)衛(wèi)星的陸續(xù)成功發(fā)射，高分辨率城區(qū)SAR圖像的獲取能力大大提升，一方面為準(zhǔn)確提取建筑區(qū)提供了前提，但是另一方面卻導(dǎo)致高分辨率條件下城區(qū)環(huán)境更為復(fù)雜，給高分辨率SAR影像應(yīng)用帶來(lái)了新的挑戰(zhàn)。

建筑區(qū)提取是城區(qū)遙感圖像解譯的重要研究課題之一。建筑區(qū)是城市環(huán)境最基本的一部分，一般由住宅區(qū)、公共建筑、道路與綠地和場(chǎng)地組成。在高分辨率SAR圖像上，灰度分布比較復(fù)雜，用傳統(tǒng)的基于圖像灰度的提取方法很容易受到居民地內(nèi)部要素以及噪聲的影響而得不到滿意的效果，針對(duì)高分辨率SAR影像中含有豐富的紋理信息特點(diǎn)，紋理分析方法是對(duì)像素間空間分布關(guān)系進(jìn)行研究，因而可利用紋理分析的方法來(lái)解決。國(guó)內(nèi)外學(xué)者也做了大量的研究，如Dell’Acqua等［1］利用灰度共生紋理度量來(lái)刻畫城區(qū)SAR圖像上不同區(qū)域的建筑覆蓋密度，從而區(qū)分城市中心、居住區(qū)和郊區(qū)。Tison等［2］利用馬爾科夫隨機(jī)場(chǎng)紋理分類方法對(duì)城區(qū)進(jìn)行精細(xì)的分類。強(qiáng)永剛等［3］利用小波變換和數(shù)學(xué)形態(tài)學(xué)提取人工建筑區(qū)。趙凌君等［4］利用變差函數(shù)的方法，實(shí)現(xiàn)了建筑區(qū)的快速提取。變差函數(shù)方法是一種較新的紋理提取方法，在建筑區(qū)和非建筑區(qū)上的特征表現(xiàn)差異大，算法簡(jiǎn)單、快速，具有重要應(yīng)用價(jià)值，然而目前常用的變差函數(shù)模型易受噪聲干擾，穩(wěn)健性差［5］。本文提出一種基于中值濾波的變差函數(shù)紋理分析方法提取建筑區(qū)，通過(guò)改進(jìn)變差函數(shù)的計(jì)算方法，使窗口的大小不再受步長(zhǎng)的約束，不僅保留了變差函數(shù)區(qū)分建筑區(qū)與非建筑區(qū)的優(yōu)勢(shì)，而且改善了高分辨率SAR提取建筑區(qū)時(shí)受強(qiáng)反射點(diǎn)、噪聲干擾大的缺點(diǎn)，提取效果更完整。

2 變差函數(shù)紋理分析

2.1 變差函數(shù)理論

變差函數(shù)理論由數(shù)學(xué)專家G.Maberon教授1962年創(chuàng)立，作為地質(zhì)統(tǒng)計(jì)學(xué)的重要工具，變差函數(shù)被應(yīng)用于空間隨機(jī)場(chǎng)的統(tǒng)計(jì)特性研究。隨著變差函數(shù)理論的不斷發(fā)展，其應(yīng)用領(lǐng)域從早期的空間數(shù)據(jù)分析、地理數(shù)據(jù)位置分析等逐漸擴(kuò)展到遙感數(shù)據(jù)分析，尤其是在20世紀(jì)90年代，Miranda等人成功將變差函數(shù)紋理分析方法應(yīng)用于SAR圖像植被識(shí)別之后［6-7］，變差函數(shù)成為了一種重要的紋理特征分析方法，被廣泛應(yīng)用于遙感圖像信息提取領(lǐng)域。

變差函數(shù)又稱半變差函數(shù)（Semivariogram Function），定義為區(qū)域化變量Z（x）和Z（x＋h）（同時(shí)包含兩點(diǎn)距離和方向信息）兩點(diǎn)之差的方差之半，即

對(duì)于離散的柵格數(shù)據(jù)，實(shí)驗(yàn)中的變差函數(shù)定義為：

其中，N（h）表示觀測(cè)數(shù)據(jù)中間距為h的點(diǎn)對(duì)數(shù)目，估計(jì)值γ＊（h）通常稱為實(shí)驗(yàn)變差函數(shù)。變差函數(shù)用于度量區(qū)域化變量的空間相關(guān)性，能充分反映圖像數(shù)據(jù)的隨機(jī)性和結(jié)構(gòu)性。圖1是一個(gè)理想的變差函數(shù)模型，當(dāng)h＝0時(shí)的變差函數(shù)值C0為塊金值（Nugget Value），代表不確定的變化，隨著h的增大，空間樣本之間的相關(guān)性越來(lái)越?。ǚ窍嚓P(guān)性越來(lái)越大），當(dāng)樣本變?yōu)橄嗷オ?dú)立時(shí)的距離a稱為變程，相應(yīng)的變差函數(shù)值逐漸增大并在間距為a時(shí)趨于平穩(wěn)值C0＋C1，這一平穩(wěn)值稱為基臺(tái)值（Sill）。

圖1 理想的變差函數(shù)模型

2.2 改進(jìn)的變差函數(shù)計(jì)算方法

變差函數(shù)用于紋理分析的常用計(jì)算方式是首先確定步長(zhǎng)h、窗口w、計(jì)算方向，然后計(jì)算窗口w內(nèi)所有間距為h的點(diǎn)對(duì)的半方差值，然后取平均（如公式（2））作為窗口中心點(diǎn)的變差函數(shù)值，遍歷全圖即得到影像的變差函數(shù)特征圖。對(duì)于高分辨率的SAR圖像而言，窗口w的選取取決于h，要保證窗口w內(nèi)間距h的點(diǎn)對(duì)數(shù)目足夠，窗口w至少應(yīng)為3h～5h［4］。然而，窗口w過(guò)大會(huì)造成圖像整體模糊、邊緣虛警率高，且窗口內(nèi)取平均的方法易受噪聲、孤立強(qiáng)反射點(diǎn)干擾，算法穩(wěn)健性差。

根據(jù)變差函數(shù)計(jì)算方法所存在的弊端，本文提出一種穩(wěn)健的變差函數(shù)計(jì)算方法，該算法繼承了變差函數(shù)和標(biāo)準(zhǔn)中值濾波方法的優(yōu)點(diǎn)，窗口w取值不再受h約束，主要思想是：在計(jì)算像素點(diǎn)的變差函數(shù)值時(shí)，不再取窗口w內(nèi)間距為h的點(diǎn)對(duì)來(lái)求取，而是計(jì)算像素點(diǎn)（x，y）與固定方向（如0°方向）上間距h的點(diǎn)（x＋h，y）的半方差值，將此值作為點(diǎn)（x，y）的半方差值；然后取窗口內(nèi)所有像素的半方差值的中值代替均值作為變差函數(shù)值。

計(jì)算公式如下：

其中，w為二維模板。

下面以窗口w＝3，h＝1，全方向?yàn)槔鐖D2所示。

圖2 變差函數(shù)紋理特征計(jì)算示意圖

3 基于改進(jìn)的變差函數(shù)紋理分析的建筑區(qū)提取方法

3.1 高分辨率SAR影像的建筑區(qū)紋理特征

遙感影像上可以分為建筑區(qū)與非建筑區(qū)兩大類。建筑區(qū)可分為城鎮(zhèn)建筑區(qū)和農(nóng)村建筑區(qū)。在高分辨率SAR中，建筑物由于角反射器及屋頂強(qiáng)反射而產(chǎn)生的亮斑，通?；叶龋炼龋┮哂谥苓叺貐^(qū)，然而城鎮(zhèn)建筑區(qū)和農(nóng)村建筑區(qū)由于建筑結(jié)構(gòu)、材質(zhì)、分布的差異，紋理特征也不盡相同。城鎮(zhèn)建筑物分布整齊，樓房之間間隔大，大多數(shù)為平頂整齊的高層樓房，使用材料大多具有良好反射率，當(dāng)屋頂在垂直于雷達(dá)觀測(cè)角的情況下產(chǎn)生的鏡面發(fā)射是強(qiáng)反射，因此多表現(xiàn)為亮斑。而高層樓房之間道路、草坪等散射回波較弱，在影像上表現(xiàn)為亮斑之間的暗斑塊，加之建筑物通常排列較為整齊，形成有一定周期規(guī)律的明暗相間的紋理特征。農(nóng)村建筑區(qū)房屋間隔小，分布相對(duì)密集，通常為房屋聚集的村莊，在影像上表現(xiàn)為棋盤式分布的亮斑，紋理特征不如城鎮(zhèn)建筑區(qū)有規(guī)律。非建筑區(qū)主要包括植被、水體等，整體亮度較低。

3.2 建筑區(qū)的變差函數(shù)曲線

高分辨率SAR圖像中建筑區(qū)、植被區(qū)和水域的實(shí)驗(yàn)變差函數(shù)曲線如圖3所示（A為建筑區(qū)，B為植被區(qū)，C為水域）。從圖中可以看出，首先，建筑區(qū)的變差函數(shù)曲線整體上都遠(yuǎn)遠(yuǎn)高于植被區(qū)和水域，且隨著間距h的增加周期性地出現(xiàn)峰值和谷點(diǎn)，但是在數(shù)值上，峰值與峰值、谷點(diǎn)與谷點(diǎn)不一定嚴(yán)格相等，隨著h的增大會(huì)表現(xiàn)出衰減的現(xiàn)象［8-9］，這是由于建筑區(qū)是富含大量強(qiáng)散射點(diǎn)，且建筑區(qū)與街道等間隔分布；植被區(qū)的變差函數(shù)曲線隨著間距h的增大，變差函數(shù)曲線先緩慢上升，隨后趨于平穩(wěn)。水體由于鏡面反射回波很弱，像素灰度值很低且非常均勻，變差函數(shù)曲線近似為一條直線。建筑區(qū)的變差函數(shù)曲線具有特殊性，與非建筑區(qū)特征有明顯差異，因此可以作為有效區(qū)分建筑區(qū)和非建筑區(qū)的依據(jù)。

圖3 實(shí)驗(yàn)區(qū)各地物的變差函數(shù)曲線

3.3 變差函數(shù)參數(shù)選取

（1）步長(zhǎng)h

如圖3（b）所示，變差函數(shù)曲線在h1＝6的位置達(dá)到第一個(gè)峰值，h2＝10的位置達(dá)到波谷。在0～h1區(qū)間內(nèi)，曲線上升非?？欤痪哂蟹€(wěn)定性，因此步長(zhǎng)h應(yīng)大于h1；另外考慮到隨著步長(zhǎng)增大，曲線逐漸衰減，相關(guān)性減小，因此本文選取稍大于h1，小于h2的位置作為步長(zhǎng)。這樣既具有穩(wěn)健性，又能保證相關(guān)性比較強(qiáng)。

（2）窗口w

一般來(lái)說(shuō)，窗口選擇太大可能包含多余的無(wú)用信息；窗口太小又不能有效、全面地描述地物的特征，而且當(dāng)窗口選擇的較小時(shí)，噪聲對(duì)紋理的計(jì)算影響較大。所以，合適的窗口大小對(duì)于紋理的計(jì)算十分重要。在實(shí)際應(yīng)用中，一般根據(jù)影像地物的實(shí)際特征來(lái)確定窗口的大小。

（3）方向

要考慮的因素是變差函數(shù)的各向同性和各向異性，即變差函數(shù)計(jì)算方向的選擇。一般采用全方向。

3.4 提取流程

建筑區(qū)提取流程如圖4所示，首先進(jìn)行變差函數(shù)實(shí)驗(yàn)，確定變差函數(shù)參數(shù)h、w，然后計(jì)算目標(biāo)影像的變差函數(shù)紋理圖，利用EM算法對(duì)得到的紋理圖進(jìn)行閾值分割，得到初步的提取結(jié)果，利用形態(tài)學(xué)方法進(jìn)行后處理，包括填洞，去除小面積區(qū)域等，最后得到建筑區(qū)提取結(jié)果。

圖4 算法提取流程

4 實(shí) 驗(yàn)

4.1 實(shí)驗(yàn)數(shù)據(jù)

為驗(yàn)證算法的有效性，選取了高分辨率圖像中農(nóng)村建筑區(qū)（實(shí)驗(yàn)區(qū)1）和城鎮(zhèn)建筑區(qū)（實(shí)驗(yàn)區(qū)2）進(jìn)行了實(shí)驗(yàn)。實(shí)驗(yàn)區(qū)1和實(shí)驗(yàn)區(qū)2均采用2007年8月19日獲取北京市升軌右視的TerraSAR－X數(shù)據(jù)，極化方式為HH，空間分辨率為1.25m，大小均為1600×1200像素。

4.2 變差函數(shù)參數(shù)的確定

根據(jù)圖3的變差曲線特征（本文所有實(shí)驗(yàn)影像均來(lái)自同一景影像），本文設(shè)定步長(zhǎng)h＝8。一般來(lái)說(shuō)，窗口w的選定根據(jù)實(shí)驗(yàn)區(qū)特征選取，目前多數(shù)依靠經(jīng)驗(yàn)選取。為了定量的分析窗口大小對(duì)提取結(jié)果的影響，本文通過(guò)固定步長(zhǎng)h，改變窗口w的大小，計(jì)算實(shí)驗(yàn)區(qū)的檢測(cè)率、虛警率、漏警率，以此選擇最優(yōu)的窗口w。

圖5為固定h＝8，檢測(cè)率、虛警率、漏警率隨窗口變化的曲線。兩實(shí)驗(yàn)區(qū)的曲線變化大致相同，檢測(cè)率曲線先上升后下降，漏警率呈下降趨勢(shì)，虛警率則是先下降后上升。當(dāng)窗口在0～20時(shí)，檢測(cè)率一直上升，大約達(dá)到23～28的時(shí)候，檢測(cè)率最高，此時(shí)虛警率和漏警率較小。因此本文選取w＝23作為最優(yōu)窗口。

4.3 實(shí)驗(yàn)結(jié)果及評(píng)價(jià)

圖5 檢測(cè)率、虛警率、漏警率隨窗口變化曲線圖

根據(jù)上節(jié)提出的方法對(duì)原始影像進(jìn)行建筑區(qū)提?。╤＝8，窗口w＝23，全方向），本文將結(jié)合光學(xué)影像人工提取的建筑區(qū)作為標(biāo)準(zhǔn)結(jié)果。實(shí)驗(yàn)區(qū)的初始提取結(jié)果如圖6（b）、圖7（b）所示，然后將初始結(jié)果進(jìn)行了填洞、去除小區(qū)域（小于300像素），與標(biāo)準(zhǔn)影像進(jìn)行疊加如圖6（c）、圖7（c）所示，并且與未改進(jìn)變差函數(shù)的方法的后處理結(jié)果進(jìn)行對(duì)比（白色部分為正確區(qū)域，藍(lán)色部分為漏警區(qū)，黃色部分為虛警區(qū)）。

從對(duì)比圖可以明顯看出，本文提出的方法更能準(zhǔn)確、完整地提取出建筑區(qū)的邊界，邊緣虛警區(qū)小，建筑區(qū)內(nèi)部的空洞很少，一定程度改善了受孤立強(qiáng)反射點(diǎn)、噪聲干擾影響的缺點(diǎn)，提取效果優(yōu)于未改進(jìn)的方法。為了使得對(duì)比結(jié)果清晰，本文在實(shí)驗(yàn)區(qū)1截取建筑區(qū)邊緣部分進(jìn)行了放大，如圖8所示，首先觀察建筑區(qū)邊界，本文的方法邊界的虛警很小，而未改進(jìn)方法明顯地?cái)U(kuò)大的邊界范圍，其次對(duì)于影像中的道路旁邊的強(qiáng)散射點(diǎn)，未改進(jìn)方法錯(cuò)誤地提取為建筑區(qū)，造成了虛警，而本文方法一定程度改善了這一缺陷。

為了定量地評(píng)價(jià)提取結(jié)果，本文采用檢測(cè)率（Detection Rate），虛警率（False Alarm Rate）［10］和漏警率（Missing Alarm Rate）3個(gè)度量值進(jìn)行結(jié)果評(píng)價(jià)，評(píng)價(jià)結(jié)果如表1所示。從結(jié)果可以看出，改進(jìn)的方法檢測(cè)率較高，虛警率和漏警率較為改進(jìn)方法均有明顯降低，總體明顯優(yōu)于未改進(jìn)方法提取的結(jié)果。

圖6 實(shí)驗(yàn)區(qū)1提取結(jié)果

圖7 實(shí)驗(yàn)區(qū)2提取結(jié)果

圖8 提取結(jié)果細(xì)節(jié)圖

表1 實(shí)驗(yàn)區(qū)精度評(píng)價(jià)表

5 結(jié)束語(yǔ)

本文提出了一種改進(jìn)變差函數(shù)的建筑區(qū)提取方法，該方法以建筑區(qū)及非建筑區(qū)的變差函數(shù)特征為依據(jù)，通過(guò)建筑區(qū)變差函數(shù)實(shí)驗(yàn)提取計(jì)算間距h，進(jìn)而確定窗口w，計(jì)算此參數(shù)下的變差函數(shù)值作為分類特征，采用EM算法實(shí)現(xiàn)建筑區(qū)和非建筑區(qū)的兩分類問(wèn)題。相比于傳統(tǒng)的變差函數(shù)計(jì)算方法，本文在計(jì)算方法上提出改進(jìn)：計(jì)算像素點(diǎn)（x，y）與固定方向（如0°方向）上間距h的點(diǎn)（x＋h，y）的半方差值，將此值作為點(diǎn)（x，y）的半方差值，然后計(jì)算窗口w內(nèi)的所有像素點(diǎn)的半方差值的中值作為變差函數(shù)值。結(jié)果表明，改進(jìn)的方法能夠更完整、更高效地提取建筑區(qū)。但是對(duì)于紋理特征不明顯，灰度值較低的建筑區(qū)，提取的效果不是很理想，下一步將嘗試引入其他特征對(duì)紋理特征不明顯的建筑區(qū)進(jìn)行提取算法研究。

［1]DELL’ACQUA F，GAMBA P.Texture－based characterization of urban environments on satellite SAR images［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41（1）：153－159.

［2]TISON C，NICOLAS J M，TUPIN F，et al.A new statistical model for Markovian classification of urban areas in highresolution SAR images［J］.IEEE Transactions on Geosicence and Remote Sensing，2004，42（10）：2046－2057.

［3]強(qiáng)永剛，殷建平，祝恩，等.基于小波變換和數(shù)學(xué)形態(tài)學(xué)的遙感圖像人工建筑區(qū)提?。跩］.中國(guó)圖象圖形學(xué)報(bào)，2008，13（8）：1459－1464.

［4]趙凌君，高貴，匡綱要.基于變差函數(shù)紋理特征的高分辨率SAR圖像建筑區(qū)提?。跩］.信號(hào)處理，2009，（25）9：1433－1442.

［5]金毅，潘懋，姚凌青，等.一種穩(wěn)健變差函數(shù)計(jì)算方法［J］.北京大學(xué)學(xué)報(bào)（自然科學(xué)版），2009（6）：1033－1038.

［6]MIRANDA F P，F(xiàn)ONSECA L E N，CARR J R，et al.Analysis of JERS－1（Fuyo－1）SAR data for vegetation discrimination in northwestern Brazil using the semivariogram textural classifier［J］.International Journal of Remote Sensing，1996，17（17）：3523－3529.

［7]MIRANDA F P，F(xiàn)ONSECA L E N，CARR J R.Semivariogram textural classification ofJERS－1（Fuyo－1）SAR data obtained over a flooded area of the Amazonrainforest［J］.International Journal of Remote Sensing，1998，19（3）：549－556.

［8]CURRAN P J.The semivariogram in remote sensing：An introduction［J］.Remote Sensing of Environment，1988（24）：493－507.

［9]OLIVER M，WEBSTER R，GERRARD J.Geostatistics in physical geography.Part I：Theory［J］.Transactions of the Institute of British Geographer，1989，41（3）：259－269.

［10]SHUFELT A.Performance evaluation and analysis of monocular building extraction from aerial imagery ［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1999，21（4）：311－326.