余連生,李智勇,文貢堅,杜 春
國防科學技術大學電子科學與工程學院,湖南長沙410073
對于有潮海域,以海岸線和零米等深線為界,海岸線以上稱為潮上帶,零米等深線以下至30 m等深線處稱為潮下帶,兩者之間,即海岸線與零米等深線之間的地帶就稱為潮間帶。潮間帶地區是海洋環境監測、近海養殖業、風電場建設的重點區域,同時也是軍事研究的重要地帶,因此,對潮間帶地形的全面掌握具有重大的民用價值和軍事意義。由于大部分潮間帶地區無法直接測量或測量難度較大,因此潮間帶的地形資料非常缺乏,海圖上大部分潮間帶的地形資料甚至是空白的[1]。目前潮間帶地形的提取,主要是通過選取不同時期大量的不同潮位的衛星遙感圖像來提取潮間帶的遙感淹水線,結合數值模式進一步推算[1-2]。這種方法不但需要有大量不同時期的衛星影像支持,而且由于海平面容易受諸如潮汐、大氣壓變化、溫室效應以及大氣海洋的相互作用而變化[3],海岸線也受水蝕等因素影響而不斷變化,因此可能會帶來較大的偏差。同時推算的結果只是潮間帶的高程和大概輪廓,并不能真實顯示潮間帶的具體地形地物特征。
綜合利用不同平臺遙感衛星對有潮海域同一潮汐周期不同時相的成像圖像,能更好地獲取到潮間帶、海岸線等地形特征。其中,SAR可以全天時、全天候地觀測目標,便于在光學傳感器無法成像的時間段獲取到潮間帶的地貌特征,同時還可以從圖像中反映出不同地形地物的差別和地勢的高低。而全色遙感圖像具有較高的空間分辨率,可以獲取到海岸線附近地物目標的清晰輪廓信息,便于識別判讀,與SAR圖像之間可以實現較好的信息互補。
目前,SAR與全色圖像的融合算法主要都是基于小波變換域方法進行。文獻[4]提出的融合方法基于二維離散小波的Mallat算法,融合效果與加權融合方法的效果相當。文獻[5]提出用átrous小波分解算法,添加全色圖像的細節特征到SAR圖像中去,但完全忽略了全色圖像的低頻信息。文獻[6—7]分別利用脊波變換(ridgelet)和曲波變換(curvelet)實現了SAR圖像與全色圖像的融合,采取低頻高頻都取大的融合規則,在保持全色圖像信息的同時,添加了SAR圖像中的重要目標信息(特別亮的點或區域),但是這會遮蓋一些比較暗的邊緣、紋理等信息。
針對上述算法存在的不足,本文將平移不變離散小波變換(shift-invariant discrete wavelet transform,SIDWT)算法和新的融合規則應用于同一潮汐周期不同時相的海岸帶地區SAR與全色遙感圖像融合,得到了既能突出潮間帶形狀,又能突出海岸線特征的更高質量的圖像,便于目視解譯及后續處理,取得了較好的結果。
針對本文的實際應用來說,由于待融合的SAR和光學圖像等光譜差異較大,而且圖像中地形差異明顯,圖像校正、配準難度較大,因此采取的融合方法是否具有平移不變性對融合結果影響較大[8],容易在融合結果中引入虛假信息,如振鈴和混疊效應。
為了克服Mallat算法不具有平移不變性的不足,對圖像進行小波變換時,需要將其分解為平移不變的小波表示,最直接的方法是計算輸入圖像在所有可能位移上的小波變換,如àtrous小波變換,但是其計算復雜度很高,而且并不必要[9]。文獻[10]提出了一種過完備的小波分解方法——離散小波框架(discrete wavelet frame,DWF),來解決這一問題?;诖朔椒?文獻[11]提出了平移不變離散小波變換(SIDWT)算法。
針對一維信號的情況,文獻[11]中SIDWT的每一步運算都將輸入信號分離為小波序列wi(n)和尺度序列si(n)

式中,si(n)作為下一層分解的輸入信號,s0(n)= f(n),f(n)為輸入信號。第 i層的低通濾波器h(2ik)和高通濾波器g(2ik)是通過在濾波器原型h(k)和 g(k)的序列間插入適量的0值得到的。其中,h(k)和 g(k)滿足

式中,H(z),G(z)分別是 h(k),g(k)的 z變換。通過上述變換,構成了SIDWT的分解。對比標準的離散小波變換,由于舍棄了下采樣,SIDWT算法導致了大量的冗余小波表示。因此,利用SIDWT算法對源圖像進行 n層分解可得到一幅逼近圖像(低頻部分)和3n幅代表不同分辨率下細節信息的細節圖像(高頻部分),且每幅圖像的尺寸大小均和源圖像一樣。
SIDWT的反變換(inverse SIDWT)是將上面得到的具有平移不變性的小波序列和尺度序列分別與相應的重建濾波器~h(2ik)和~g(2ik)進行卷積

通過SIDWT的反變換,從而實現輸入信號f(n)的重建。對于二維圖像信號,其分解過程是上述一維信號流程分別在圖像的行與列上所進行的處理過程

式中,s0(m,n)=f(m,n),f(m,n)為輸入信號。描述二維圖像每一層分解后的水平高頻分量、垂直高頻分量、對角高頻分量,各自對應相應三個方向的小波序列。
綜上所述,利用SIDWT算法進行融合的過程與利用離散小波變換相同:輸入圖像被分解成具有平移不變性的小波表示,然后通過一定融合規則得到融合后的小波表示,最后通過SIDWT的反變換得到最終融合圖像。
本文將平移不變離散小波變換算法用于SAR與全色圖像的融合,相應的小波基選擇haar小波。這主要是因為haar小波計算簡單,而且它是一切具有緊支集的規范正交小波基中唯一具有對稱性(或反對稱性)的[12],具有平移不變性。
相對全色圖像來說,SAR圖像的整體圖像亮度較低,但對比度也較好。針對實際應用,本文的融合目的是既保留SAR與光學圖像的邊緣、輪廓、紋理等信息,又能在融合圖像中添加SAR圖像和光學圖像所特有的地物特征描述,從而能夠更好地確定潮間帶特征,以便于后續的圖像處理。
圖像經過SIDWT分解后,得到的低頻逼近圖像包含了圖像的基本輪廓、形狀等大部分信息,高頻部分則包含圖像的邊緣細節信息。結合融合目的,本文主要針對低頻部分進行處理。在參考現有的相關算法[13-15]基礎上,本文提出一種基于邊緣信息的加權融合規則,具體實現如下:
首先,使用ROA算子對SAR圖像進行邊緣提取,用Canny算子對全色圖像進行邊緣提取,分別得到SAR與全色圖像的邊緣圖像。然后,在下列四種可能存在的情況下進行逐像素點的比較:
(1)該像素點只是SAR圖像的邊緣點;
(2)該像素點只是全色圖像PAN的邊緣點;
(3)該像素點同時是兩圖像的邊緣點;
(4)該像素點都不是兩圖像的邊緣點。
針對不同的情況,對SAR與全色圖像經過SIDWT分解后的低頻部分進行逐像素值的加權。設融合圖像的低頻圖像像素值為

式中,(i,j)為圖像中對應位置的像素值坐標;ω1, ω2分別為SAR與全色圖像中的權值,滿足ω1+ ω2=1。
依據實際應用,為了更好地突出邊緣細節信息,并保持較好的圖像對比度,基于如下原則選擇權值:
對于情況1,為更好地添加SAR圖像的特征信息,則SAR圖像所占的比例較大,選擇ω1= 0.75,ω2=0.25;
對于情況2,相應地,全色圖像所占的比例較大,選擇ω1=0.25,ω2=0.75;
對于情況3,則兩者所占的比例相同,選擇ω1=0.5,ω2=0.5;
對于情況4,由于一幅圖像中非邊緣點個數占絕大部分,為了能更好地保持SAR圖像的特征,取ω1=0.75,ω2=0.25。
由于圖像的邊緣細節信息表現為SIDWT分解后高頻部分中的絕對值較大的系數,因此,對于高頻部分,采用基于絕對值最大準則進行融合。
得到新的低頻部分和高頻部分后,再通過平移不變離散小波反變換,最終得到融合結果。算法流程如圖1所示。

圖1 算法流程圖Fig.1 Flow chart of the proposed method
本文試驗數據采用配準好的2007年10月18日拍攝的某島嶼的SAR和全色遙感圖像,如圖2所示。結合兩幅圖像的成像時間和該地區當天的氣象資料可知,SAR圖像中呈現的是晚上退潮后的情景,而全色圖像則是反映了港口上午漲潮時刻情形。從圖中可以明顯看出該地區在不同時刻的地形差異。

圖2 原始SAR與全色圖像Fig.2 Original SAR and panchromatic images
圖3給出了基于不同方法的融合結果。作為對比,圖3(a)是像素等比加權融合算法得到的融合結果,圖3(b)為利用Mallat小波變換方法(小波基為Bior2.2小波)及本文融合規則得到的融合結果,圖3(c)為利用àtrous小波變換方法[4]得到的融合結果,圖3(d)為利用本文算法得到的融合結果。本試驗中,后三種算法中的小波變換分解層數都是三層。

圖3 融合結果圖Fig.3 Fusion results
可以看出,各種算法都體現了潮間帶和近岸地區的地形地貌,但是本文算法結果亮度較高,對比度也較好,潮間帶、潮上帶分界明顯,目視解譯可以很好地區分兩者,同時還較好地實現了SAR圖像與全色圖像的地物信息互補,可以為繪制地圖提供圖像支持。其中,由于àtrous小波變換算法完全忽略了全色圖像的低頻信息,導致融合結果中光學特征丟失嚴重。
為了更客觀地評價融合效果,本文采用信息熵、平均梯度、標準差等指標對融合圖像進行評價。信息熵是衡量圖像信息豐富程度的重要指標,信息熵越大,說明信息量增加得越多;平均梯度反映了圖像中的微小細節反差和紋理變化等結構特征,同時也反映了圖像的清晰度,其值越大,表示圖像的清晰度越好;標準差反映了圖像灰度值相對于平均值的離散狀況,標準差越大,圖像灰度級分布越分散,反差越大。
表1給出了不同融合算法的客觀指標評價結果??梢钥闯?本文算法得到的圖像熵值最大,清晰度也最高,標準差也最大,這與目視評價一致。綜合所有評價指標,說明本文算法的效果最好。
為了消除潮汐和天氣等因素的不良影響,更好地探測港口近岸地區及潮間帶地區的地形地物等全面狀況,本文利用同一天不同時相獲取的同一地區的SAR與全色圖像,基于兩者之間較強的互補性,采用平移不變離散小波變換(SIDWT)算法和新的融合規則實現圖像融合。從融合結果看,潮間帶地形和近岸地區的地貌等較好地得到了綜合、增強,在后續的目視解譯中,更容易區分不同地形的特征,可以為繪制海圖和電子海圖提供支持,具有一定的應用意義。

表1 不同算法的客觀指標評價結果Tab.1 Quantitative indices values of fusion results
[1] ZHAN G Wenjing,ZHU Shouxian,HUAN G Weigen. Application ofSatellite Remote Sensing Data to Storm Surge Flooding Calculation[J].Journal of PLA University of Science and Technology:Natural Science Edition,2009, 10(5):501-506.(張文靜,朱首賢,黃韋艮.衛星遙感資料在湛江港風暴潮漫灘計算中的應用[J].解放軍理工大學學報:自然科學版,2009,10(5):501-506.)
[2] ZHU Shouxian,HUANG Weigen,DING Pingxing,et al. Application of Satellite Remote Sensing Data to the Calculation ofTopography and Flooding at Littoral Zone[J]. Journal of East China Normal University:Natural Science, 2007(4):1-8.(朱首賢,黃韋艮,丁平興,等.衛星遙感資料在潮間帶地形和漫灘計算中的應用[J].華東師范大學學報:自然科學版,2007(4):1-8.)
[3] ZHAN Jingang,WANG Yong,XU Houze,et al.The Wavelet Analysis of Sea Level Change in China Sea during 1992~2006[J].Acta Geodaetica et Cartographica Sinica, 2008,37(4):438-443.(詹金剛,王勇,許厚澤,等.我國近海1992~2006年海平面變化的小波分析[J].測繪學報,2008,37(4):438-443.)
[4] SONG Jianshe,ZHENG Yongan,LIU Yingchun.Image Fusion Algorithm of SAR and Optical Images Based on Wavelet Transform[J].Application Reasearch of Computers, 2004,21(10):110-113.(宋建社,鄭永安,劉迎春.基于小波變換的SAR與可見光圖像融合算法[J].計算機應用研究,2004,21(10):110-113.)
[5] CHIBANI Y.Selective Synthetic Aperture Radar and Panchromatic Image Fusion by Using theàtrous Wavelet Decomposition[J].EURASIP Journal on Applied Signal Processing,2005,14:2207-2214.
[6] LI Huihui,GUO Lei,LI Guoxin.Is Ridgelet Transform Better than Wavelet Transform in SAR and Optical Image Fusion?[J].Journalof Northwestern Polytechnical University,2006,24(4):418-422.(李暉暉,郭雷,李國新.基于脊波變換的SAR與可見光圖像融合研究[J].西北工業大學學報,2006,24(4):418-422.)
[7] L I Huihui,GUO Lei,L IU Hang.Research on Image Fusion Based on the Second Generation Curvelet Transform [J].Acta Optica Sinica,2006,26(5):657-662.(李暉暉,郭雷,劉航.基于二代Curvelet變換的圖像融合研究[J].光學學報,2006,26(5):657-662.)
[8] L IU Kun,GU O L ei,L I H uihui,et al.Image Fusion Algorithm Using Stationary Wavelet Transform[J]. Computer Engineering and Applications,2007,43(12):59-61.(劉坤,郭雷,李暉暉,等.一種基于靜態小波變換的圖像融合算法[J].計算機工程與應用,2007,43(12):59-61.)
[9] BEYKIN G.On the Representation of Operators in Bases of Compactly Supported Wavelets[J].SIAM Journal on Numerical Analysis,1992,29(6):1716-1740.
[10] UNSER M.Texture Classification and Segmentation Using Wavelet Frames[J].IEEE Transactions on Image Processing,1995,4(11):1549-1560.
[11] ROCKINGER O.Image Sequence Fusion Using a Shiftinvariant WaveletTransform[C]∥Proceedings of the IEEE InternationalConference on Image Processing. Santa Barbara,CA:IEEE,1997:288-291.
[12] DAUBECHIES I.Ten Lecture on Wavelets[M]. Philadelphia:Capital City Press,1992.
[13] L I H uihui.Study on M ultisensor Images Fusion Algorithm[D]. Xi’an:Northwestern Polytechnical University,2006.(李暉暉.多傳感器圖像融合算法研究[D].西安:西北工業大學,2006.)
[14] WU Zhiguo,WAN G Yanjie.Image Fusion Algorithm Using Curvelet Transform Based on the Edge Detection [J].Optical Technique,2009,35(5):682-686.(武治國,王延杰.一種基于邊緣的曲波變換圖像融合方法[J].光學技術,2009,35(5):682-686.)
[15] ZHOU Liya,QIN Zhiyuan.Research on Enhancing Edge Feature Image Fusion Method[J].Journal of Geomatics Science and Technology,2008,25(6):448-450.(周麗雅,秦志遠.突出邊緣特征的圖像融合方法研究[J].測繪科學技術學報,2008,25(6):448-450).