同震地表破裂帶空間分布形態的自動無損測量

鄧德貝爾 劉小利 高天琪 樂子揚

1 中國地震局地震研究所，武漢市洪山側路40號，430071 2 防災科技學院，河北省三河市學院街465號，065201

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

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

為高效、精細地測量、統計和分析大量同震地表破裂數據，本文設計并實現了同震地表破裂空間分布形態的無損測量方法和基于Python平臺的自動測量工具，并以2021年瑪多MW7.4地震和1983年Idaho Lost River MW6.9地震數據進行性能測試，驗證該方法的可信度和運算效率。

1 同震地表破裂帶分布形態的無損測量方法

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

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

此外，當地表破裂分布在數百米甚至更寬的范圍時，由于現場觀測視野的局限性或者距離主破裂較遠處次級破裂規模較小、變形微弱而被忽視(P1點南側的C、D區域)，使得部分破裂未被納入測量范圍，進而導致破裂帶寬度被低估[12]。從統計意義上來說，應量取測量點位處包含所有地表破裂的最大垂向距離作為該點位地表破裂帶的寬度。以圖1為例，A區、C區代表不同規模的破裂跡線，將A區、C區及中間空白區域合并測量，得到該段破裂帶的總體寬度。

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

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

2 算法設計與實現

基于上述無損測量方法，設計并實現了相應的自動無損測量和變尺度統計算法，詳細流程見圖2。

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

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

2)求包絡形跡點。以主破裂斷層為基準線，設置搜索窗口(圖1中E區)，調節窗口步長(scale)，等間隔掃描主斷層正交方向兩側最遠的點。此時，需要根據主斷層的走向(式(1))確定其正交方向，并進行包絡形跡點的判定。假設A、B兩點為主斷層上測量點位兩側相鄰的點，這兩點間地理方位角即為該測量點處斷層的走向值(α)：

(1)

式中，y為緯度，x為經度。

3)計算寬度。以測量點位主斷層走向為依據進行正交方向的限定，主斷層兩側包絡形跡點到主斷層的垂向距離即為該點處主斷層兩盤地表破裂的寬度，兩盤寬度之和即為該點處地表破裂帶總體寬度。當地表破裂帶存在空區(部分段落沒有破裂)時，程序會給定寬度值為0。此外，程序還允許以震中為水平距離零點位置，僅需調整起始端點的距離參數值即可。

4)統計分析與平滑。基于無損測量值，可根據實際需要靈活選擇分析窗口尺寸(如原始數據分辨率或其任意倍數)，進行最值、中位數、平均數等統計分析或變尺度平滑處理。

基于上述步驟獲得的結果，可進一步展示破裂帶的可視化空間分布形態。基于Paython的shapefile、pandas、scipy、matplotlib等工具包實現該算法。

3 算法測試與震例分析

3.1 實驗數據

為測試算法的可靠性和運算效率，分別選取1983年Idaho Lost River MW6.9地震[14]和2021年青海瑪多MW7.4地震的同震地表破裂數據進行實驗。

3.2 可靠性分析

3.2.1 Idaho Lost River地震

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

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

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

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

3.2.2 瑪多震例

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

圖5 分支破裂Fig.5 Branch rupture

圖6 破裂分布形態Fig.6 Rupture distribution pattern

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

3.3 運算效率測試

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

圖7 耗時結果Fig.7 Time consuming results

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

圖8 不同大小測量窗口對應的地表破裂空間分布形態Fig.8 Spatial distribution pattern of surface rupture corresponding to different measurement windows

較高的數據分辨率可緩解測量窗口選取對測量結果的干擾，但隨著數據分辨率的增加，數據量也越大，處理耗時越高。破裂趨勢線是基于地表破裂的整體趨勢確定的，實際上其走向變化頻率遠小于破裂線文件的最高分辨率精度。此外，根據破裂趨勢線的方向確定各個測量點位處正交方向的過程耗時最高。因此，基于相同的數據集，可根據實際需要設置不同的測量窗口因子，以窗口內破裂寬度的某一統計數值表示該區域破裂寬度，近似獲取整條破裂帶寬度的測量結果。總之，對于精細的地表破裂數據，無損數據處理必然導致時間消耗增大，為了在保留地表破裂精細特征的同時降低耗時，可以適當增大計算窗口尺寸或選擇并行處理。

4 討 論

4.1 測量結果的細分

利用本文算法可測量破裂帶的整體空間分布、破裂分支的空間分布及破裂空區，但對于規模較小的次級破裂或破裂分布較為彌散的情況，無法自動細分處理，需要人為干預。

4.2 運算性能的優化

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

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

4.3 最佳窗口步長的選取

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

5 結 語

本文基于Python平臺，針對高密集度的同震地表破裂數據，設計并實現了一種自動、無損的地表破裂空間形態測量方法和變尺度統計、平滑、制圖工具，可實現自動、高效、精細地測量及繪制地表破裂的空間分布形態，提高地震地質調查數據處理與分析效率和精細程度，可為活動斷層避讓帶設定、重點基礎設施抗震設防與風險防范提供更加準確的參考。