基于廣義條件譜的某核電廠安全殼多元地震易損性分析

王曉磊，閻衛東，呂大剛

(1.沈陽建筑大學 土木工程學院，遼寧 沈陽 110168；2.河北省地震災害防御與風險評價重點實驗室，河北 三河 065201；3.哈爾濱工業大學 土木工程學院，黑龍江 哈爾濱 150090)

日本福島核事故后，核電廠地震安全評估受到全社會高度關注。抗震裕量評估[1]和概率地震風險評估[2]是核電廠地震安全評估的主要方法。由于概率地震風險評估方法具有分析結果更為精細化等優點，近年來得到了大量應用研究。概率地震風險評估主要由概率地震危險性分析、地震易損性分析、系統分析和損失分析等組成[2]，各組成部分精細化分析結果是精細化風險評估的基礎。

場地相關譜生成與地震動記錄選取是地震危險性分析與地震易損性分析的關鍵連接步驟，基于標量型概率地震危險性分析與分解的場地相關譜已得到廣泛應用，包括一致危險譜和條件均值譜等。一致危險譜的各譜加速度強度參數具有一致超越概率，譜型較為保守；條件均值譜考慮了譜型參數相關性信息，與實際地震動更為符合，但條件均值譜譜型較窄，作為場地相關譜輸入，可能偏于不保守。為了克服上述場地相關譜的不足，Kishida[3]提出了具有多個條件周期的廣義條件譜理論。目前，基于我國場地地震環境自身特點的廣義條件譜及地震動記錄選取研究還較為缺乏。

地震易損性分析是概率地震風險評估的重要步驟之一，對核電廠地震安全評估結果的精確性與合理性至關重要。安全系數法是核電廠地震易損性分析的常用方法，近年來，基于解析易損性數據的核電廠易損性分析得到了應用研究[4]。目前核電廠地震易損性通常基于單個地震動強度參數進行分析，單個地震動強度參數包含信息有限，得到的易損性分析結果存在較大不確定性。Cai等[5-6]提出了考慮多個地震強度參數的核電廠地震易損性分析方法和應用，但其提出的多參數地震易損性分析方法采用的是安全系數法，基于增量動力分析等解析地震易損性分析方法的核電廠多元易損性研究還較為缺乏。

本文考慮我國地震環境特點，提出適用于中國場地的廣義條件譜理論方法，生成我國某核電廠廠址廣義條件譜，采用基于精細化增量動力分析的多元易損性分析方法，進行我國某核電廠安全殼多元易損性分析計算，生成我國某核電廠安全殼地震易損性曲面，為未來我國核電廠進行精細化的地震風險評估提供易損性分析研究基礎。

1 我國場地廣義條件譜與地震動記錄選取

條件均值譜的條件周期只有1個，通常譜型較窄，分析結果可能偏于不保守。Kishida[3]提出了多個條件周期的廣義條件譜理論，該譜譜型克服了條件譜譜型較窄的不足。中國場地廣義條件譜生成需要基于適用于中國場地的危險性分析與向量型概率地震危險性分解結果。

1.1 向量型概率地震危險性分析與分解理論

1.1.1向量型概率地震危險性分析 向量型概率地震危險性分析考慮了強度參數間的相關性，在中國標量型概率地震危險性分析基礎上[7]，考慮譜加速度相關性，中國向量型(兩參數)概率地震危險性分析理論可表示為：

MRDIM1,IM2(im1,im2)=

(1)

式中：MRDIM1，IM2(im1，im2)為IM1和IM2的平均發生率密度；IM1、IM2為兩個強度參數變量；im1、im2為兩個強度參數大小；i為潛在震源區個數變量；M為震級變量；R為距離變量；Θ為方向角變量；m為震級大小；r為距離大小；θ為方向角大小；νi為地震年平均發生率；fIM1,IM2(im1,im2|m,r,θ)為地震動強度參數聯合發生概率密度函數；fM,R,Θ(m,r,θ)為震級、距離和方向角聯合發生概率密度函數。

IM1和IM2的聯合平均發生超越概率λIM1，IM2為：

(2)

式中：MRDIM1，IM2(u1，u2)為IM1和IM2的平均發生率密度；u1、u2分別為兩個強度參數大小。

中國向量型概率地震危險性分析中的強度參數相關性需要考慮中國場地地震動發生時空不均勻性特點，同時中國地震動預測方程通常采用長短軸預測方程的形式，不同潛在震源區中方向角不同，所以除了震級和距離之外，還需要對方向角進行積分運算。對于適用于中國場地的譜型相關性模型，Ji等[8]運用中國地震動數據生成了適用于中國場地的譜加速度相關性模型，王曉磊等[9]給出了基于蒙特卡羅模擬方法的中國向量型概率地震危險性分析程序。

1.1.2向量型概率地震危險性分解 適用于中國場地條件的向量型概率地震危險性分解理論可表示為單位區間震級、距離和方向角下向量型概率地震危險性與總的向量型概率地震危險性之比，其中單位震級、距離和方向角下的向量型強度參數年平均發生率MRDIM1,IM2,x,y,z(im1,im2)可表示為：

MRDIM1,IM2,x,y,z(im1,im2)=

(3)

式中：θz-1～θz為方向角的單位區間；ry-1～ry為距離的單位區間；mx-1～mx為震級的單位區間；x、y、z為單位區間。

單位震級、距離和方向角條件下的向量型概率地震危險性可表示為：

λIM1,IM2,x,y,z=

(4)

式中：λIM1，IM2，x，y，z為單位震級、距離和方向角條件下的IM1和IM2聯合發生超越概率；MRDIM1,IM2,x,y,z(u1,u2)為單位震級、距離和方向角條件下的IM1和IM2平均發生率密度。

參照標量型概率地震危險性分解理論[7]，兩個參數聯合發生條件下，單位區間震級、單位區間距離和單位區間方向角發生概率，即向量型概率地震危險性分解可表示為：

(5)

式中：λIM1，IM2為IM1和IM2聯合發生超越概率；P為事件發生概率；xN、yN、zN為單位區間。

基于上述分解結果，可得到設定地震為：

(6)

1.2 廣義條件譜

1.2.1廣義條件譜理論 Kishida[3]給出了廣義條件譜理論，基本假設是地震動強度參數符合多元正態分布，具體理論可總結如下。

強度參數向量IM可表示為：

IM={IMc,IMs}

(7)

式中：IMc為條件強度參數；IMs為預測強度參數。

強度參數均值向量μ可表示為：

μ=(μc,μs)

(8)

式中：μc為條件強度參數均值；μs為預測強度參數均值。

強度參數協方差矩陣Σ可表示為：

(9)

式中：Σcc為條件強度參數協方差矩陣；Σcs為條件強度參數和預測強度參數協方差矩陣；Σsc為Σcs的對稱矩陣；Σss為預測強度參數協方差矩陣。

(10)

式中：ε為譜型參數向量；σc為譜型參數方差。

(11)

譜型參數向量與譜型參數方差乘積可表示為：

εσc=lnimc-μc

(12)

式中，imc為條件強度參數。

1.2.2我國場地廣義條件譜生成步驟 根據中國場地地震動環境和危險性分析的特點，中國場地廣義條件譜生成需要考慮中國場地地震動預測方程和地震動環境特點，主要包括：1) 基于中國場地概率地震危險性分解結果，生成適用于中國場地地震動預測方程求解的向量型設定地震，包括震級、距離和方向角；2) 考慮中國地震動環境特點，采用適用于中國場地的譜加速度相關系數模型。

中國場地廣義條件譜可基于以下步驟生成：1) 基于場地危險性信息，進行中國場地概率地震危險性分析，計算指定超越概率強度參數；2) 基于場地危險性信息，進行中國場地向量型概率地震危險性分解，生成指定超越概率強度參數向量型設定地震；3) 采用適用于中國場地的譜型相關性模型，利用中國廣義條件譜理論公式，基于概率地震危險性分析和向量型概率地震危險性分解結果，生成場地廣義條件譜。

1.3 基于廣義條件譜理論的地震動記錄選取

地震動初步挑選的誤差平方和SSE可表示為：

(13)

平均值誤差ERRmean和標準差誤差ERRstd[10]可分別表示為：

(14)

(15)

式中：mln imj為選取地震動強度參數對數的平均值；μln imj為目標強度參數對數的平均值；sln imj為選取地震動強度參數對數的標準差；σln imj為目標強度參數對數的標準差。

地震動選取的平均值和標準差誤差組合值SSEs為：

w(sln imj-σln imj)2]

(16)

式中，w為誤差權重系數。

基于廣義條件譜的地震動選取方法步驟為：1) 選定廣義條件譜為目標譜(式(8)、(9))；2) 以步驟1中廣義條件譜為目標譜，采用蒙特卡羅模擬方法生成模擬譜；3) 確定備選數據庫(如美國NGA-West2)；4) 基于震級等地震學參數，對地震動數據進行初篩；5) 采用式(13)，基于最小誤差原則初步選取地震動；6) 基于式(14)、(15)判斷步驟5所選地震動是否滿足選取誤差要求；7) 如果步驟6滿足要求，輸出選取結果；8) 如果步驟6不滿足要求，基于式(16)采用優化算法重復上述選取過程，直到滿足步驟6誤差要求，輸出最終選取地震動記錄。

2 多元地震易損性分析

2.1 多元地震需求分析

結構多元概率地震需求模型可表示為：

GD|IM(d|im2,…,imn)=

P[D≥d|IM1=im1,…,IMn=imn]=

(17)

式中：GD|IM(d|im2,…,imn)為多元概率需求模型；P[D≥d|IM1=im1,…,IMn=imn]為多個地震動強度參數條件下，需求參數超越概率；D為需求參數變量；d為需求參數大小；mD|IM為結構地震需求中位值；βD|IM為地震需求模型對數標準差；Φ為標準正態累積分布函數。mD|IM與強度參數IM之間服從冪指數回歸關系，可表示為：

mD|IM=aIMb

(18)

式中，a和b為系數。

對式(18)兩邊取自然對數：

lnmD|IM=β0+β1ln IM

(19)

式中，β0、β1為回歸系數，一般基于回歸分析得到。

式(19)為一元強度參數IM與工程需求參數EDP對數線性關系模型，二元強度參數IM與EDP對數線性關系模型可表示為：

lnmD|IM=β0+β1ln IM1+β2ln IM2

(20)

式中，β2為回歸系數。

地震需求標準差可表示為：

(21)

式中，D84%與D16%分別為條帶84%與16%的分位值。

2.2 多元地震易損性分析

多元地震易損性函數可表示為：

FR(im1,…,imn)=

P[D≥C|IM1=im1,…,IMn=imn]=

(22)

式中，FR(im1,…,imn)為多強度參數的易損性函數；C為抗震能力；mR為能力中位值；βR為地震易損性函數對數標準差；N為某地震動強度下結構響應總數；n為達到或超過某極限狀態點的個數。式(22)為基于蒙特卡羅模擬方法的多元易損性分析的基礎。

3 我國某核電廠安全殼多元地震易損性分析

3.1 我國某核電廠廠址地震危險性分析與分解

3.1.1算例廠址地震危險性信息 算例廠址信息列于表1，潛在震源區分布圖如圖1所示，地震統計區范圍為東經109°～116°、北緯19°～24°，主要包含32個潛在震源區。

表1 地震統計區參數值Table 1 Parameter value of seismic statistical zone

采用文獻[11]給出的華南地區地震動預測方程，可表示為：

lgY=C1+C2M+

C3lg(R+C4exp(C5M))+σlg Yε

(23)

式中：C1、C2、C3、C4和C5分別為方程系數，具體數值參照文獻[11]；σlg Y為不確定性標準差；ε為譜型參數。

3.1.2算例廠址標量型和向量型概率地震危險性分析與分解 針對算例廠址，基于算例廠址地震危險性信息，采用中國標量型概率地震危險性分析[7]和向量型概率地震危險性分析方法，可得到地震危險性曲線和地震危險性曲面，由于本文安全殼前兩階平動周期與0.24 s和0.07 s接近，所以本文選擇Sa(0.24 s)和Sa(0.07 s)為危險性曲面組合強度參數，危險性曲線和曲面如圖2所示，可發現相同強度Sa(0.07 s)的超越概率比Sa(0.24 s)的大，Sa(0.07 s)和Sa(0.24 s)的聯合危險性曲面上，相同Sa(0.24 s)的危險性曲線比相同Sa(0.07 s)的更陡。圖2中，PGA為峰值地面加速度。

圖1 潛在震源區分布圖Fig.1 Distribution map of potential seismic zone

圖2 算例廠址地震危險性曲線(a)和曲面(b)Fig.2 Seismic hazard curve (a) and surface (b) for example site

基于標量型概率地震危險性分解方法[7]，可得到算例廠址標量型概率地震危險性分解結果，分別給出強度參數Sa(0.24 s)和Sa(0.07 s)年超越概率為萬分之一的地震危險性分解結果，同時基于式(5)可得到算例廠址向量型概率地震危險性分解結果，強度參數分別選擇對應Sa(0.24 s)和Sa(0.07 s)萬年一遇強度的向量型分解結果，如圖3所示。由圖3可見，3個分解結果都不相同。同時基于標量型設定地震[7]和向量設定地震(基于式(6))計算方法生成了場地向量型設定地震(表2)，同樣發現3個設定地震也不相同。

a——標量型Sa(0.07 s)分解；b——標量型Sa(0.24 s)分解；c——向量型Sa(0.07 s)和Sa(0.24 s)分解圖3 地震危險性分解結果Fig.3 Seismic hazard disaggregation result

表2 設定地震Table 2 Scenario earthquake

3.2 算例廠址廣義條件譜

基于式(7)～(12)生成了算例廠址廣義條件譜，如圖4所示，可發現兩個條件強度參數大小與一致危險譜一致，其他部分周期強度參數和一致危險譜不同；兩個條件強度參數的不確定性標準差為0，距離條件周期越遠，不確定性標準差越大。

圖4 生成的廣義條件譜與一致危險譜Fig.4 Uniform hazard spectrum and generalized conditional spectrum generated

3.3 基于廣義條件譜的算例廠址安全殼多元地震易損性分析

3.3.1算例廠址安全殼模型信息 本文以我國華南地區某核電廠安全殼為算例，安全殼集中質量梁單元模型如圖5所示，節點和單元信息列于表3，材料參數列于表4，該安全殼模型前5階模態分析結果列于表5。

圖5 集中質量梁單元模型Fig.5 Lumped mass and beam element model

表3 安全殼和筏板基礎集中質量梁單元模型節點和單元Table 3 Node and element of lumped mass and beam element model of containment and raft foundation

表4 材料參數Table 4 Material parameter

表5 模態分析結果Table 5 Modal result

安全殼是核電廠最后一道防線，在極端災害作用下，其完整性能力可有效保障放射性物質的泄漏。所以通常假設安全殼在災害作用下達到塑性狀態時，安全殼達到失效狀態。安全殼結構可由三線性骨架曲線表示(圖6)[12]，假設當安全殼剪應力達到第1條線拐點時，安全殼達到塑性狀態，即安全殼失效，相應極限狀態對應的剪應力τ1可表示為：

(24)

圖6 三線性骨架曲[12]Fig.6 Trilinear skeleton curve[12]

式中：FC為混凝土抗壓強度；σV為豎向壓應力。

3.3.2基于算例廠址廣義條件譜的地震動記錄選取 基于生成的算例廠址廣義條件譜，采用1.3節地震動記錄選取步驟，可選取30組地震動記錄，地震動選取結果如圖7所示，可發現選取地震動能夠匹配目標譜。

a——模擬譜；b——選取地震動記錄反應譜；c——選取地震動反應譜均值與目標譜均值；d——選取地震動標準差與目標譜標準差圖7 選取地震動結果Fig.7 Result of selected ground motion

3.3.3基于廣義條件譜的算例廠址安全殼多元地震易損性分析 基于以廣義條件譜為目標譜選取的30組地震動，進行多元地震需求分析，基于式(20)進行回歸分析，得到相關系數，結果列于表6。

表6 二元回歸系數Table 6 Binary regression coefficient

基于增量動力分析方法，采用式(22)可計算得到安全殼結構多元地震易損性曲面，如圖8所示。圖8中，Pf為失效概率，多元地震動強度參數分別選用Sa(0.07 s)和Sa(0.24 s)。

圖8 安全殼三維易損性曲面Fig.8 3D fragility surface of containment

指定易損性失效概率，分別做地震易損性曲面等高線，如圖9所示。通過等高線可發現，安全殼模型對于Sa(0.24 s)更為敏感。同時以Sa(0.24 s)為條件參數，做Sa(0.07 s)的地震易損性曲線，如圖10所示；以Sa(0.07 s)為條件參數，做Sa(0.24 s)的地震易損性曲線，如圖11所示。由圖10、11可見：安全殼結構對于Sa(0.24 s)較Sa(0.07 s)更為敏感；安全殼地震易損性結果對Sa(0.24 s)和Sa(0.07 s)都較為敏感，對地震易損性評估結果都有較大影響。不同強度參數條件下地震易損性失效概率列于表7，可發現Sa(0.24 s)較Sa(0.07 s)更為敏感，如Sa(0.24 s)為0.5g且Sa(0.07 s)為1g對應的失效概率為0.008 7，而Sa(0.24 s)為1g且Sa(0.07 s)為0.5g對應的失效概率為0.058 9，即Sa(0.24 s)強度增加0.5g同時Sa(0.07 s)減少0.5g對應的失效概率增加了5.7倍。

圖9 安全殼易損性曲面等高線Fig.9 Contour of fragility surface of containment

圖10 分別以4個強度參數Sa(0.24 s)為條件的Sa(0.07 s)地震易損性曲線Fig.10 Seismic fragility curve of Sa(0.07 s) conditional at four levels of Sa(0.24 s)

圖11 分別以4個強度參數Sa(0.07 s)為條件的Sa(0.24 s)地震易損性曲線Fig.11 Seismic fragility curve of Sa(0.24 s) conditional at four levels of Sa(0.07 s)

4 結論

本文提出了中國場地廣義條件譜理論，生成了算例廠址廣義條件譜，選取地震動記錄，采用多元地震易損性分析方法，生成核電廠安全殼地震易損性曲面，得到以下結論。

表7 不同強度參數條件下的地震易損性失效概率Table 7 Seismic fragility failure probability under different intensity parameters

1) 中國廣義條件譜生成需要適用于中國場地的譜型相關性系數與中國場地地震動預測方程求解的震級、距離和方向角分解結果。

2) 兩個譜加速度強度參數(第1平動周期和第2平動周期)對核電廠安全殼地震易損性都有較大影響，在核電廠安全殼地震易損性分析中應考慮多個強度參數。

3) 基于增量動力分析等解析地震易損性方法，能夠得到更為精細化易損性分析結果，考慮多個地震動強度參數的地震易損性分析結果，可為更為精細化核電廠地震風險提供研究基礎。