瀘定Ms6.8地震房屋損失快速評估

趙登科,王自法,,李兆焱,周陽,高曹珀,WANG Jianming,位棟梁,張昕

(1.中國地震局工程力學研究所 中國地震局地震工程與工程振動重點實驗室,黑龍江 哈爾濱 150080;2.地震災害防治應急管理部重點實驗室,黑龍江 哈爾濱 150080;3.中震科建(廣東)防災減災研究院,廣東 韶關 512000;4.河南大學 土木建筑學院,河南 開封 475004)

0 引言

地震是地殼快速釋放能量過程中引起的地表劇烈振動,通常會引起大量的人員傷亡和經濟損失[1],地震發生時大多數損失是由房屋的倒塌與破壞引起的,人員傷亡與商業中斷也與房屋的破壞程度有極大的關系。因此利用科學的方法快速準確地評估震后房屋的損失后果以及可能的分布,對于震后指揮、緊急救援、分配救援人員和物資有著重要意義[2]。

目前常用的震后損失評估方法主要有：(a)基于歷史震害資料的震害矩陣方法[3-6]。(b)基于性能的易損性分析方法[7-8]。(c)以現場調查為依據的抽樣統計方法。(d)基于遙感影像和航片判讀的識別法[9-11]。文獻[5]詳細闡述了上述幾種評估方法的流程,并指出后3種方法難以應用于破壞性地震的建筑群體損失快速估計中。建筑物群體快速損失評估目前應用最多的仍為震害矩陣方法[12],該方法的主要是依據震后發布的地震烈度圖,承災體資料以及所建立的房屋震害矩陣,分別計算各烈度區的房屋破壞比例,最終定量給出影響場內房屋的整體損失結果[13]。例如徐國棟等[5]對汶川地震損失估計;李成帥等[14]對蘆山地震的損失估計等。

傳統的震害矩陣方法雖然能夠快速地提供定量的震害損失評估,但是以烈度為基礎的評價方法精度有限,而且評估結果僅代表了損失的平均水平,忽略了地震損失中各類不確定性的影響[15]。由于震源參數、傳播路徑、場地效應和結構反應的多重復雜性,房屋的地震損失有著極大的不確定性[16-17],其對于量化尾部風險尤為重要。PAGER是美國聯邦地質調查局開發的全球地震損失快速估計系統[18],對于不同國家或地區,PAGER都定義了歸一化的標準差常數ζ(我國為1.895),用以描述總體損失的不確定性[19]。需要說明的是：該方法雖然一定程度上能夠刻畫整體損失的變異性,但同時意味著,對于不同震源位置和不同震級的地震,PAGER所估計的損失分布形狀都是固定的,這顯然過于簡化和失準,最終導致計算結果及其分布的誤差。事實上,由于各類不確定性的相互傳遞,地震影響場內每個空間位置點的損失都是隨機的[20],需要均值、方差和具體的分布形狀來共同確定。另外,在城市尺度或者區域層次地震估計損失中,各空間位置點的損失相關關系也是影響地震損失結果概率分布的重要因素[21-23]。因此,震后房屋損失評估是一個涉及多維相關變量的隨機過程,僅用均值或者固定的分布形狀無法描述房屋地震損失的客觀特征,若要得到接近實際的高精度損失結果,必須進行大量的樣本計算和模擬分析。

本文基于Copula理論,提出了一種適用于地震巨災風險分析的相關隨機變量模擬方法,能夠更準確地刻畫地震損失評估過程中的一系列不確定性和損失相關性,并利用所提方法對2022年9月5日瀘定6.8級地震的損失進行快速評估,給出了損失估計結果及其概率分布。希望研究方法和結果為震后損失快速評估技術提供參考,也為此次地震的應急管理提供了依據。

1 地質背景

2022年9月5日,四川省甘孜州瀘定縣發生了6.8級地震,震源深度16 km,震源機制解為走滑型破裂,震中位于北緯29.59°,東經102.08°附近的磨西鎮海螺溝公園內。磨西鎮整體位于四川西部地區的鮮水河斷裂帶上[24],鮮水河斷裂帶是四川境內最長也是地震活動性最強的一條斷裂帶。記載表明[24]：自1700年以來發生7級以上地震高達9次,平均每百年發生3次大地震,此次地震是鮮水河斷裂上40年來發生的最大地震,震中附近活動斷層分布如圖1所示。鮮水河斷裂帶西起甘孜東谷北,向東南延伸,經爐霍、道孚和康定,南達石棉,呈北西-南東走向,長約350 km,其位于巴顏喀拉塊體與川滇塊體的走滑活動邊界,與龍門山斷裂帶和安寧河斷裂帶交匯構成了川西地區著名的“Y”字形斷裂帶[25]。2008年汶川地震和2013年蘆山地震使得鮮水河斷裂附近庫侖力明顯增加,此次地震釋放了康定-石棉段積累的應變能[26]。

圖1 瀘定地震震中附近活動斷層分布 圖2 各單元網格房屋價值的空間分布Fig. 1 Distribution of active faults near the epicenter of the Luding earthquake Fig. 2 Spatial distribution of the building values of each cell grid

2 數據準備

2.1 建筑庫存和單元網格劃分

全國房屋普查統計數據顯示,震中周圍縣區的建筑類型以砌體結構和框架結構為主,其他結構類型占比極少,為實現損失的快速估算,本文僅考慮砌體和框架結構這2種主要的房建筑類型。在震后損失評估中,通常需要對風險評估區進行地理信息編碼[27],例如單元網格化處理,良好的地理信息編碼能夠在不影響準確度的前提下,快速地確定每個地理位置可能產生的損失大小。該研究通過網格化處理方法,將各縣(區)級別的行政區劃建筑財產數據轉換為多等級非均勻的網格數據。網格尺寸按精度降序分別為：10 km×10 km、5 km×5 km、1 km×1 km、500 m×500 m和250 m×250 m,行政區劃的人口密度、GDP以及占地面積等因素越高(大),網格劃分精度也就越高。依據轉換得到的網格數據,統計得出各單元網格的房屋價值空間分布,如圖2所示,發現自西向東各縣區網格精度呈階梯狀升高,符合實際情況。

2.2 地震動強度

震后第一時間,研究團隊基于瀘定地震的震中位置、震級和震源機制解等信息,并參考余言祥等[28]的地震動衰減關系見式(1),對此次地震事件的地震動強度分布情況進行模擬。

lgy=A+BM+Clg(R+DeEM)+ε

(1)

式中：y為地面峰值加速度PGA,M為震級,R為震中距,A、B、C、D和E為回歸系數,ε為標準差,回歸系數的系數值見表1。

表1 所用衰減模型的回歸系數[28]Table 1 Regression coefficients of the ground motion prediction model in use

由于缺少可靠的臺站場地信息,因此將臺站觀測記錄按照震源距進行均值化處理,將模擬得到的結果與臺站記錄均值對比,如圖3所示,發現模擬結果與實際觀測結果較為接近,說明本文的地震動強度模擬結果是合理的。此外,為了與應急管理部發布的烈度圖對比,參考《中國地震烈度表》(GB/T 17742—2020)[29],將本文模擬得到的PGA換算為地震烈度,得到地震烈度分布如圖4所示,震中地區的地震烈度達IX度以上,整體災害空分布大致呈橢圓形,可能受災范圍約18 115 km2。其中：Ⅸ度區面積為168 km2,主要涉及瀘定縣;Ⅷ度區面積為1 320 km2,主要涉及瀘定縣及石棉縣;Ⅶ度區面積為4 463 km2,主要涉及康定市、瀘定縣、九龍縣、漢源縣、石棉縣;Ⅵ度區面積為12 164 km2,主要涉及康定市、九龍縣、石棉縣、冕寧縣、甘洛縣、漢源縣、滎經縣和天全縣。根據結果可知：本文模擬得到的烈度分布以及各烈度區的面積與應急管理部發布的結果基本一致[30]。

圖3 衰減結果對比臺站記錄 圖4 本文計算得到的地震烈度分布Fig. 3 Comparison of ground motion prediction model results with actual station records Fig. 4 Earthquake intensity distribution obtained in this study

2.3 結構易損性

結構易損性描述了地震動與結構破壞之間的關系,相比于震害矩陣,脆弱性曲線,易損性曲線將地震動和結構破壞參數均連續化,因此最為準確[31]。因為缺少大量翔實的震害資料,基于統計方法得到的易損性曲線較少,目前最常見的是基于分析方法的易損性曲線。圖5中連續的曲線為文獻[32-33]中的模型易損性曲線,散點為汶川地震統計得到的實際破壞數據,可以發現：在地震動強度較高時,無論是框架結構還是砌體結構,分析方法得到的易損性曲線明顯高于實際調查結果,這與日本311地震中的發現是一致的[34]。

圖5 模型易損性曲線與汶川地震實際調查的對比Fig. 5 Comparison of model vulnerability curves with the actual survey in the Wenchuan earthquake

顯然,基于分析方法的易損性曲線不能直接應用于地震損失估計中,一個簡單有效的修正方法是：針對有限的震害調查點數據,計算各地震動強度下模型損失與實際損失的修正比例系數η,其余地震動強度處的修正系數通過線性內插計算。通過圖5中模型易損性與實際調查點的損失比例,得到各地震動強度下的修正系數η以及修正后的易損性曲線,如圖6-7所示。可以發現修正后的易損性曲線與實際的調查結果比較接近,能夠客觀地反映建筑物的抗震能力。需要說明的是：修正后得到的易損性結果符合仝文博等[35]和WESSON等[36]發現的“震害飽和”現象,即隨著地震動的增大,結構的損失不會明顯上升。

圖6 修正后的框架結構易損性曲線 圖7 修正后的砌體結構易損性曲線Fig. 6 Corrected vulnerability curve for RC structureFig. 7 Corrected vulnerability curve for masonry structure

3 房屋地震損失的模擬方法

3.1 不確定性與空間相關性

依據所得到的地震動與結構易損性,能夠獲得每個結構的平均損失程度。但實際震害調查發現：即使是相同地震動水平下,結構損失仍包含著高度的不確定性,因此僅利用平均損失程度無法描述房屋損失的客觀規律,需要進一步考慮結構地震損失的分布特征。在以往的研究中,我們針對新西蘭和日本大量震害資料,發現相同地震動水平下結構的損失近似服從Beta分布[37]。胡少卿等[38]基于寧蒗地震的實際調查數據,發現我國房屋的地震損失特征也服從Beta分布。由于Beta分布的非線性不可累加特征,在一個區域的損失計算中,無法直接利用Beta分布描述整體損失的不確定,需要首先對空間各單元網格點分別采樣確定損失后,再將具體數值累加得到整體損失的分布特征。

當各網格點的損失相互獨立時,可以借助蒙特卡洛模擬方法進行隨機模擬采樣。周陽等[21]基于實際的震害資料,發現各網格點的地震損失是存在相關性的,地震損失的相關性會對損失估計結果產生影響,忽略地震損失空間相關性會高估小震損失或低估大震損失。綜上,各網格點的損失就可以看作具有特定相關性結構的多維隨機變量X1,X2,… ,Xn,其中：n為空間位置的數量,對于任意的Xi,其邊緣分布都為Beta分布。若能夠構建出各網格點損失的聯合分布,則能夠容易地獲得到損失樣本,但多維隨機變量聯合分布的構造在相關理論推導和計算中都是極為繁瑣的,尤其是當隨機變量的數量比較多時,因此需要尋找更好的方法模擬損失。Copula理論提供了一種能夠處理隨機變量相關性問題的方法,此處引出Copula理論。

3.2 Copula理論

Copula理論基于Sklar定理[39]：假設H是n維隨機變量X1,X2,… ,Xn的聯合分布函數,與其對應的邊際分布分別是F1,F2,… ,Fn,那么就存在一個n元Copula函數C使得對于全部的x1,x2,… ,xn,有：

H(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn))=C(u1,u2,…,un)

(2)

若F1,F2,… ,Fn是連續的,則C唯一。因此有以下推理：若C是邊緣分布分別為F1,F2,… ,Fn的Copula函數,則式(2)定義的H(x1,x2,… ,xn)是隨機變量X1,X2,… ,Xn的聯合分布函數。Sklar定理給出了一種利用邊際分布對多元聯合分布建模的方法,由式(2)可知：多維隨機變量X1,X2,… ,Xn的概率密度函數為：

(3)

式(3)說明：一個多維的聯合概率密度函數可以拆解為Copula函數與邊緣密度概率函數乘積表達的形式,其核心思想是將邊緣分布和隨機變量間的相關結構分開進行研究,其中隨機變量的隨機性由各自的邊緣分布描述,隨機變量之間的耦合特性由Copula函數描述,因此Copula函數能夠用于構建具有相關性的隨機變量概率模型。

Copula函數總體上可以分為阿基米德型、橢圓型和二次型[40],由于橢圓型高斯Copula函數具有對稱性和計算簡單的特點,因此在實際中被廣泛應用,多維高斯Copula的分布函數如下：

(4)

式(4)中：Φ∑(·)為相關性矩陣滿足∑的高斯聯合概率分布函數,Φ(·)為標準的高斯分布函數,dx1,dx2,… ,dxn為積分變量。依據Copula理論,式(4)中N元隨機數向量u1,u2,… ,un保留了與所研究隨機變量相同的相關結構,即同為∑,上述過程很好地解決了隨機變量中相關性的問題。

3.3 地震損失采樣流程

由概率積分變換可知：對于任意一個已知的連續隨機變量X,其累計概率分布函數U=F(x),U都服從[0,1]上的均勻分布：

P(U≤y)=P(F(x)≤y)=P(F-1(F(x))≤F-1(y))
=P(x≤F-1(y))=F(F-1(y))=y

(5)

式中：y∈[0,1]。由式(5)可知：通過Copula函數得到的n維隨機數向量U1,U2,… ,UN是服從[0,1]上的均勻分布。由分位數轉換的通用性方法可知：只需要以U1,U2,… ,Un為累積分布概率值,對目標邊緣分布(我們所用為Beta)求逆累積概率運算,即可得到即既滿足特定的相關結構,又滿足各自的目標邊緣分布的多維隨機變量。

隨機問題通常需要借助模擬方法采樣得到具體的損失結果隨機數,本節以蒙特卡洛模擬方法為基礎,結合高斯Copula函數,建立一種適用于地震損失的采樣方法。在實際應用過程中,問題具體轉化為建立符合相關矩陣為R的N元隨機變量向量ZBeta,其中：R為房屋地震損失的經驗相關矩陣,可以利用地震損失相關性與空間距離的經驗公式建立,R∈[0,1]n×n,且滿足Rij=Rji,Rii=1,N為地震影響場內的空間位置數量,ZBeta為結構的損失比。采樣過程可以分解為以下3個步驟：

1)利用高斯Copula函數生成相關矩陣為R的N元隨機向量XGauss1,…,XGaussN,該步可以結合蒙特卡洛抽樣高斯分布獨立隨機數與Cholesky分解的方法實現[41]。

2)計算上步得到的各隨機向量XGauss1,…,XGaussN的高斯累積分布函數,將N元隨機數向量XGauss轉換為服從均勻分布[0,1]的YUniform1,…,YUniformN。

圖8展示了所述方法的采樣流程,圖8(c)中,每個單元網格點的損失邊緣分布ZBeta是根據該點的地震動強度與易損性曲線中的損失均值方差共同確定的,且會隨著Beta分布控制參數a和b的變化呈現出不同的形狀,滿足不同地震動強度下結構損失分布形狀改變的需要。

圖8 基于Copula理論的地震損失采樣流程Fig. 8 Seismic loss sampling process based on Copula theory

當多元變量的維數較低時,可以使用該方法直接抽取樣本。然而當維度很高時,尤其是估計破壞性地震的損失時,網格點的數量可以輕易地超過十萬級別,即使是對每個網格點進行一次模擬采樣,都需要消耗大量的算力,影響震后應急工作的開展。一個較好的方法是：先選擇少量位置生成隨機樣本,然后利用克里金插值方法得到其他位置的采樣值,WANG[42]已經證明了克里金插值對于巨災響應方法的有效性。

4 瀘定地震損失估計結果

4.1 采樣次數的影響

方法中的損失隨機數是基于蒙特卡洛模擬實現的,理論上,只有在足夠大的采樣規模下才能得到較高的精度。但樣本數增加會帶來計算量的大幅上升,影響震后應急工作的開展,模擬次數過少,則隨機數的分布不均勻,影響模擬結果的可靠性,因此采樣次數對結果的影響不可忽略。

本文通過式(6)衡量不同模擬次數下的相對誤差ε：

(6)

式中：μk為K次采樣后得到的損失均值,ELoss為損失期望基準值,可利用所建立網格數據和易損性曲線得到：

(7)

式(7)中：N為網格數量,M為結構分類數量,計算得到的損失期望(54.29億元人民幣)。針對此次地震,分別進行1～3 000次采樣,分析隨著抽樣次數增加時相對誤差ε的變化,以確定合適的采樣次數,結果如圖9所示。可以發現：當模擬次數達到2 000次時,損失結果趨于基準值,且誤差能夠控制在0.5%以內,綜合考慮計算的精度和速度,認為2 000次模擬結果滿足要求。

圖9 模擬次數對采樣精度的影響Fig. 9 Effect of simulation times on sampling accuracy

4.2 瀘定地震損失估計結果

基于上述流程,利用多維相關隨機變量的抽樣方法,對所建立的各網格點損失模擬2 000次得到了瀘定地震損失分析結果。圖10為框架結構和砌體結構的損失概率分布,可以發現雖然框架結構和砌體結構的損失最大可能在10～20億元人民幣,但是兩者在該損失區間的概率值相差較大,差值達到了23.8%。圖10給出了砌體結構和框架結構在不同損失區間下的概率差值,在損失小于20億元人民幣的區間內,框架結構發生損失的概率都遠大于砌體結構,在損失大于20億元人民幣的各區間內,框架結構發生損失的概率均小于等于砌體結構,意味著砌體結構較大概率會發生更高的損失。兩種結構類型損失的均值分別為35.82億元人民幣和18.63億元人民幣,砌體結構的損失均值幾乎是框架結構的2倍。原因可歸結于兩方面：其一,砌體結構抗震能力相比于框架結構較差;其二,地震影響場內,尤其是震中附近(鄉鎮)的砌體結構房屋占比較大。

圖10 不同結構類型的損失分布對比 圖11 不同縣區地震損失的均值和變異系數Fig. 10 Loss distribution comparison of different structure types Fig. 11 Mean and coefficient of variation of earthquake losses in different counties

依據單元網格,分別計算了各縣區的地震損失評估結果,由于受地震影響的縣區較多,此處僅展示損失較高的15個縣區的房屋損失均值和變異系數,如圖11所示。可以看到：瀘定縣房屋損失最大,其次是石棉縣和滎經縣,損失均值分別達到了17.81億元人民幣、12.07億元人民幣和6.4億元人民幣。從損失均值來看：損失的估計結果與圖4的烈度空間分布具有相似性,符合客觀規律,需要說明的是：雖然漢源縣的地震烈度比滎經縣更高,但圖2顯示滎經縣各單元網格的房屋價值高于漢源縣,因此損失估計結果顯示滎經縣的損失均值(6.4億元人民幣)略高于漢源縣(5.3億元人民幣)。此外我們發現：不同縣區的損失變異系數存在較大差異,意味著僅用一個恒定的標準差系數無法描述不同空間位置的損失的不確定性特征。圖12(a)和圖12(b)給出了損失較高的瀘定縣和石棉縣的損失概率分布,圖12(c)和圖12(d)給出了相同均值下,利用PAGER處理損失不確定性方法得到的結果。對比圖12(a)和圖12(b)和圖12(c)和圖12(d)可以發現：兩種方法損失分布結果有較大差異,由于不提前設置標準差常數,本文方法給出的損失概率分布較為靈活,能夠較好地反映不同縣區損失的特征。而PAGER對于不同縣區給出的損失分布形狀卻相對固定,無法反映出不同地區房屋空間分布的差異性。兩種方法的對比再次驗證了地震損失中相關隨機采樣方法的重要性。

圖12 本文方法與PAGER方法計算損失分布的對比Fig. 12 Comparison of earthquake building loss distribution between this study and PAGER

地震影響場內房屋總體損失的分布如圖13所示,損失的概率大致呈現出右偏態的長尾分布,損失最大的可能出現在20～30億元人民幣之間,但概率僅為18%。損失在20～50億元人民幣中的各區間的概率均大于15%。雖然分布中的峰值并不明顯,但從估計的結果中能夠看出：損失超過89.8%的概率處于10～100億元人民幣之間,能夠基本確定損失會處于該量級水平,且最有可能為20～50億元人民幣。

圖13 瀘定地震整體損失概率分布Fig. 13 Building loss probability distribution of Luding earthquake in this study

PAGER在震后發布了此次地震的損失估計結果,如圖14所示(單位為百萬美元),若按照美元對人民幣匯率為7.0計算,則計算結果顯示：此次地震損失在7～70億元人民幣和70-700億元人民幣的概率較大,分別為33%和31%,并且有13%的概率超過700億元人民幣。通過圖14可以計算出PAGER所估計的損失均值約為10億美元(70億元人民幣),該值在一定程度上與本文的估計結果較為近似,同時可以發現：PAGER所提供結果的離散型更大,從某種意義上來說缺乏精度。

圖14 PAGER提供的瀘定地震損失概率分布Fig. 14 Loss probability distribution of the Luding earthquake by PAGER

5 結論

針對地震損失中各網格位置點損失的隨機性和相關性,本文提出了考慮隨機變量相關性的隨機變量模擬方法,能夠用于未來的地震風險分析和震后損失快速估計中。利用所提方法評估了2022年9月瀘定6.8級地震的損失,得到了不同結構類型和不同縣區的損失概率分布結果,并與PAGER的估計結果進行了對比,得到了以下結論：

1)提出了一種適用于地震巨災風險分析的相關隨機模擬方法,并在單元網格,克里金插值方法和控制采樣次數3個方面進行數據的降維,大幅度地縮減了所需的計算量,能夠實現損失結果的快速估算。

2)瀘定地震損失估算結果顯示,相較于框架結構,砌體結構的損失更高,約是框架結構損失的2倍。此次地震中損失最高的幾個縣區分別為瀘定縣、石棉縣、滎經縣和漢源縣,其中前兩個縣的損失均值超過了10億元人民幣,分別為17.81億元人民幣和12.07億元人民幣。

3)利用PAGER計算損失不確定性的方法,分別得到了瀘定縣和石棉縣的損失概率分布,并與該方法得到的結果對比,發現PAGER對于不同縣區給出的損失分布形狀相對固定,無法反映不同地區損失特征的差異,驗證了所提方法的重要性。

4)估算結果顯示,瀘定地震造成的房屋損失在10～100億元人民幣量級水平的概率超過89%,并最有可能處于20～50億元人民幣之間。損失的均值為54.45億元人民幣,該結果在一定程度上接近于PAGER的估計結果。