趙志鵬,張 錦
(太原理工大學 礦業工程學院,山西 太原 030024)
我國是發生滑坡災害較多的國家之一,滑坡的發生常造成嚴重的人員傷亡和財產損失[1]。隨著空間遙感技術的發展,遙感圖像已經成為人們獲取滑坡信息的重要數據源。遙感影像滑坡的輪廓特征具有直觀性強,且不受光譜、紋理、地形和種類不同等因素的影響,因此準確快速地獲取滑坡輪廓特征,對應急救援和災后重建等工作具有重要意義。
國內外研究者對滑坡災害進行了多方面的研究,在滑坡特征提取方面,文獻[2]對滑坡影像利用Gabor變換來提取滑坡的紋理特征。文獻[3]從滑坡目標和背景顏色差異建立特征模型來提取滑坡。文獻[4-5]利用滑坡影像多種特征實現了滑坡自動識別提取。文獻[6-7]根據滑坡特征建立規則集實現滑坡的分層提取。顯然研究者對滑坡的輪廓特征大都是結合光譜和紋理等特征來共同識別提取滑坡,對于滑坡單一的輪廓特征在遙感影像上的表現形式和規律性研究不是很充分。
目標輪廓是圖像中用于識別的非常重要的一種特征,文獻[8]采用傅里葉描述子來提取光斑輪廓特征。文獻[9]利用傅里葉描述子來表達高分辨率遙感影像4種不同地物的形狀特征。文獻[10]利用人工神經網絡和傅里葉形狀描述子來識別造船部件。文獻[11]使用傅里葉描述子描述服裝輪廓特征。以上研究證明了傅里葉描述子提取目標輪廓特征的優越性,但針對提取滑坡遙感影像輪廓特征,傅里葉描述子是否具有實用性還有待研究。
本文提出基于傅里葉描述子的遙感影像滑坡輪廓特征提取方法,可以通過少量描述子重構滑坡輪廓。實驗結果表明,傅里葉描述子可以有效提取滑坡輪廓特征,計算量小,具有良好的普適性。
本文研究區位于滑坡災害頻發的四川省綿竹市北部山區,影像來源于谷歌地球。滑坡影像數據獲取日期為2014年12月23日,空間分辨率為0.25 m。
滑坡平面有舌形、橢圓形、長椅形、倒犁形、牛角形和平行四邊形等11種[12]形狀,本文實驗數據選取倒犁形、舌形和長條形3種不同形狀的滑坡高分辨率遙感影像。本次實驗采用MatlabR2018b編程實現,運行環境為Intel(R)-Core(TM)-i5-1135G7,2.00 GHz,內存8 GB,64位Window10操作系統。
在對滑坡輪廓特征提取之前,需要對滑坡圖像進行預處理。預處理的目的是消除圖像中與滑坡目標無關的信息,去除噪聲,最大限度地簡化數據,增強滑坡信息的可檢測性,提高滑坡輪廓提取精度。預處理的一般過程有灰度化、去噪、圖像平滑、形態學處理和邊緣檢測等。預處理過程如圖1所示。

圖1 滑坡遙感影像預處理流程Fig.1 Flow chart of landslide remote sensing image preprocessing
獲取滑坡輪廓的過程如圖2所示,圖2(a)為滑坡原始遙感影像,可以看出滑坡在影像上的色調與周圍環境有明顯的差別,呈灰白色,滑坡體表面幾乎沒有植被覆蓋。圖2(b)將滑坡遙感影像轉化為灰度圖像,利用中值濾波來去除圖像中的噪聲。圖2(c)選擇OTSU(最大類間方差)閾值分割法對滑坡圖像進行二值化,選擇最佳閾值使滑坡目標與背景進行分割,但同時受周圍地物的干擾,使滑坡目標與無關信息分離的不夠充分。圖2(d)中首先對分割后的圖像刪除小面積像素信息的干擾,其次對圖像進行孔洞填充,然后構造結構元素,進行形態學開運算,使滑坡目標與背景完全分離開。圖2(e)利用Canny邊緣檢測,獲得滑坡輪廓。

(a) 滑坡原始遙感影像

(a) 100個描述子
2.1.1 傅里葉描述子的原理
傅里葉描述子就是將圖像的輪廓從空間域轉換到頻域,提取頻域信息來表示目標形狀特征的方法。首先提取目標的邊界函數,然后將其展開后的系數幅度值作為目標輪廓的特征向量[13]。本文將遙感影像中滑坡的邊界信息通過傅里葉變換將滑坡的輪廓特征從空間域變換到頻域內,提取滑坡的頻域信息作為特征向量,將滑坡輪廓數字化,從而區分滑坡的不同輪廓,進而達到識別不同形狀滑坡的目的。
假設由K個點構成圖形邊界,從起始點(x0,y0)開始,沿逆時針方向追蹤至邊界點,將這些點的坐標當作復數來處理,可表示為:
z(k)=x(k)+y(k),k=0,1,2,…,K-1。
(1)
通過一維離散傅里葉變換,變換系數z(k)可表示為:
(2)
Z(u)稱為邊界的傅里葉描述子,通過傅里葉反變換可以重構z(k),表示為:

(3)
2.1.2 傅里葉描述子的歸一化
通過上述方法得到的傅里葉描述子為了使其具有平移、旋轉和尺度不變性,需要對其進行歸一化處理。根據傅里葉變換的性質,將目標輪廓邊界起始點位置平移α長度,放大r倍,旋轉角度φ和平移(x0,y0)后,新的傅里葉變換系數Z′(u)為:
Z′(u)=F[(x′+iy′)rejφ+(x0+iy0)]=
rejφF(x(l+α)+iy(l+α))+F(x0+iy0)=


(4)

u=0,1,2,…,K-1,
(5)

(6)
d(u)對于旋轉、位移、尺度和起始點位置不存在依賴關系,即得到歸一化傅里葉描述子,且具有平移、旋轉和尺度不變性,是一個魯棒性較好的圖像特征[14]。
2.1.3 傅里葉描述子提取滑坡輪廓特征
利用傅里葉變換得到滑坡邊界的傅里葉描述子,在計算傅里葉反變換進行圖像還原時,高頻分量決定滑坡細節部分,低頻分量決定滑坡總體輪廓。如果描述子過多就會造成特征冗余,增加計算量。隨著特征向量的減少,如果描述子過少則無法表達滑坡的總體輪廓特征。因此,需要選擇合適的描述子數量來提取滑坡的輪廓特征。
從圖3可以看出,用100個描述子重構得到的圖3(a)產生的輪廓與原始滑坡輪廓十分相似。圖3(b)顯示了50個描述子描述的結果,可以看出丟失了少量的細節信息,存在一些稍平滑的部分。由于高頻分量決定滑坡細節信息,低頻分量決定滑坡總體輪廓。隨著描述子的減少,描述的滑坡輪廓會丟失許多細節信息。圖3(c)顯示了30個描述子描述的結果,丟失了大量的細節信息,但保留了滑坡基本的細節特征。圖3(d)顯示了20個描述子描述的結果,保留了滑坡總體輪廓特征。圖3(e)是10個描述子描述的結果,已經產生了無法接受的失真,產生的邊界越來越平滑,丟失了關鍵的輪廓信息。因此,使用30個描述子可以基本描述滑坡輪廓的細節特征,使用20個描述子可以基本描述滑坡總體輪廓特征。

(a) 原始圖像

(a) 漏勺形滑坡
矩特征主要表征了圖像區域的幾何特征,又稱為幾何矩。Hu幾何矩具有平移、旋轉和尺度不變性,7個不變矩M1~M7是由二階和三階歸一化中心矩組成的[15-16]。在本文中,通過Hu不變矩對滑坡圖像計算得到的7個不變矩構成一組特征向量,可以驗證滑坡輪廓的平移、旋轉和尺度的不變性。其公式定義如下:
M1=η20+η02,
M2=(η20-η02)2+4η112,
M3=(η30-3η12)2+(3η21-η03)2,
M4=(η30+3η12)2+(3η21+η03)2,
M5=(η30-3η12)(η30+η12)[(η30+η12)2-3(η21+η03)2]+
(3η21-η03)(η21+η03)[3(η30+η12)2-(η21+η03)2],
M6=(η20-η02)[(η30+η12)2-(η21+η03)2]+
4η11(η30+3η12)(η21+η03),
M7=(3η21-η03)(η30+η12)[(η30+η12)2-3(η21+η03)2]+
(3η12-η30)(η03+η12)[3(η30+η12)2-(η21+η03)2]。
(7)
選取長條形滑坡為原始圖像,對滑坡進行平移、旋轉、鏡像、放大和縮放,驗證Hu不變矩描述滑坡輪廓具有平移、旋轉和縮放不變性,圖像變換如圖4所示。
通過Hu不變矩計算圖像變換后和另外2種形狀的滑坡輪廓形狀的矩不變量,其結果如表1所示。從表1可以看出,通過對滑坡圖像平移、旋轉和縮放變換后的Hu不變矩特征值與原始滑坡圖像的Hu不變矩特征值基本一致,同時也容易看出,原始滑坡圖像與不同形狀的滑坡圖像的Hu不變矩特征值具有較大差異,特別是M2,M4,M5,M6和M7可以明顯區別出不同形狀的滑坡。

表1 滑坡輪廓的矩不變量Tab.1 Moment invariants of landslide contour
通過傅里葉描述子和Hu不變矩提取滑坡輪廓特征后,選取支持向量機作為滑坡輪廓特征的分類器。本文采用LibSVM工具箱[17-18],對滑坡輪廓進行支持向量機多分類。首先,將描述滑坡輪廓特征的傅里葉描述子、Hu不變矩和滑坡不同形狀類別標簽輸入到支持向量機分類器中。然后,選擇線性核函數,其中Svmstrain參數選用“懲罰因子c=2,Gamma值g=0.07”。最后,將測試集的特征矩陣輸入到支持向量機分類器中得到滑坡形狀類別的分類準確率。
基于研究區的滑坡遙感影像,創建了一個小型滑坡遙感影像樣本庫,其中包含倒犁形、舌形和長條形滑坡影像各30幅。為提高運算速度,將影像尺寸壓縮為600 pixel×500 pixel,隨機抽取樣本庫中40%作為訓練集,剩余60%作為測試集,滑坡不同形狀的典型樣本如圖5所示。
為選取合適數量的傅里葉描述子描述滑坡輪廓,分別采用不同數量的傅里葉描述子(n取10,20,30,50,100)作為特征向量。利用LibSVM工具箱對樣本集進行分類和測試,其分類準確率如圖6所示。20個傅里葉描述子識別準確率最高,因此選擇特征向量長度為20來提取滑坡輪廓特征。

圖6 不同數量傅里葉描述子滑坡分類準確率Fig.6 Classification accuracy of sub-landslides described by different number of Fourier descriptors
由上文可知,選擇20個傅里葉描述子作為支持向量機分類器特征向量的輸入。融合特征是將傅里葉描述子的20個特征向量和Hu不變矩的7個特征向量采用串聯的方式進行融合。由于2種特征取值范圍存在較大的差異,且不同特征表達的意義也不同,因此融合之后對特征向量進行歸一化,將每一個特征向量歸一化至[0,1][19]。為了驗證傅里葉描述子可以很好地提取滑坡輪廓特征,對比分析了Hu不變矩和其二者的融合特征的識別效果。
對樣本集分別提取傅里葉描述子、Hu不變矩和二者的融合特征,利用支持向量機對樣本集進行分類訓練與測試,其識別結果如表2所示。結果顯示,3種方法消耗時間相近,Hu不變矩的分類準確率最低(87.04%),融合特征(傅里葉描述子和Hu不變矩)的分類準確率(90.74%)低于僅使用傅里葉描述子的識別效果(92.59%)。因此,傅里葉描述子可以有效地描述滑坡影像的輪廓特征,并具有較好的性能,故傅里葉描述子更適用于滑坡遙感影像輪廓特征提取。

表2 不同特征提取方法實驗對比Tab.2 Experimental comparison of different feature extraction methods
本文針對遙感影像不同形狀的滑坡,提出了一種基于傅里葉描述子的滑坡影像輪廓特征提取方法,實驗結果表明:
① 運用少量描述子就可以很好地描述滑坡遙感影像的輪廓特征,減少了數據冗余,提高了計算效率。
② 利用Hu不變矩驗證了滑坡輪廓具有平移、旋轉和尺度不變性,并且不同形狀滑坡之間的不變矩具有較大差異。
③ 采用多分類支持向量機對3種不同形狀的滑坡進行分類,通過對90幅滑坡影像的樣本庫進行實驗測試發現:傅里葉描述子比Hu不變矩和融合特征能更好地描述滑坡的輪廓特征。因此,傅里葉描述子是一種穩定有效、適用于滑坡遙感影像滑坡輪廓特征提取的方法。
④ 本文不足之處在于利用部分傅里葉描述子代替滑坡輪廓存在誤差,選擇多少個傅里葉描述子能夠表達滑坡的輪廓特征為經驗性值,缺乏有效的依據。優化傅里葉描述子數量重構滑坡輪廓和增加樣本庫容量提高分類精度是接下來需要繼續研究的內容。