同震地表破裂帶空間分布形態(tài)的自動(dòng)無損測(cè)量

鄧德貝爾 劉小利 高天琪 樂子揚(yáng)

1 中國(guó)地震局地震研究所，武漢市洪山側(cè)路40號(hào)，430071 2 防災(zāi)科技學(xué)院，河北省三河市學(xué)院街465號(hào)，065201

地震地表破裂帶包含了重要的幾何結(jié)構(gòu)關(guān)系，其幾何形態(tài)、同震變形及其在空間上的分布特征是理解大陸地殼變形模式和地震破裂行為的關(guān)鍵[1-5]，也是構(gòu)造活動(dòng)區(qū)各種重要基礎(chǔ)設(shè)施抗震設(shè)防的依據(jù)[6-8]。已有學(xué)者根據(jù)歷史震例統(tǒng)計(jì)規(guī)律建立了同震地表位錯(cuò)與震級(jí)大小之間的函數(shù)映射關(guān)系，稱之為地震危險(xiǎn)性概率模型[6-7,9]，用于活動(dòng)斷層避讓帶的限定。但目前對(duì)同震地表破裂帶寬度及其影響因素缺少足夠的關(guān)注，不僅影響對(duì)地震破裂行為的全面認(rèn)識(shí)，還制約了活動(dòng)斷層工程避讓帶寬度及其邊界的準(zhǔn)確限定[10]。

活動(dòng)斷層避讓帶寬度的測(cè)定方法主要有2種：跨斷層地質(zhì)探槽剖面分析法和同震地表破裂帶寬度統(tǒng)計(jì)法[11]。前者需要實(shí)地開挖探槽，價(jià)格昂貴，一般適用于工地建設(shè)的精細(xì)探測(cè)；后者則是在震后應(yīng)急調(diào)查中通過實(shí)地抽樣調(diào)查和測(cè)量，獲得主破裂斷層沿線若干點(diǎn)位處破裂帶寬度測(cè)量值，經(jīng)過空間插值擬合得到整條破裂帶的寬度分布。隨著亞m級(jí)甚至cm級(jí)高空間分辨遙感技術(shù)在地震地質(zhì)調(diào)查中的應(yīng)用，快速獲得大空間范圍、詳盡的、精細(xì)的同震地表破裂信息成為可能，以往地面調(diào)查難以覆蓋到的地表破裂或現(xiàn)場(chǎng)肉眼易忽略的微弱地表裂縫可被再現(xiàn)[12]，通過室內(nèi)精細(xì)解譯和填圖即可獲得空間上幾乎連續(xù)的同震地表破裂的分布形態(tài)。

為高效、精細(xì)地測(cè)量、統(tǒng)計(jì)和分析大量同震地表破裂數(shù)據(jù)，本文設(shè)計(jì)并實(shí)現(xiàn)了同震地表破裂空間分布形態(tài)的無損測(cè)量方法和基于Python平臺(tái)的自動(dòng)測(cè)量工具，并以2021年瑪多MW7.4地震和1983年Idaho Lost River MW6.9地震數(shù)據(jù)進(jìn)行性能測(cè)試，驗(yàn)證該方法的可信度和運(yùn)算效率。

1 同震地表破裂帶分布形態(tài)的無損測(cè)量方法

已有研究認(rèn)為，地震地表破裂帶一般由一系列性質(zhì)不同的次級(jí)斜列斷層、鼓包、裂縫、溝槽等組合而成[1,13]。目前，地表破裂帶寬度主要指集中在發(fā)震斷層上的同震地表破裂在斷層正交方向上的延展寬度[13]，根據(jù)測(cè)量對(duì)象的空間尺度不同，又可分為單條地表破裂寬度和地表破裂帶組合寬度。

準(zhǔn)確限定同震地表破裂帶的幾何學(xué)寬度并不容易。在野外現(xiàn)場(chǎng)調(diào)查中，一般沿主破裂斷層走向每間隔一定距離抽樣測(cè)量一個(gè)破裂寬度值作為該段落地表破裂帶的寬度。如圖1，該段破裂帶總體上沿aa′展布，走向大致為NEE，長(zhǎng)約1 km，由多支近平行或小角度斜交的斜列式張剪性破裂組成。在實(shí)地調(diào)查中，首先根據(jù)破裂的延展性、規(guī)模大小和總體走向，確定主破裂位置作為基準(zhǔn)，測(cè)量其正交方向bb′兩側(cè)破裂的最遠(yuǎn)位置，作為測(cè)量點(diǎn)處地表破裂帶的寬度值。例如，在P1點(diǎn)的測(cè)量值為48 m，作為該點(diǎn)兩側(cè)一定范圍(P1點(diǎn)處陰影區(qū))的破裂帶寬度值，但該點(diǎn)東南側(cè)A區(qū)仍有地表破裂未被覆蓋；P2點(diǎn)的測(cè)量值為166 m，但該點(diǎn)東南側(cè)B區(qū)較大范圍內(nèi)沒有地表破裂?？梢?，在不同位置上抽樣測(cè)量所獲得的測(cè)量值往往會(huì)偏大或偏小，具有較大的不確定性。

此外，當(dāng)?shù)乇砥屏逊植荚跀?shù)百米甚至更寬的范圍時(shí)，由于現(xiàn)場(chǎng)觀測(cè)視野的局限性或者距離主破裂較遠(yuǎn)處次級(jí)破裂規(guī)模較小、變形微弱而被忽視(P1點(diǎn)南側(cè)的C、D區(qū)域)，使得部分破裂未被納入測(cè)量范圍，進(jìn)而導(dǎo)致破裂帶寬度被低估[12]。從統(tǒng)計(jì)意義上來說，應(yīng)量取測(cè)量點(diǎn)位處包含所有地表破裂的最大垂向距離作為該點(diǎn)位地表破裂帶的寬度。以圖1為例，A區(qū)、C區(qū)代表不同規(guī)模的破裂跡線，將A區(qū)、C區(qū)及中間空白區(qū)域合并測(cè)量，得到該段破裂帶的總體寬度。

基于上述分析，為客觀、真實(shí)地再現(xiàn)同震地表破裂帶的空間分布形態(tài)，滿足不同研究的差異性需要，本文提出一種無損失逐點(diǎn)掃描、變尺度統(tǒng)計(jì)的測(cè)量方法，即以主破裂斷層為基準(zhǔn)線，由一端向另一端等間隔移動(dòng)掃描與主破裂斷層正交方向上兩側(cè)破裂的最大垂向距離，作為該測(cè)量點(diǎn)位破裂帶的最遠(yuǎn)包絡(luò)形跡點(diǎn)，所有測(cè)量點(diǎn)位處的最遠(yuǎn)包絡(luò)形跡點(diǎn)即可形成地表破裂帶的包絡(luò)多邊形，測(cè)量點(diǎn)位到主破裂斷層起始端點(diǎn)的長(zhǎng)度即為測(cè)量點(diǎn)到主斷層起始端點(diǎn)的實(shí)際距離。因此，地表破裂空間分布形態(tài)的測(cè)量過程，就是繪制地表破裂包絡(luò)多邊形的過程，也是將空間坐標(biāo)系統(tǒng)轉(zhuǎn)換到以主破裂斷層走向方向及其正交方向?yàn)榛鶞?zhǔn)坐標(biāo)軸的投影轉(zhuǎn)換過程。

顯然，該包絡(luò)多邊形的精度取決于原始數(shù)據(jù)的分辨精度和測(cè)量密度，因此沿主破裂斷層掃描時(shí)，必須以最高分辨率設(shè)置移動(dòng)窗口尺寸才能實(shí)現(xiàn)無損測(cè)量。在此基礎(chǔ)上，取移動(dòng)窗口尺寸的任意倍數(shù)對(duì)數(shù)據(jù)進(jìn)行最值、中位數(shù)、平均值等統(tǒng)計(jì)分析或制圖，滿足變尺度研究需要。

2 算法設(shè)計(jì)與實(shí)現(xiàn)

基于上述無損測(cè)量方法，設(shè)計(jì)并實(shí)現(xiàn)了相應(yīng)的自動(dòng)無損測(cè)量和變尺度統(tǒng)計(jì)算法，詳細(xì)流程見圖2。

圖2 算法流程Fig.2 Algorithm flow chart

1)數(shù)據(jù)重采樣。將地表破裂線、主破裂斷層線文件按最高分辨率大小離散為點(diǎn)文件，以實(shí)現(xiàn)無損轉(zhuǎn)換與測(cè)量。如基于5 cm像素分辨率的航空影像獲得同震地表破裂信息，其最小分辨能力即為5 cm，將破裂線離散為間隔5 cm的一組點(diǎn)，可無損地再現(xiàn)原始數(shù)據(jù)包含的破裂信息。線文件可以是地理投影或平面投影，算法自動(dòng)檢測(cè)并轉(zhuǎn)換為平面投影坐標(biāo)。

2)求包絡(luò)形跡點(diǎn)。以主破裂斷層為基準(zhǔn)線，設(shè)置搜索窗口(圖1中E區(qū))，調(diào)節(jié)窗口步長(zhǎng)(scale)，等間隔掃描主斷層正交方向兩側(cè)最遠(yuǎn)的點(diǎn)。此時(shí)，需要根據(jù)主斷層的走向(式(1))確定其正交方向，并進(jìn)行包絡(luò)形跡點(diǎn)的判定。假設(shè)A、B兩點(diǎn)為主斷層上測(cè)量點(diǎn)位兩側(cè)相鄰的點(diǎn)，這兩點(diǎn)間地理方位角即為該測(cè)量點(diǎn)處斷層的走向值(α)：

(1)

式中，y為緯度，x為經(jīng)度。

3)計(jì)算寬度。以測(cè)量點(diǎn)位主斷層走向?yàn)橐罁?jù)進(jìn)行正交方向的限定，主斷層兩側(cè)包絡(luò)形跡點(diǎn)到主斷層的垂向距離即為該點(diǎn)處主斷層兩盤地表破裂的寬度，兩盤寬度之和即為該點(diǎn)處地表破裂帶總體寬度。當(dāng)?shù)乇砥屏褞Т嬖诳諈^(qū)(部分段落沒有破裂)時(shí)，程序會(huì)給定寬度值為0。此外，程序還允許以震中為水平距離零點(diǎn)位置，僅需調(diào)整起始端點(diǎn)的距離參數(shù)值即可。

4)統(tǒng)計(jì)分析與平滑?；跓o損測(cè)量值，可根據(jù)實(shí)際需要靈活選擇分析窗口尺寸(如原始數(shù)據(jù)分辨率或其任意倍數(shù))，進(jìn)行最值、中位數(shù)、平均數(shù)等統(tǒng)計(jì)分析或變尺度平滑處理。

基于上述步驟獲得的結(jié)果，可進(jìn)一步展示破裂帶的可視化空間分布形態(tài)?；赑aython的shapefile、pandas、scipy、matplotlib等工具包實(shí)現(xiàn)該算法。

3 算法測(cè)試與震例分析

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

為測(cè)試算法的可靠性和運(yùn)算效率，分別選取1983年Idaho Lost River MW6.9地震[14]和2021年青?，敹郙W7.4地震的同震地表破裂數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。

3.2 可靠性分析

3.2.1 Idaho Lost River地震

如圖3所示，基于文獻(xiàn)[15]提供的1983年Idaho Lost River MW6.9地震空間最小單元為1 m的地表破裂數(shù)據(jù)，測(cè)量窗口分別設(shè)置為1 km、100 m、1 m，獲得的地表破裂帶寬度變化特征見圖4。

圖3 文獻(xiàn)[15]選用的斷層及寬度(WGS 1984 UTM Zone 12N坐標(biāo)系)Fig.3 Fault and width used in the literature[15](WGS 1984 UTM Zone 12N coordinates)

圖4顯示，測(cè)量窗口為100 m或1 m時(shí)，測(cè)量結(jié)果幾乎一致，但前者耗時(shí)更少；當(dāng)測(cè)量窗口大于1 km時(shí)，測(cè)量結(jié)果失真較為嚴(yán)重。圖3中，8～12 km區(qū)間，地表破裂帶寬度達(dá)到第1個(gè)峰值，約為6～7 km；17～18 km區(qū)間，地表破裂帶寬度達(dá)到第2個(gè)峰值，約為5～6 km；27～34 km區(qū)間，地表破裂帶集中分布于3～5 km范圍內(nèi)；13～17 km區(qū)間為一段空區(qū)，表示該區(qū)間沒有地表破裂。圖4對(duì)應(yīng)位置測(cè)量值的整體趨勢(shì)與圖3結(jié)果一致，展現(xiàn)了地表破裂沿主斷層走向變化的更多細(xì)節(jié)和總體寬度變化。值得注意的是，在8～12 km區(qū)間出現(xiàn)的最大峰值并不代表在6～7 km范圍內(nèi)廣泛分布著同震地表破裂。結(jié)合地表破裂的空間分布(圖3)，此處存在兩支斜交的同震地表破裂帶，使破裂帶整體寬度較大。

圖4 地表破裂空間分布形態(tài)Fig.4 Spatial distribution pattern of surface rupture

3.2.2 瑪多震例

為進(jìn)一步考察數(shù)據(jù)分辨率對(duì)測(cè)量結(jié)果的影響，選擇2021年瑪多MW7.4地震空間最小單元為3 cm的地表破裂數(shù)據(jù)，并根據(jù)光學(xué)遙感影像解譯和D-InSAR反演形變場(chǎng)指示的發(fā)震斷層[12]勾繪主破裂斷層趨勢(shì)線進(jìn)行測(cè)試。如圖5所示，該段地表破裂由北(AB)、南(CD)兩支破裂帶組成。為無損失地保留原始數(shù)據(jù)的精度，線文件離散為點(diǎn)文件時(shí)，以原始數(shù)據(jù)分辨率(3 cm)作為移動(dòng)窗口的尺寸。如前文所述，趨勢(shì)線的選取(大致沿破裂帶整體走向)對(duì)破裂帶整體分布寬度的量取影響較小，但會(huì)影響兩側(cè)破裂帶分布寬度的測(cè)量結(jié)果。如果將圖5中所有地表破裂基于AB或CD任意一條趨勢(shì)線進(jìn)行處理，其結(jié)果都無法準(zhǔn)確表征破裂帶的細(xì)節(jié)特征。因此，當(dāng)破裂分叉時(shí)，需要對(duì)破裂分支分別進(jìn)行趨勢(shì)線擬合，分支破裂分別擬合時(shí)的分布形態(tài)結(jié)果見圖6(a)和6(b)，整條破裂帶的分布形態(tài)見圖6(c)。

圖5 分支破裂Fig.5 Branch rupture

圖6 破裂分布形態(tài)Fig.6 Rupture distribution pattern

圖6(a)中，AB段破裂西側(cè)的寬度明顯大于東側(cè)；以趨勢(shì)線為參考，其北側(cè)地表破裂寬度明顯小于南側(cè)，很好地展示了圖5中AB、CD段破裂的空間分布形態(tài)。圖5中地表破裂僅基于AB或CD任意一條趨勢(shì)線處理，會(huì)對(duì)C～B區(qū)間大片空白區(qū)域進(jìn)行計(jì)算，造成較大的測(cè)量誤差。顯然，該震例中同震地表破裂空間分布形態(tài)的無損測(cè)量是基于原始數(shù)據(jù)的分辨率來完成的，數(shù)據(jù)分辨率越高，測(cè)量結(jié)果越精細(xì)，可反映的破裂細(xì)節(jié)特征越多。

3.3 運(yùn)算效率測(cè)試

為測(cè)試算法的運(yùn)算性能，選取瑪多地震一段長(zhǎng)約1 km的地表破裂進(jìn)行實(shí)驗(yàn)(圖1)，其空間分辨率為3 cm，密集分布了約3 900條破裂線，累積長(zhǎng)度約為5.95 km，離散后點(diǎn)集數(shù)達(dá)1.98×105個(gè)，按不同測(cè)量窗口(3 cm、10 cm、50 cm、1 m、5 m、10 m)進(jìn)行線轉(zhuǎn)點(diǎn)離散化、寬度計(jì)算、平滑及繪圖測(cè)試，耗時(shí)結(jié)果見圖7。可以看出，算法耗時(shí)主要取決于寬度計(jì)算。隨著數(shù)據(jù)分辨率的降低，時(shí)間消耗大致呈下降趨勢(shì)，但因數(shù)據(jù)分辨率不同而引起的寬度耗時(shí)差別不明顯。在實(shí)際應(yīng)用中，可根據(jù)對(duì)數(shù)據(jù)精細(xì)程度的要求在不同階段選擇合適的處理尺度。

圖7 耗時(shí)結(jié)果Fig.7 Time consuming results

圖8為利用不同大小測(cè)量窗口得到的地表破裂空間分布形態(tài)。可以看出，隨著測(cè)量窗口尺度的增大，地表破裂空間形態(tài)變得愈加平滑，部分細(xì)節(jié)特征丟失，破裂帶寬度的離散程度逐漸減小，尤其是當(dāng)測(cè)量窗口大于5 m時(shí)，破裂帶離散程度趨于穩(wěn)定；當(dāng)測(cè)量窗口在5～50 m時(shí)，測(cè)量寬度變化顯著，特征仍有保留；當(dāng)測(cè)量窗口超過150 m時(shí)，測(cè)量特征大量丟失。

圖8 不同大小測(cè)量窗口對(duì)應(yīng)的地表破裂空間分布形態(tài)Fig.8 Spatial distribution pattern of surface rupture corresponding to different measurement windows

較高的數(shù)據(jù)分辨率可緩解測(cè)量窗口選取對(duì)測(cè)量結(jié)果的干擾，但隨著數(shù)據(jù)分辨率的增加，數(shù)據(jù)量也越大，處理耗時(shí)越高。破裂趨勢(shì)線是基于地表破裂的整體趨勢(shì)確定的，實(shí)際上其走向變化頻率遠(yuǎn)小于破裂線文件的最高分辨率精度。此外，根據(jù)破裂趨勢(shì)線的方向確定各個(gè)測(cè)量點(diǎn)位處正交方向的過程耗時(shí)最高。因此，基于相同的數(shù)據(jù)集，可根據(jù)實(shí)際需要設(shè)置不同的測(cè)量窗口因子，以窗口內(nèi)破裂寬度的某一統(tǒng)計(jì)數(shù)值表示該區(qū)域破裂寬度，近似獲取整條破裂帶寬度的測(cè)量結(jié)果?？傊瑢?duì)于精細(xì)的地表破裂數(shù)據(jù)，無損數(shù)據(jù)處理必然導(dǎo)致時(shí)間消耗增大，為了在保留地表破裂精細(xì)特征的同時(shí)降低耗時(shí)，可以適當(dāng)增大計(jì)算窗口尺寸或選擇并行處理。

4 討 論

4.1 測(cè)量結(jié)果的細(xì)分

利用本文算法可測(cè)量破裂帶的整體空間分布、破裂分支的空間分布及破裂空區(qū)，但對(duì)于規(guī)模較小的次級(jí)破裂或破裂分布較為彌散的情況，無法自動(dòng)細(xì)分處理，需要人為干預(yù)。

4.2 運(yùn)算性能的優(yōu)化

在Idaho Lost River地震中，根據(jù)破裂帶走向分段批處理測(cè)量破裂分布形態(tài)，同時(shí)基于破裂分布密度剔除空白階區(qū)內(nèi)稀少破裂的誤差和粗差干擾，從而可獲得分布更集中的破裂跡線。這些手段有效減少了算法單次計(jì)算量，提高了大數(shù)據(jù)量文件的快速處理，但海量地表破裂數(shù)據(jù)的自動(dòng)化批處理策略還有待優(yōu)化。此外，實(shí)驗(yàn)發(fā)現(xiàn)設(shè)置較小的測(cè)量窗口因子有利于實(shí)現(xiàn)無損或相對(duì)更高精度的測(cè)量，但耗時(shí)過多。實(shí)際測(cè)量中宜采取“兩步法”的操作，即先用較大測(cè)量窗口因子粗測(cè)并劃分破裂密集區(qū)域，再用較小測(cè)量窗口因子實(shí)際測(cè)量，可有效兼顧無損測(cè)量與效率的平衡。但測(cè)量窗口因子的選取嚴(yán)重依賴經(jīng)驗(yàn)判斷，對(duì)精度與效率影響較大，測(cè)量窗口因子的自適應(yīng)選取將在后續(xù)研究中優(yōu)化。

瑪多地震展示了基于重采樣后的破裂數(shù)據(jù)密度分布跟蹤高密度區(qū)作為破裂趨勢(shì)線的形跡，結(jié)合密度閾值參數(shù)，有助于自動(dòng)判定規(guī)模較大的分支破裂，但密度閾值的確定需要人為干預(yù)，因此分支破裂的處理仍需人工二次判讀。

4.3 最佳窗口步長(zhǎng)的選取

Idaho Lost River地震和瑪多地震均顯示，隨著測(cè)量窗口的增大，破裂帶寬度的離散程度逐漸減小。因此，對(duì)于10 km長(zhǎng)的破裂帶，測(cè)量窗口初始值取50～100 m為宜；對(duì)于20 km以上長(zhǎng)度的破裂帶，測(cè)量窗口初始值取100～200 m為宜。原始數(shù)據(jù)分辨率對(duì)測(cè)量窗口取值有一定影響，但無直接關(guān)系，測(cè)量窗口的最佳取值與破裂帶變化分布特征關(guān)聯(lián)度更高。

5 結(jié) 語

本文基于Python平臺(tái)，針對(duì)高密集度的同震地表破裂數(shù)據(jù)，設(shè)計(jì)并實(shí)現(xiàn)了一種自動(dòng)、無損的地表破裂空間形態(tài)測(cè)量方法和變尺度統(tǒng)計(jì)、平滑、制圖工具，可實(shí)現(xiàn)自動(dòng)、高效、精細(xì)地測(cè)量及繪制地表破裂的空間分布形態(tài)，提高地震地質(zhì)調(diào)查數(shù)據(jù)處理與分析效率和精細(xì)程度，可為活動(dòng)斷層避讓帶設(shè)定、重點(diǎn)基礎(chǔ)設(shè)施抗震設(shè)防與風(fēng)險(xiǎn)防范提供更加準(zhǔn)確的參考。