李劍南,魏永明,陳玉,高錦風(fēng)

(1 中國(guó)科學(xué)院空天信息創(chuàng)新研究院, 北京 100094; 2 中國(guó)科學(xué)院大學(xué)電子電氣與通信工程學(xué)院, 北京 100049)

地震裂縫是受地震活動(dòng)影響,沿震源斷層錯(cuò)動(dòng)在地表形成的具有一定長(zhǎng)度、寬度和深度的裂隙[1],是同震地表破裂的表現(xiàn)形式之一,發(fā)育于各類型地震同震地表破裂帶內(nèi)[2]。除地震活動(dòng)本身造成的災(zāi)害外,地面永久變形和破裂產(chǎn)生的地震裂縫同樣會(huì)造成房屋等建筑物的開裂倒塌[3],嚴(yán)重威脅人民生命財(cái)產(chǎn)安全。同時(shí),地震裂縫的空間分布特征及與其他破裂類型(地震鼓包、地震凹陷等)的組合形式是震源斷層活動(dòng)性質(zhì)在地表最直接的體現(xiàn),查明地震裂縫分布規(guī)律,有助于深入理解震源斷層性質(zhì)、構(gòu)造運(yùn)動(dòng)狀態(tài)和地震破裂過程等基礎(chǔ)信息,在震后進(jìn)行發(fā)震構(gòu)造動(dòng)力學(xué)、幾何學(xué)等相關(guān)研究。

受制于地形和自然條件等因素,地表破裂野外調(diào)查工作往往難以開展且存在一定的盲目性[4]。遙感技術(shù)具有宏觀性強(qiáng)、時(shí)效性高、覆蓋范圍廣等優(yōu)勢(shì)[5],在地表破裂調(diào)查研究中發(fā)揮了不可替代的作用。前人借助遙感手段,及時(shí)準(zhǔn)確地獲取了近20年7級(jí)以上地震同震地表破裂的空間位置并測(cè)量了同震位移量,為震災(zāi)評(píng)估、地質(zhì)調(diào)查提供了重要的科學(xué)依據(jù)[6-9]。然而,上述工作大都通過目視解譯完成,工作量大且主觀性較強(qiáng),如何通過計(jì)算機(jī)結(jié)合專家知識(shí)實(shí)現(xiàn)自動(dòng)、半自動(dòng)提取,降低目視解譯工作量,提高提取精度與效率,成為目前國(guó)內(nèi)外亟需解決的重要問題。

諸多地表破裂類型中,地震裂縫結(jié)構(gòu)簡(jiǎn)單,在遙感影像上通常表現(xiàn)為低亮度,紋理與光譜特征明顯[10-11],最易于提取。楊進(jìn)生等[12]結(jié)合多種卷積濾波器以及紋理信息提取方法對(duì)影像進(jìn)行增強(qiáng),有效識(shí)別了華北平原地裂縫信息。王婭娟等[13]根據(jù)地裂縫的方向性特征,通過方向?yàn)V波有效識(shí)別提取了礦區(qū)地裂縫。然而上述方法僅僅停留在圖像增強(qiáng)處理后再進(jìn)行目視解譯的階段,并未實(shí)現(xiàn)自動(dòng)、半自動(dòng)提取。湯伏全等[14]利用最大似然法結(jié)合隨機(jī)森林算法,通過二次分類實(shí)現(xiàn)了采動(dòng)地裂縫的精細(xì)提取。張興航等[10]利用Canny邊緣檢測(cè)算子提取線性特征,通過分形維數(shù)區(qū)分地裂縫與其他線性要素,實(shí)現(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ū)無人機(jī)影像,結(jié)果表明MF-FDOG算法對(duì)黃土地區(qū)地裂縫提取具有較強(qiáng)的適應(yīng)性。在前人研究中,受巖土體性質(zhì)及其強(qiáng)度影響,采動(dòng)地裂縫、滑坡地裂縫的分布規(guī)律性不強(qiáng),而震源斷層錯(cuò)動(dòng)形成的地震裂縫與發(fā)震斷裂帶位置基本重合,除光譜特征外,空間分布也具有明顯的方向性特征。本文選取昆侖山口西地震同震地表破裂帶為研究區(qū),引入國(guó)產(chǎn)高分二號(hào)(GF-2)作為影像數(shù)據(jù)源,對(duì)比分析目前地裂縫提取常用方法,建立了一套基于地震裂縫光譜、空間特征的地震裂縫提取方法,以期在震后第一時(shí)間快速獲取并分析地震裂縫空間分布及組合關(guān)系,指導(dǎo)地表破裂地面詳細(xì)調(diào)查、災(zāi)后應(yīng)急救援以及地震科學(xué)相關(guān)研究。

1 研究區(qū)概況以及數(shù)據(jù)源

1.1 東昆侖斷裂帶

2001年11月14日,昆侖山口西發(fā)生MS8.1地震,震源斷層為東昆侖斷裂帶。該斷裂以左旋走滑為主要特征,同時(shí)兼有部分?jǐn)D壓分量,自第四紀(jì)以來活動(dòng)強(qiáng)烈,曾多次發(fā)生7級(jí)以上地震[18-21]。地震形成的地表破裂帶西起庫(kù)水浣,東至卡巴紐爾多湖,全長(zhǎng)達(dá)426 km,地表同震最大左旋位移量為7.6 m, 最大垂直位移量為4 m[2,22]。地表破裂延伸之長(zhǎng),規(guī)模之大,為世界之罕見。其中,構(gòu)造性地震裂縫在整個(gè)破裂帶上發(fā)育最為廣泛。此次地震由于地處藏北高海拔無人區(qū),氣候寒冷干燥,無植被覆蓋,破裂形跡至今大都得以完整保留,為地震裂縫提取研究提供了一個(gè)得天獨(dú)厚的試驗(yàn)場(chǎng)。

1.2 數(shù)據(jù)源

研究區(qū)共收集3景原始GF-2數(shù)據(jù),覆蓋布喀達(dá)坂峰冰舌至紅水河口西破裂帶(表1)。以15 m分辨率Landsat ETM+ 波段8全色影像作為參考影像,結(jié)合有理多項(xiàng)式系數(shù)(rational polynomial coefficient,RPC)數(shù)據(jù)和90 m分辨率SRTM DEM,使用ENVI5.3流程化工具 RPC Orthorectification Using Reference Image進(jìn)行正射校正,該工具可自動(dòng)從參考影像上尋找控制點(diǎn)并應(yīng)用于正射校正,極大提高了校正精度。采用Pansharp算法進(jìn)行圖像融合,融合后影像分辨率為0.8 m,色彩信息豐富,影紋特征清晰,完全滿足昆侖山口西地震同震地震裂縫的識(shí)別及提取(圖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

2 地震裂縫提取方法

基于地震裂縫與其他地物的光譜差異,先后采用自適應(yīng)直方圖均衡化和閾值分割對(duì)地震裂縫進(jìn)行粗提取,再通過數(shù)學(xué)形態(tài)學(xué)算法對(duì)地震裂縫孔隙進(jìn)行填充和連接,提取出地震裂縫完整形態(tài);最后根據(jù)地震裂縫的空間特征,排除其他地表信息(如陰影、紋溝、河床)的干擾,得到地震裂縫的精確提取結(jié)果。

2.1 地震裂縫粗提取

光譜信息是遙感影像中對(duì)象的顏色及灰度或者波段間的亮度比等,是組成地物成分、結(jié)構(gòu)等屬性的直接反映,是遙感影像特征提取區(qū)別于另一種地物的本質(zhì)特征[23]。地震裂縫因其具有一定的深度而導(dǎo)致光線反射率低,在遙感影像中表現(xiàn)為低亮度值,可通過閾值分割將其完整形態(tài)提出(圖2)。同時(shí)也可以看出,地震裂縫兩側(cè)存在明顯的階躍型邊緣,傳統(tǒng)邊緣檢測(cè)算子如Sobel算子、Canny算子會(huì)對(duì)階躍型邊緣做出響應(yīng),檢測(cè)出偽邊緣,因此該方法在此并不適用[14,16]。

圖2 地震裂縫灰度剖面圖Fig.2 Gray profile of earthquake fissure

為避免光照等外界條件影響造成的明部或暗部細(xì)節(jié)丟失,捕捉更多的圖像細(xì)節(jié),首先采用限制對(duì)比度自適應(yīng)直方圖均衡化(contrast limited adaptive histgram equalization,CLAHE)對(duì)影像進(jìn)行增強(qiáng)處理[24]。均衡化之后裂縫特征更為明顯,在此基礎(chǔ)上進(jìn)行閾值分割,單閾值往往無法得到清晰有效的閾值分割圖像,同樣采用自適應(yīng)閾值(adaptive threshold)分割[25],通過高斯方程計(jì)算得到每個(gè)像素點(diǎn)的權(quán)重值,并使用該閾值對(duì)當(dāng)前像素點(diǎn)進(jìn)行處理,提取出地震裂縫完整形態(tài)。

2.2 地震裂縫精提取

當(dāng)?shù)卣鹆芽p較窄或由于部分被遮擋時(shí),需要對(duì)提取結(jié)果進(jìn)行形態(tài)學(xué)濾波,以填充裂縫之間的間隙及裂縫內(nèi)的孔洞,使其形狀趨于規(guī)則。本研究使用形態(tài)學(xué)閉運(yùn)算對(duì)裂縫孔隙進(jìn)行填充和連接。閉運(yùn)算是先進(jìn)行膨脹再進(jìn)行腐蝕,具有平滑圖像、連接間隙、消除孔洞的作用[26]。

地震裂縫輪廓雖已被完整詳盡的提取,但提取結(jié)果中包含大量背景干擾信息(如噪聲點(diǎn)、紋溝、河床以及陰影等),需根據(jù)地震裂縫的空間特征實(shí)現(xiàn)精確提取。經(jīng)高分影像解譯,昆侖山口西地震裂縫主要由走向100°～110°的剪切裂縫、走向80°～90°的張剪裂縫以及走向45°～60°的張裂縫構(gòu)成(圖3)。不同性質(zhì)地震裂縫幾何形狀均表現(xiàn)為直線線性特征,空間組合呈平行排列、雁列狀斜列分布特征。因此,本文引入面積、長(zhǎng)寬比和方向3個(gè)形狀特征進(jìn)行描述。

面積為對(duì)象區(qū)域內(nèi)所有的像元個(gè)數(shù)之和乘以單個(gè)像元對(duì)應(yīng)的地面面積,計(jì)算公式為

(1)

式中:A為單個(gè)像元對(duì)應(yīng)的地面面積,xkj和xki為第k行上對(duì)象邊界上的x坐標(biāo),且xkj大于xki;n為該對(duì)象的總行數(shù)。

長(zhǎng)寬比為對(duì)象長(zhǎng)軸與短軸的比值,計(jì)算公式為

(2)

式中:l1、l2分別為對(duì)象長(zhǎng)軸和短軸,該參數(shù)反映對(duì)象的延伸性,可進(jìn)一步將裂縫與其他對(duì)象區(qū)分開來。

方向定義為對(duì)象長(zhǎng)軸與水平方向的夾角,取值范圍為0°～180°。

不同類型裂縫的描述和閾值范圍如表2所示。地震裂縫表現(xiàn)為連續(xù)延伸,具有一定面積的8連通區(qū)域,因此可將小于面積閾值的區(qū)域判定為噪點(diǎn)去除。剪切裂縫平直延伸,其長(zhǎng)寬比最大;張裂縫波折狀斷續(xù)延伸,長(zhǎng)寬比最小,張剪裂縫受剪切和拉張共同作用;其長(zhǎng)寬比介于上述二者之間。受構(gòu)造運(yùn)動(dòng)控制,地震裂縫延伸具有明顯的方向性特征,可以有效與紋溝、河床等線性特征進(jìn)行區(qū)分,河床雖可能與地震裂縫方向一致,但其長(zhǎng)度應(yīng)大于地震裂縫,且具有地震裂縫所不具備的延伸性,可通過設(shè)定地震裂縫長(zhǎng)寬比這一空間特征的閾值去除其干擾,提取結(jié)果如圖4所示。

在本研究中,CLCHE處理、自適應(yīng)閾值分割、形態(tài)學(xué)處理均基于 Python 語(yǔ)言開發(fā)環(huán)境實(shí)現(xiàn),基于空間特征的誤分類剔除通過ImageJ軟件實(shí)現(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

3 地震裂縫提取精度評(píng)價(jià)

本研究的提取結(jié)果體現(xiàn)了地震裂縫的空間分布與延伸特征。受制于地理位置和惡劣自然條件,地震裂縫精度評(píng)價(jià)很難開展野外實(shí)地驗(yàn)證工作。為對(duì)提取精度進(jìn)行定量評(píng)估,將專家目視解譯圖作為真值與提取結(jié)果進(jìn)行比較,并與改進(jìn)的邊緣檢測(cè)[28]以及MF-FDOG算法[16]對(duì)比分析。根據(jù)檢測(cè)和參考圖中是否存在裂縫,將每個(gè)像素分配為裂縫或非裂縫,從提取結(jié)果和專家目視解譯圖兩組數(shù)據(jù)中計(jì)算真陽(yáng)性率(true positive rate,RTP)和假陽(yáng)性率(false positive rate,RFP)。

(3)

(4)

式中:TP-檢測(cè)為裂縫,實(shí)際為裂縫;FN-檢測(cè)為非裂縫,實(shí)際為裂縫;FP-檢測(cè)為裂縫,實(shí)際為非裂縫;TN-檢測(cè)為非裂縫,實(shí)際為非裂縫。

各方法提取結(jié)果如圖5所示,提取精度評(píng)價(jià)如表3所示,由于LoG-Canny邊緣檢測(cè)算法檢測(cè)出的是裂縫兩側(cè)邊緣,得到的結(jié)果中含有大量孔洞,采用閉運(yùn)算時(shí),較大尺寸的濾波核雖能將其完全填充,但會(huì)造成相鄰地震裂縫堆疊。同樣,該算法無法區(qū)別地裂縫與其他線性特征,導(dǎo)致部分河床的線性結(jié)構(gòu)與地裂縫相連,從而造成RFP較大。MF-FDOG算法則因?yàn)椴襟E復(fù)雜,且需要設(shè)置參數(shù)閾值較多,對(duì)于不同分辨率、不同類型的地震裂縫提取缺乏適應(yīng)性,提取結(jié)果未能反映地震裂縫的延伸性特征,從而造成RTP較小。而本研究方法提取精度最優(yōu),RTP最高,RFP最低,較好地反映了地震裂縫的形狀特征。

圖5 其他常用裂縫提取方法結(jié)果Fig.5 Results of other common fissure extraction methods

4 地震裂縫空間分布特征及組合形式

查明地震裂縫的空間分布特征及組合形式,有助于掌握分析震源斷層的最新幾何學(xué)、運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征,同時(shí)也能為理解里德爾剪切模型提供良好示例。

4.1 分布特征

幾何上,里德爾剪切破裂[29]包括R 和共軛剪切破裂R′、 P破裂、Y破裂和共軛剪切X破裂、局部張性破裂(T)和局部擠壓結(jié)構(gòu)(圖6)。在走滑型地震中,當(dāng)?shù)乇韼r性松散、強(qiáng)度較低時(shí),如洪積扇(臺(tái)地)、湖積臺(tái)地,震源斷層錯(cuò)動(dòng)常常發(fā)育純剪切裂縫(Y破裂);當(dāng)位于山前陡傾位置,受拉張力和剪切力共同作用,往往發(fā)育張剪裂縫(R破裂),其小角度相交于主剪切帶(交角一般小于15°),走滑方向與主剪切方向一致;張裂縫(T破裂)與主剪切帶大致呈45°交角,開口較大,平面上主要呈楔形或平行四邊形,呈雁列狀排列(圖3)。在調(diào)查同震地震裂縫的力學(xué)性質(zhì)、空間分布及延伸特征之后,將裂縫幾何形狀與里德爾剪切結(jié)構(gòu)相關(guān)聯(lián),可確定震源斷層的主剪切帶方向和區(qū)域應(yīng)力場(chǎng)[30],有助于進(jìn)行發(fā)震構(gòu)造動(dòng)力學(xué)、幾何學(xué)等相關(guān)研究。

表3 地震裂縫提取精度評(píng)價(jià)Table 3 Evaluation of earthquake fissure extraction accuracy %

圖6 里德爾剪切破裂模式Fig.6 Comparison of Riedel shear pattern

4.2 組合形式

經(jīng)遙感解譯和實(shí)地調(diào)查,地表破裂帶雖在宏觀上表現(xiàn)為一條直線,但往往由地震裂縫與其他破裂單元組合而成,地震裂縫呈羽列狀排列,右階階區(qū)發(fā)育表征地表縮短構(gòu)造,左階階區(qū)發(fā)育表征地表伸展構(gòu)造。

4.2.1 地震裂縫與伸展構(gòu)造

在地震破裂過程中,受兩條左行左階破裂相對(duì)運(yùn)動(dòng),地表松散堆積或風(fēng)化殘坡積層上因局部拉張而發(fā)育伸展構(gòu)造,如張裂縫、地震凹陷和拉分盆地(圖7)。松散覆蓋層發(fā)育的拉分盆地從構(gòu)造上與黏土物質(zhì)中的雁列狀張裂縫類似,其規(guī)模與剪切位移呈正相關(guān)的關(guān)系[31-32]。演化模式如圖7(d),隨著變形的繼續(xù),雁列狀斜列的張裂縫發(fā)生疊覆,單個(gè)張裂縫或小型拉分盆地逐漸變寬、變深,最后聯(lián)合起來形成規(guī)模更大的盆地。圖7(c)中拉分盆地規(guī)模及形態(tài)發(fā)育程度高,可能對(duì)應(yīng)多次地震事件累計(jì)而成。

圖7 里德爾剪切盆地演化模式Fig.7 Evolution model of Riedel shear basin

4.2.2 地震裂縫與縮短構(gòu)造

在地震破裂過程中,受兩條左行右階破裂相對(duì)運(yùn)動(dòng),在地表松散堆積或風(fēng)化殘坡積層上因局部壓扭形成縮短構(gòu)造,如地震鼓包、鼓梁(圖8)。前者為寬度略小于長(zhǎng)度的小幅度隆起, 一般在平面上呈橢圓形狀,后者隆起長(zhǎng)度遠(yuǎn)大于寬度,平面上呈長(zhǎng)條狀。地震鼓梁長(zhǎng)軸方向一般在N45 °～ 70 °W,在地表破裂帶各段落均有不同程度發(fā)育。張(剪)性裂縫段與壓扭性鼓包段相間排列,是整個(gè)破裂帶分布最為普遍的破裂形態(tài)組合[33]。相鄰鼓包呈左階斜列展布, 其間以右階斜列展布的張(剪)性裂縫首尾相連。前述圖3(a)中張剪裂縫呈波折狀不連續(xù)延伸,其階區(qū)位置實(shí)質(zhì)上也為地震鼓包(梁)相連。

地震裂縫的空間分布特征及組合形式反映地表破裂的幾何學(xué)、運(yùn)動(dòng)學(xué)特征,揭示了構(gòu)造活動(dòng)下地震裂縫的動(dòng)態(tài)發(fā)育過程。需要指出,地震裂縫的空間分布對(duì)于其他伸展、縮短型地表破裂的圈定具有重要的指示意義,后者雖然形態(tài)不一,無法通過遙感手段直接提取,但可根據(jù)地震裂縫階區(qū)間接圈劃出其可能發(fā)育位置。

圖8 張剪裂縫與縮短構(gòu)造的組合形式Fig.8 Combination of shear fracture and compression textures

5 結(jié)論與展望

震源斷層錯(cuò)動(dòng)形成的地震裂縫與發(fā)震斷裂帶位置基本重合,除光譜特征外,空間分布也具有明顯的方向性特征。本文選取昆侖山口西地震同震地表破裂帶為研究區(qū),建立了一套基于光譜、空間特征的地震裂縫提取方法。對(duì)比目前常用地裂縫提取方法,結(jié)果表明該方法結(jié)果最優(yōu),MF-FDOG算法次之,改進(jìn)的邊緣檢測(cè)算法效果較差。該方法適用于重點(diǎn)區(qū)域地震裂縫的精細(xì)刻畫描述,本文選取試驗(yàn)區(qū)各類型地震裂縫形態(tài)完整,具備其典型的光譜特征和空間特征,在地表破裂提取方法研究和應(yīng)用中能起到一定的參考價(jià)值。如何在更為復(fù)雜的背景中提取地震裂縫,將是今后的研究重點(diǎn)。

運(yùn)用里德爾剪切理論在遙感影像中識(shí)別地震裂縫,分析走滑型地震同震地表破裂的運(yùn)動(dòng)特征與空間展布關(guān)系。將裂縫幾何形狀及延伸方向與里德爾剪切結(jié)構(gòu)相關(guān)聯(lián),能快速確定主剪切帶方向和區(qū)域應(yīng)力場(chǎng),有助于進(jìn)行發(fā)震構(gòu)造動(dòng)力學(xué)、幾何學(xué)等相關(guān)研究。同時(shí)指出,地震裂縫的空間分布對(duì)于其他伸展、縮短型地表破裂的圈定具有重要的指示意義。

本研究引入國(guó)產(chǎn)高分衛(wèi)星GF-2作為影像數(shù)據(jù)源,顯示了國(guó)產(chǎn)高分辨率遙感數(shù)據(jù)在地表破裂考察等地震科學(xué)研究方面的應(yīng)用價(jià)值。提取結(jié)果證明,GF-2衛(wèi)星完全滿足大地震后地表破裂考察的需求,可為災(zāi)后應(yīng)急與災(zāi)情評(píng)估提供有效的空間信息依據(jù)。