李劍南,魏永明,陳玉,高錦風(fēng)
(1 中國科學(xué)院空天信息創(chuàng)新研究院, 北京 100094; 2 中國科學(xué)院大學(xué)電子電氣與通信工程學(xué)院, 北京 100049)
地震裂縫是受地震活動影響,沿震源斷層錯動在地表形成的具有一定長度、寬度和深度的裂隙[1],是同震地表破裂的表現(xiàn)形式之一,發(fā)育于各類型地震同震地表破裂帶內(nèi)[2]。除地震活動本身造成的災(zāi)害外,地面永久變形和破裂產(chǎn)生的地震裂縫同樣會造成房屋等建筑物的開裂倒塌[3],嚴(yán)重威脅人民生命財產(chǎn)安全。同時,地震裂縫的空間分布特征及與其他破裂類型(地震鼓包、地震凹陷等)的組合形式是震源斷層活動性質(zhì)在地表最直接的體現(xiàn),查明地震裂縫分布規(guī)律,有助于深入理解震源斷層性質(zhì)、構(gòu)造運動狀態(tài)和地震破裂過程等基礎(chǔ)信息,在震后進(jìn)行發(fā)震構(gòu)造動力學(xué)、幾何學(xué)等相關(guān)研究。
受制于地形和自然條件等因素,地表破裂野外調(diào)查工作往往難以開展且存在一定的盲目性[4]。遙感技術(shù)具有宏觀性強、時效性高、覆蓋范圍廣等優(yōu)勢[5],在地表破裂調(diào)查研究中發(fā)揮了不可替代的作用。前人借助遙感手段,及時準(zhǔn)確地獲取了近20年7級以上地震同震地表破裂的空間位置并測量了同震位移量,為震災(zāi)評估、地質(zhì)調(diào)查提供了重要的科學(xué)依據(jù)[6-9]。然而,上述工作大都通過目視解譯完成,工作量大且主觀性較強,如何通過計算機結(jié)合專家知識實現(xiàn)自動、半自動提取,降低目視解譯工作量,提高提取精度與效率,成為目前國內(nèi)外亟需解決的重要問題。
諸多地表破裂類型中,地震裂縫結(jié)構(gòu)簡單,在遙感影像上通常表現(xiàn)為低亮度,紋理與光譜特征明顯[10-11],最易于提取。楊進(jìn)生等[12]結(jié)合多種卷積濾波器以及紋理信息提取方法對影像進(jìn)行增強,有效識別了華北平原地裂縫信息。王婭娟等[13]根據(jù)地裂縫的方向性特征,通過方向濾波有效識別提取了礦區(qū)地裂縫。然而上述方法僅僅停留在圖像增強處理后再進(jìn)行目視解譯的階段,并未實現(xiàn)自動、半自動提取。湯伏全等[14]利用最大似然法結(jié)合隨機森林算法,通過二次分類實現(xiàn)了采動地裂縫的精細(xì)提取。張興航等[10]利用Canny邊緣檢測算子提取線性特征,通過分形維數(shù)區(qū)分地裂縫與其他線性要素,實現(xiàn)了地裂縫的精細(xì)提取。受文獻(xiàn)[15]啟發(fā),Stumpf等[16]發(fā)現(xiàn)地裂縫的灰度剖面圖可近似擬合為高斯曲線,并通過高斯匹配濾波(matched filter with first-order derivative of Gaussian,MF-FDOG)算法進(jìn)行提取。韋博文等[17]將該方法應(yīng)用于西部某山體及礦區(qū)無人機影像,結(jié)果表明MF-FDOG算法對黃土地區(qū)地裂縫提取具有較強的適應(yīng)性。在前人研究中,受巖土體性質(zhì)及其強度影響,采動地裂縫、滑坡地裂縫的分布規(guī)律性不強,而震源斷層錯動形成的地震裂縫與發(fā)震斷裂帶位置基本重合,除光譜特征外,空間分布也具有明顯的方向性特征。本文選取昆侖山口西地震同震地表破裂帶為研究區(qū),引入國產(chǎn)高分二號(GF-2)作為影像數(shù)據(jù)源,對比分析目前地裂縫提取常用方法,建立了一套基于地震裂縫光譜、空間特征的地震裂縫提取方法,以期在震后第一時間快速獲取并分析地震裂縫空間分布及組合關(guān)系,指導(dǎo)地表破裂地面詳細(xì)調(diào)查、災(zāi)后應(yīng)急救援以及地震科學(xué)相關(guān)研究。
2001年11月14日,昆侖山口西發(fā)生MS8.1地震,震源斷層為東昆侖斷裂帶。該斷裂以左旋走滑為主要特征,同時兼有部分?jǐn)D壓分量,自第四紀(jì)以來活動強烈,曾多次發(fā)生7級以上地震[18-21]。地震形成的地表破裂帶西起庫水浣,東至卡巴紐爾多湖,全長達(dá)426 km,地表同震最大左旋位移量為7.6 m, 最大垂直位移量為4 m[2,22]。地表破裂延伸之長,規(guī)模之大,為世界之罕見。其中,構(gòu)造性地震裂縫在整個破裂帶上發(fā)育最為廣泛。此次地震由于地處藏北高海拔無人區(qū),氣候寒冷干燥,無植被覆蓋,破裂形跡至今大都得以完整保留,為地震裂縫提取研究提供了一個得天獨厚的試驗場。
研究區(qū)共收集3景原始GF-2數(shù)據(jù),覆蓋布喀達(dá)坂峰冰舌至紅水河口西破裂帶(表1)。以15 m分辨率Landsat ETM+ 波段8全色影像作為參考影像,結(jié)合有理多項式系數(shù)(rational polynomial coefficient,RPC)數(shù)據(jù)和90 m分辨率SRTM DEM,使用ENVI5.3流程化工具 RPC Orthorectification Using Reference Image進(jìn)行正射校正,該工具可自動從參考影像上尋找控制點并應(yīng)用于正射校正,極大提高了校正精度。采用Pansharp算法進(jìn)行圖像融合,融合后影像分辨率為0.8 m,色彩信息豐富,影紋特征清晰,完全滿足昆侖山口西地震同震地震裂縫的識別及提取(圖1)。
表1 研究區(qū)所選GF-2影像參數(shù)Table 1 Parameters of selected GF-2 image in study area
圖1 研究區(qū)GF-2遙感影像Fig.1 GF-2 remote sensing image of study area
基于地震裂縫與其他地物的光譜差異,先后采用自適應(yīng)直方圖均衡化和閾值分割對地震裂縫進(jìn)行粗提取,再通過數(shù)學(xué)形態(tài)學(xué)算法對地震裂縫孔隙進(jìn)行填充和連接,提取出地震裂縫完整形態(tài);最后根據(jù)地震裂縫的空間特征,排除其他地表信息(如陰影、紋溝、河床)的干擾,得到地震裂縫的精確提取結(jié)果。
光譜信息是遙感影像中對象的顏色及灰度或者波段間的亮度比等,是組成地物成分、結(jié)構(gòu)等屬性的直接反映,是遙感影像特征提取區(qū)別于另一種地物的本質(zhì)特征[23]。地震裂縫因其具有一定的深度而導(dǎo)致光線反射率低,在遙感影像中表現(xiàn)為低亮度值,可通過閾值分割將其完整形態(tài)提出(圖2)。同時也可以看出,地震裂縫兩側(cè)存在明顯的階躍型邊緣,傳統(tǒng)邊緣檢測算子如Sobel算子、Canny算子會對階躍型邊緣做出響應(yīng),檢測出偽邊緣,因此該方法在此并不適用[14,16]。
圖2 地震裂縫灰度剖面圖Fig.2 Gray profile of earthquake fissure
為避免光照等外界條件影響造成的明部或暗部細(xì)節(jié)丟失,捕捉更多的圖像細(xì)節(jié),首先采用限制對比度自適應(yīng)直方圖均衡化(contrast limited adaptive histgram equalization,CLAHE)對影像進(jìn)行增強處理[24]。均衡化之后裂縫特征更為明顯,在此基礎(chǔ)上進(jìn)行閾值分割,單閾值往往無法得到清晰有效的閾值分割圖像,同樣采用自適應(yīng)閾值(adaptive threshold)分割[25],通過高斯方程計算得到每個像素點的權(quán)重值,并使用該閾值對當(dāng)前像素點進(jìn)行處理,提取出地震裂縫完整形態(tài)。
當(dāng)?shù)卣鹆芽p較窄或由于部分被遮擋時,需要對提取結(jié)果進(jìn)行形態(tài)學(xué)濾波,以填充裂縫之間的間隙及裂縫內(nèi)的孔洞,使其形狀趨于規(guī)則。本研究使用形態(tài)學(xué)閉運算對裂縫孔隙進(jìn)行填充和連接。閉運算是先進(jìn)行膨脹再進(jìn)行腐蝕,具有平滑圖像、連接間隙、消除孔洞的作用[26]。
地震裂縫輪廓雖已被完整詳盡的提取,但提取結(jié)果中包含大量背景干擾信息(如噪聲點、紋溝、河床以及陰影等),需根據(jù)地震裂縫的空間特征實現(xiàn)精確提取。經(jīng)高分影像解譯,昆侖山口西地震裂縫主要由走向100°~110°的剪切裂縫、走向80°~90°的張剪裂縫以及走向45°~60°的張裂縫構(gòu)成(圖3)。不同性質(zhì)地震裂縫幾何形狀均表現(xiàn)為直線線性特征,空間組合呈平行排列、雁列狀斜列分布特征。因此,本文引入面積、長寬比和方向3個形狀特征進(jìn)行描述。
面積為對象區(qū)域內(nèi)所有的像元個數(shù)之和乘以單個像元對應(yīng)的地面面積,計算公式為
(1)
式中:A為單個像元對應(yīng)的地面面積,xkj和xki為第k行上對象邊界上的x坐標(biāo),且xkj大于xki;n為該對象的總行數(shù)。
長寬比為對象長軸與短軸的比值,計算公式為
(2)
式中:l1、l2分別為對象長軸和短軸,該參數(shù)反映對象的延伸性,可進(jìn)一步將裂縫與其他對象區(qū)分開來。
方向定義為對象長軸與水平方向的夾角,取值范圍為0°~180°。
不同類型裂縫的描述和閾值范圍如表2所示。地震裂縫表現(xiàn)為連續(xù)延伸,具有一定面積的8連通區(qū)域,因此可將小于面積閾值的區(qū)域判定為噪點去除。剪切裂縫平直延伸,其長寬比最大;張裂縫波折狀斷續(xù)延伸,長寬比最小,張剪裂縫受剪切和拉張共同作用;其長寬比介于上述二者之間。受構(gòu)造運動控制,地震裂縫延伸具有明顯的方向性特征,可以有效與紋溝、河床等線性特征進(jìn)行區(qū)分,河床雖可能與地震裂縫方向一致,但其長度應(yīng)大于地震裂縫,且具有地震裂縫所不具備的延伸性,可通過設(shè)定地震裂縫長寬比這一空間特征的閾值去除其干擾,提取結(jié)果如圖4所示。
在本研究中,CLCHE處理、自適應(yīng)閾值分割、形態(tài)學(xué)處理均基于 Python 語言開發(fā)環(huán)境實現(xiàn),基于空間特征的誤分類剔除通過ImageJ軟件實現(xiàn)[27]。
圖3 地震裂縫GF-2影像圖Fig.3 GF-2 image of earthquake fissure
表2 不同類型地震裂縫描述及空間特征閾值范圍Table 2 Characteristic description and spatial characteristic threshold range of different types of earthquake fissure
圖4 地震裂縫提取結(jié)果Fig.4 Extraction results of seismic fractures
本研究的提取結(jié)果體現(xiàn)了地震裂縫的空間分布與延伸特征。受制于地理位置和惡劣自然條件,地震裂縫精度評價很難開展野外實地驗證工作。為對提取精度進(jìn)行定量評估,將專家目視解譯圖作為真值與提取結(jié)果進(jìn)行比較,并與改進(jìn)的邊緣檢測[28]以及MF-FDOG算法[16]對比分析。根據(jù)檢測和參考圖中是否存在裂縫,將每個像素分配為裂縫或非裂縫,從提取結(jié)果和專家目視解譯圖兩組數(shù)據(jù)中計算真陽性率(true positive rate,RTP)和假陽性率(false positive rate,RFP)。
(3)
(4)
式中:TP-檢測為裂縫,實際為裂縫;FN-檢測為非裂縫,實際為裂縫;FP-檢測為裂縫,實際為非裂縫;TN-檢測為非裂縫,實際為非裂縫。
各方法提取結(jié)果如圖5所示,提取精度評價如表3所示,由于LoG-Canny邊緣檢測算法檢測出的是裂縫兩側(cè)邊緣,得到的結(jié)果中含有大量孔洞,采用閉運算時,較大尺寸的濾波核雖能將其完全填充,但會造成相鄰地震裂縫堆疊。同樣,該算法無法區(qū)別地裂縫與其他線性特征,導(dǎo)致部分河床的線性結(jié)構(gòu)與地裂縫相連,從而造成RFP較大。MF-FDOG算法則因為步驟復(fù)雜,且需要設(shè)置參數(shù)閾值較多,對于不同分辨率、不同類型的地震裂縫提取缺乏適應(yīng)性,提取結(jié)果未能反映地震裂縫的延伸性特征,從而造成RTP較小。而本研究方法提取精度最優(yōu),RTP最高,RFP最低,較好地反映了地震裂縫的形狀特征。
圖5 其他常用裂縫提取方法結(jié)果Fig.5 Results of other common fissure extraction methods
查明地震裂縫的空間分布特征及組合形式,有助于掌握分析震源斷層的最新幾何學(xué)、運動學(xué)和動力學(xué)特征,同時也能為理解里德爾剪切模型提供良好示例。
幾何上,里德爾剪切破裂[29]包括R 和共軛剪切破裂R′、 P破裂、Y破裂和共軛剪切X破裂、局部張性破裂(T)和局部擠壓結(jié)構(gòu)(圖6)。在走滑型地震中,當(dāng)?shù)乇韼r性松散、強度較低時,如洪積扇(臺地)、湖積臺地,震源斷層錯動常常發(fā)育純剪切裂縫(Y破裂);當(dāng)位于山前陡傾位置,受拉張力和剪切力共同作用,往往發(fā)育張剪裂縫(R破裂),其小角度相交于主剪切帶(交角一般小于15°),走滑方向與主剪切方向一致;張裂縫(T破裂)與主剪切帶大致呈45°交角,開口較大,平面上主要呈楔形或平行四邊形,呈雁列狀排列(圖3)。在調(diào)查同震地震裂縫的力學(xué)性質(zhì)、空間分布及延伸特征之后,將裂縫幾何形狀與里德爾剪切結(jié)構(gòu)相關(guān)聯(lián),可確定震源斷層的主剪切帶方向和區(qū)域應(yīng)力場[30],有助于進(jìn)行發(fā)震構(gòu)造動力學(xué)、幾何學(xué)等相關(guān)研究。
表3 地震裂縫提取精度評價Table 3 Evaluation of earthquake fissure extraction accuracy %
圖6 里德爾剪切破裂模式Fig.6 Comparison of Riedel shear pattern
經(jīng)遙感解譯和實地調(diào)查,地表破裂帶雖在宏觀上表現(xiàn)為一條直線,但往往由地震裂縫與其他破裂單元組合而成,地震裂縫呈羽列狀排列,右階階區(qū)發(fā)育表征地表縮短構(gòu)造,左階階區(qū)發(fā)育表征地表伸展構(gòu)造。
4.2.1 地震裂縫與伸展構(gòu)造
在地震破裂過程中,受兩條左行左階破裂相對運動,地表松散堆積或風(fēng)化殘坡積層上因局部拉張而發(fā)育伸展構(gòu)造,如張裂縫、地震凹陷和拉分盆地(圖7)。松散覆蓋層發(fā)育的拉分盆地從構(gòu)造上與黏土物質(zhì)中的雁列狀張裂縫類似,其規(guī)模與剪切位移呈正相關(guān)的關(guān)系[31-32]。演化模式如圖7(d),隨著變形的繼續(xù),雁列狀斜列的張裂縫發(fā)生疊覆,單個張裂縫或小型拉分盆地逐漸變寬、變深,最后聯(lián)合起來形成規(guī)模更大的盆地。圖7(c)中拉分盆地規(guī)模及形態(tài)發(fā)育程度高,可能對應(yīng)多次地震事件累計而成。
圖7 里德爾剪切盆地演化模式Fig.7 Evolution model of Riedel shear basin
4.2.2 地震裂縫與縮短構(gòu)造
在地震破裂過程中,受兩條左行右階破裂相對運動,在地表松散堆積或風(fēng)化殘坡積層上因局部壓扭形成縮短構(gòu)造,如地震鼓包、鼓梁(圖8)。前者為寬度略小于長度的小幅度隆起, 一般在平面上呈橢圓形狀,后者隆起長度遠(yuǎn)大于寬度,平面上呈長條狀。地震鼓梁長軸方向一般在N45 °~ 70 °W,在地表破裂帶各段落均有不同程度發(fā)育。張(剪)性裂縫段與壓扭性鼓包段相間排列,是整個破裂帶分布最為普遍的破裂形態(tài)組合[33]。相鄰鼓包呈左階斜列展布, 其間以右階斜列展布的張(剪)性裂縫首尾相連。前述圖3(a)中張剪裂縫呈波折狀不連續(xù)延伸,其階區(qū)位置實質(zhì)上也為地震鼓包(梁)相連。
地震裂縫的空間分布特征及組合形式反映地表破裂的幾何學(xué)、運動學(xué)特征,揭示了構(gòu)造活動下地震裂縫的動態(tài)發(fā)育過程。需要指出,地震裂縫的空間分布對于其他伸展、縮短型地表破裂的圈定具有重要的指示意義,后者雖然形態(tài)不一,無法通過遙感手段直接提取,但可根據(jù)地震裂縫階區(qū)間接圈劃出其可能發(fā)育位置。
圖8 張剪裂縫與縮短構(gòu)造的組合形式Fig.8 Combination of shear fracture and compression textures
震源斷層錯動形成的地震裂縫與發(fā)震斷裂帶位置基本重合,除光譜特征外,空間分布也具有明顯的方向性特征。本文選取昆侖山口西地震同震地表破裂帶為研究區(qū),建立了一套基于光譜、空間特征的地震裂縫提取方法。對比目前常用地裂縫提取方法,結(jié)果表明該方法結(jié)果最優(yōu),MF-FDOG算法次之,改進(jìn)的邊緣檢測算法效果較差。該方法適用于重點區(qū)域地震裂縫的精細(xì)刻畫描述,本文選取試驗區(qū)各類型地震裂縫形態(tài)完整,具備其典型的光譜特征和空間特征,在地表破裂提取方法研究和應(yīng)用中能起到一定的參考價值。如何在更為復(fù)雜的背景中提取地震裂縫,將是今后的研究重點。
運用里德爾剪切理論在遙感影像中識別地震裂縫,分析走滑型地震同震地表破裂的運動特征與空間展布關(guān)系。將裂縫幾何形狀及延伸方向與里德爾剪切結(jié)構(gòu)相關(guān)聯(lián),能快速確定主剪切帶方向和區(qū)域應(yīng)力場,有助于進(jìn)行發(fā)震構(gòu)造動力學(xué)、幾何學(xué)等相關(guān)研究。同時指出,地震裂縫的空間分布對于其他伸展、縮短型地表破裂的圈定具有重要的指示意義。
本研究引入國產(chǎn)高分衛(wèi)星GF-2作為影像數(shù)據(jù)源,顯示了國產(chǎn)高分辨率遙感數(shù)據(jù)在地表破裂考察等地震科學(xué)研究方面的應(yīng)用價值。提取結(jié)果證明,GF-2衛(wèi)星完全滿足大地震后地表破裂考察的需求,可為災(zāi)后應(yīng)急與災(zāi)情評估提供有效的空間信息依據(jù)。