姚令侃, 黃藝丹

(1. 西南交通大學土木工程學院, 四川 成都 610031; 2. 高速鐵路線路工程教育部重點實驗室, 四川 成都 610031; 3. 抗震工程技術(shù)四川省重點實驗室道路與鐵道工程抗震技術(shù)研究所, 四川 成都 610031)

山地系統(tǒng)災變行為自組織臨界性研究

姚令侃1,2,3, 黃藝丹1,2,3

(1. 西南交通大學土木工程學院, 四川 成都 610031; 2. 高速鐵路線路工程教育部重點實驗室, 四川 成都 610031; 3. 抗震工程技術(shù)四川省重點實驗室道路與鐵道工程抗震技術(shù)研究所, 四川 成都 610031)

為了研究山地系統(tǒng)宏觀動力學的整體規(guī)律，基于戴維斯地貌發(fā)育理論,提出處于地貌發(fā)育階段幼年晚期的河谷及壯年期的山地具有自組織臨界性(self organized criticality, SOC)的內(nèi)稟屬性,并建立了基于斯特拉勒積分的大流域地貌發(fā)育階段判別方法.以地震觸發(fā)崩塌滑坡為切入點,通過震區(qū)實震資料分析、元胞自動機模擬、振動臺沙堆模型實驗,提出不同烈度區(qū)地震觸發(fā)崩塌滑坡分布的演化規(guī)律：在Ⅶ度、Ⅷ度、Ⅸ度區(qū),崩塌滑坡規(guī)模與出現(xiàn)頻率之間存在良好的負冪律關(guān)系，在Ⅹ度區(qū),冪律關(guān)系弱化，在Ⅺ度區(qū),這一關(guān)系轉(zhuǎn)為對數(shù)正態(tài)分布. 通過3個案例介紹了SOC理論在巖土體地震擾動深度評估、泥石流防治工程設(shè)計徑流量極值計算、基于地震活動性參數(shù)b值在地應力評估中的應用.

山地系統(tǒng);自組織臨界性;地震觸發(fā)崩塌滑坡;元胞自動機

在非線性物理學中,作為災變理論提出的前沿理論——自組織臨界狀態(tài)(self-organized criticality, SOC)理論,是P Bak等為解釋無序的、非線性復雜系統(tǒng)的行為特征而提出的新概念.這類系統(tǒng)包含著眾多發(fā)生短程相互作用的組元,并自發(fā)地向著一種臨界狀態(tài)進化.在臨界狀態(tài)下,小事件引起的連鎖反應能對系統(tǒng)中大量數(shù)目的組元產(chǎn)生影響,從而導致大規(guī)模事件的發(fā)生.雖然發(fā)生的小事件比大事件多,但遍及所有規(guī)模的連鎖反應是動態(tài)特性的一個必不可少的部分,所有時空關(guān)聯(lián)函數(shù)都是冪律(power-law)的,因此冪律可以作為自組織臨界狀態(tài)的證據(jù)[1].

沙堆模型是SOC的范例.G A Held等采用在圓盤上逐粒加沙的方式構(gòu)造沙堆,當沙堆的傾角在臨界角附近時沙堆停止增長,此時,對新添加沙粒的響應是無法預測的,沙?？赡芄潭ㄔ谏扯焉?也可能引起小范圍沙粒的滑動,還可能導致更大規(guī)模的雪崩(avalanche),但總是呈現(xiàn)崩塌規(guī)模與出現(xiàn)頻率成反比的冪律關(guān)系[2].

SOC的概念已被用于解釋從山脈形成到股市波動等許多復雜現(xiàn)象.地震學是SOC最早關(guān)注的領(lǐng)域,如古登堡-里克特定律表明地震的頻度與震級之間存在冪律關(guān)系.P Bak認為這種冪律來源于地殼系統(tǒng)的SOC,反之,也可把冪律作為地殼被鎖定在永久臨界狀態(tài)的證據(jù)[3].在地學領(lǐng)域研究SOC是具有學科交叉性質(zhì)的前沿研究,從宏觀層面,於崇文院士以完整和獨立的命題提出了固體地球系統(tǒng)的復雜性與自組織臨界性[4],認為地質(zhì)系統(tǒng)是自然界中的一種異常復雜的開放、遠離平衡、相互作用的巨大耗散動力系統(tǒng),具有自組織臨界性的內(nèi)稟屬性,其時空行為服從地質(zhì)作用的自組織臨界過程動力學.從微觀層面,許強等[5]論述了巖石破裂過程的自組織臨界特征.

運用SOC理論認識和描述山地系統(tǒng)災變動力學行為是本團隊的主攻方向.國家自然科學基金委先后立項資助本團隊開展了“散粒體自組織臨界性”(2003)、“大尺度散粒體自組織臨界性及其判據(jù)”(2005)、“地震觸發(fā)崩塌滑坡自組織臨界性”(2012)項目研究,這些題目也反映了本團隊從散粒體材料層面至山地系統(tǒng)層面SOC性質(zhì)的研究歷程.本文主要以地震觸發(fā)崩塌滑坡作為研究山地系統(tǒng)災變行為SOC的切入點,開展一項專題研究,并介紹在應用領(lǐng)域開展的一些探索.

1 山地系統(tǒng)具備SOC性質(zhì)的條件

1.1 SOC系統(tǒng)的必要條件

從現(xiàn)象學的角度,SOC系統(tǒng)應具有的必要條件是[6]:

(1) 系統(tǒng)是耗散的,包含大量(數(shù)目應在數(shù)百萬個以上)發(fā)生短程相互作用的組元;大量自由度以某種均衡態(tài)勢存在,不出現(xiàn)優(yōu)勢自由度而使系統(tǒng)僅有幾個集體自由度來表征的現(xiàn)象.

(2) 組元之間最近鄰位置相互作用,但存在著長程相關(guān)的關(guān)系.

(3) 系統(tǒng)自發(fā)地朝著臨界狀態(tài)進化,并將會永久性地或在一個有意義的時段內(nèi)被鎖定在這個狀態(tài).

(4) 系統(tǒng)同時具備兩種動力學特征:一是敏感性,系統(tǒng)處于臨界狀態(tài),微小的擾動也可能產(chǎn)生遍及全局的連鎖反應;二是魯棒性,系統(tǒng)的整體不受破壞.

(5) 所有的時空關(guān)聯(lián)函數(shù)均為冪次.

從必要條件可看出,SOC系統(tǒng)的基本特征是以臨界狀態(tài)為動力學吸引子.即在外界持續(xù)的物質(zhì)(或能量)供給條件下,廣延耗散動力系統(tǒng)自發(fā)地向臨界狀態(tài)演化,一旦到達臨界狀態(tài),將通過與外界能量的交換,永久性地或在一個有意義的時段內(nèi)能被鎖定在這個狀態(tài).這是從自然系統(tǒng)演化趨勢角度提出的SOC基本條件,其特點是強調(diào)自然系統(tǒng)SOC的時段性.

但以上5條尚不能涵蓋SOC系統(tǒng)需要滿足的所有條件,還需要補充一些能夠引起系統(tǒng)產(chǎn)生自組織到臨界態(tài)的趨勢的特定特征.

從構(gòu)造性的角度,SOC行為只會在緩慢驅(qū)動的、相互作用占主導地位的閾值系統(tǒng)中出現(xiàn).閾的存在使系統(tǒng)積累能量、應力和物質(zhì)直至達到不穩(wěn)定閾值.SOC 是慢驅(qū)動機制,即能量積累速率緩慢,而能量耗散的時間尺度遠小于前者.以地震為例,地殼處于運動中,因而經(jīng)受緩慢的形變,形變將導致巖石中應力的積累.如果在某些地方巖石不能再承受應力,就會發(fā)生突發(fā)事件,并向周圍環(huán)境釋放能量.如果環(huán)境可以承擔這種能量的增加,那么這是一個孤立事件;如果環(huán)境也處在一種臨界狀態(tài),就可能發(fā)生鏈式反應,地震就相當于應力釋放的鏈式反應,涉及到比較大的區(qū)域.這兩種過程的時間標度差別巨大,如“5·12”汶川大地震,專家估計應力在龍門山推覆構(gòu)造帶上的積累過程可能長達千年,地震時應力釋放的主要過程僅在約100 s內(nèi)完成[7].

1.2 基于戴維斯地貌發(fā)育理論的山地系統(tǒng)臨界狀態(tài)思辨

美國著名地理學家戴維斯(W M Davis)于1899年首次創(chuàng)立了侵蝕循環(huán)學說(theory of the cycle of erosion),認為地塊從開始上升到被逐漸剝蝕夷平,直至降低到起伏不大的地面或者接近基準面的準平原之間,存在著連續(xù)的、同時又有階段性的剝蝕過程和地表形態(tài).在地表發(fā)育的過程中,Davis強調(diào)構(gòu)造、作用和時間(侵蝕階段)這3個要素之間的相互作用影響[8],進而將循環(huán)過程中的地形發(fā)展分為:地形起伏不大,河間地廣闊平坦的幼年期;地面主要由谷坡和狹窄的分水嶺組成的壯年期;具有殘丘的準平原的老年期(圖1).

(a)最初,地形起伏和緩,流水不暢(b)幼年早期,溝緣狹窄,高地寬闊平坦(c)幼年晚期,巖坡為主,仍有溝緣,平坦高地(d)壯年期,多為巖坡與狹窄的分水嶺(e)壯年晚期,地形起伏較緩,谷底寬展(f)老年期,成為具有蝕余殘山的準平原(g)再次構(gòu)造抬升,進入第二循環(huán),重現(xiàn)幼年早期圖1 侵蝕循環(huán)示意[9]Fig.1 Schematicoferosioncycle[9]

按照侵蝕循環(huán)理論,地貌發(fā)育循環(huán)過程中的地形變化如圖2所示，圖中，曲線L1、L2分別代表原始地面河谷的平均高度和分水嶺的平均高度隨時間的變化.虛線表示地殼短暫快速的抬升階段,H1表示原始地形的起伏程度.

由圖2可知,整個幼年期階段直至壯年中期,河谷強烈迅速下切,河谷平均海拔持續(xù)下降,河谷之間由比較寬廣的分水地發(fā)展為分水嶺,平均海拔高度沒有顯著變化,因此,原始地面的起伏程度迅速增加,至壯年中期,起伏程度達到最大,分水嶺與河谷之間的高差達到最大值H2.壯年期晚期和老年期階段,河谷已經(jīng)接近侵蝕基準面,高程基本沒有變化,分水嶺受到外營力的侵蝕而不斷降低,分水嶺變得渾圓低矮,直至準平原.

圖2 Davis河流地貌發(fā)育圖式(根據(jù)W. M. Davis)Fig.2 The sequence in the developmental changes of landforms of an ideal geographical cycle by W. M. Davis’s theory

地貌的形成和發(fā)展是內(nèi)、外營力相互作用的結(jié)果.內(nèi)營力趨向于使山體隆升,增強區(qū)域起伏程度,即使山體愈發(fā)陡峻;外營力趨向于使山體高度降低、削平,減弱區(qū)域起伏程度,即使山體愈發(fā)渾圓低矮.內(nèi)營力是系統(tǒng)的,來自全球性的板塊構(gòu)造,可以認為其速度是確定的,而外營力是變化的.在新構(gòu)造運動強烈的山區(qū),山地隆升速度每年可達幾毫米至十幾毫米,山體剝蝕降低的速度小于地塊隆升速度;但隨著山體高度的上升,帶來侵蝕基準面的降低,進而加大地表物質(zhì)的重力作用和水流的下切侵蝕與搬運作用,外營力效應隨之增強.外營力效應絕對量的增加必將導致內(nèi)外營力作用的相對差距縮小,當?shù)孛惭莼幱趦?nèi)外營力相當?shù)碾A段時,由于二者對山地地貌塑造的反向效應,山體坡度最大只能達到一個特定值,即所謂的臨界坡度(對應于圖2的H2時期),按照SOC理論,臨界坡度成為該階段山地斜坡系統(tǒng)演化的吸引子.

對照圖1,雖然在一個侵蝕循環(huán)內(nèi),地貌將經(jīng)歷幼年期、壯年期、老年期,流域內(nèi)的坡體相應經(jīng)歷向臨界坡度發(fā)展、達到臨界坡、偏離臨界坡的演變過程,但處于地貌發(fā)育階段幼年晚期地塊的河谷(圖1(c)),以及壯年期的山地(圖1(d)),其斜坡系統(tǒng)能夠保持在臨界坡度,即在這一時段內(nèi)已經(jīng)演化到了臨界狀態(tài).

若從構(gòu)造性的角度,在漫長的地質(zhì)過程中,地表形態(tài)的變化是組成地表的物質(zhì)在漸變與突變的不斷轉(zhuǎn)化中進行的.漸變是指地形在漫長時間完成的變化過程,突變則是在急促的短時間內(nèi)完成的變化過程.由于板塊運動,在長達千百萬年的時間尺度上使山脈隆升,將能量(勢能)緩慢而持續(xù)地輸入斜坡系統(tǒng),使斜坡系統(tǒng)能量積累到臨界狀態(tài),這一過程是非隨機的;而在外營力隨機擾動條件下,斜坡物質(zhì)可能失穩(wěn)并在重力作用下向下輸移,即斜坡系統(tǒng)能量又通過突發(fā)性的山地災害的方式耗散,這一過程是隨機的.山地隆升驅(qū)動斜坡系統(tǒng)構(gòu)造的速度(cm/a)和崩塌滑坡滑動速度(m/s)之間出現(xiàn)時間標度的巨大分離,相互作用為主導、閾動力學和慢驅(qū)動的有機結(jié)合,使得斜坡系統(tǒng)的演化成為自然界中SOC的良好實例.

綜上所述,無論從現(xiàn)象學的角度,還是從構(gòu)造性的角度,對照SOC的基本定義,可以認為當且僅當處于地貌發(fā)育階段幼年晚期的河谷,如橫斷山三江并流區(qū)的高山峽谷地段(圖3、圖4),以及壯年期的山地,如龍門山(圖5),山地斜坡系統(tǒng)成為這些區(qū)域地貌的主體,斜坡系統(tǒng)已經(jīng)演化到了臨界狀態(tài),SOC成為系統(tǒng)的內(nèi)稟屬性,它的災變行為將服從SOC過程動力學.特別是崩塌、滑坡、泥石流等斜坡重力作用類地表過程,其共同特征是能量耗散以斜坡物質(zhì)失穩(wěn)下滑實現(xiàn)的,從物理概化角度,可認為沙堆模型是一種能抓住真實斜坡系統(tǒng)基本特征的簡單理想模型,因而它們的動力學特征都應能在SOC 的概念框架下得到解釋.

圖3 三江并流區(qū)怒江河谷已演化到臨界狀態(tài)Fig.3 Nujiang river valleys in the three parallel rivers region having evolved to the critical state

圖4 三江并流區(qū)瀾滄江河谷已演化到臨界狀態(tài)Fig.4 Lancang river valleys in three parallel rivers region having evolved to the critical state

圖5 龍門山已演化到臨界狀態(tài),在強震作用下大面積崩塌Fig.5 Longmen mountains having evolved to the critical state, with lots of landslides triggered by earthquake

1.3 基于斯特拉勒積分的山地系統(tǒng)發(fā)育階段判別方法

以上提出了以地貌發(fā)育階段作為界定山地系統(tǒng)具有SOC性質(zhì)的適用范圍,應用中還必須建立地貌發(fā)育階段的定量判別方法.

1952年,斯特拉勒(A N Strahler)提出侵蝕流域的面積-高程分析方法[10],可以定量推求Davis的地貌發(fā)育階段.面積-高程曲線分析法是描述一定高度范圍內(nèi)的面積隨相對高度變化所表示的曲線及其所圍成的面積,相對高度可以確定侵蝕過程的強度,而殘留的面積可以代表這種強度下地貌的保持能力,因此,可以說面積-高程曲線提供了地貌的發(fā)育信息[11].

記流域內(nèi)等高線的值和最低點之間的高差為h,每條等高線以上的面積為a,全流域面積為A0,流域內(nèi)最高點和最低點之間的高差為H,分別以

X=a/A0,

(1)

為橫坐標和縱坐標畫圖,可以得到曲線

Y=h/H,

(2)

也就是面積-高程曲線,稱為斯特拉勒曲線(the Strahler’s curve).

設(shè)定積分

(3)

式中:

S為斯特拉勒曲線與坐標軸包圍的面積,稱為斯特拉勒積分(the Strahler’s integral).

用S值推求侵蝕流域地貌演化階段,即:

S>0.60,幼年期;

0.35≤S≤0.60,壯年期;

S<0.35,老年期.

由于斯特拉勒積分為無量綱參數(shù),因此,曲線可以描述和比較不同規(guī)模的流域,但必須是在流域內(nèi)地貌處于同一發(fā)育階段的前提下,所以,斯特拉勒積分應用于小流域一般沒有問題.然而,對于大型山區(qū)河流,不同河段可能處于不同的地貌單元,這與Davis循環(huán)理論中流域內(nèi)的各個部分同時發(fā)育的假定是相悖的,不能直接應用該方法.

作者提出包含不同地貌單元的大流域斯特拉勒積分計算方法:

首先確定計算區(qū)域,通過選擇河道控制點(也就是計算流域的出口點),確定適用于Davis侵蝕循環(huán)理論的范圍.選擇河道控制點的原則為:

① 控制點以上的流域是災勢評估的研究區(qū);

② 河道控制點以上的流域應基本上處于同一地貌單元.

然后利用DEM和ArcGIS技術(shù),獲得斯特拉勒曲線的橫、縱坐標值.具體方法是在具有規(guī)則格網(wǎng)的DEM上,對研究區(qū)域內(nèi)所有大于某個高程值的柵格單元個數(shù)進行累加,再乘以柵格單元的面積A0,即得到該高程值以上的面積[12],如式(4)所示.

(4)

式中:

Ah0為等高線值為h0以上的面積;

Nh為研究區(qū)域內(nèi)高程值為h的柵格個數(shù).

以(Ah0/A,h0/H)為橫、縱坐標值,繪制斯特拉勒曲線,可計算得到斯特拉勒積分.

現(xiàn)舉例說明基于斯特拉勒積分的山地系統(tǒng)發(fā)育階段判別方法.2010年4月14日,我國青海省玉樹縣發(fā)生Ms7.1級地震，影響流域主要水系包括通天河、扎曲、巴曲等,地形以高海拔、低起伏為主.震區(qū)位于通天河和金沙江的交匯地段,圈定的研究區(qū)域的流域包括沱沱河、通天河、金沙江上游在內(nèi)的長江上游河段,研究區(qū)域的水系格局如圖6所示.

圖6 研究流域及周邊流域、裂點位置、計算流域區(qū)域圖Fig.6 Region graphs of the study basin and its surrounding basin, the knick point location, and the computation basin

由于青藏高原的分階段、非均一隆升,使高原東緣的外流河產(chǎn)生溯源侵蝕,并形成新的河谷,它與未被溯源侵蝕的老河谷交替的地方,河床坡度突然增加,形成裂點[13],裂點上下游的河谷往往處于不同的地貌發(fā)育狀態(tài).通過提取河道的縱剖面圖(圖7)可以看出,在海拔4 000 m的地方存在裂點.目前青藏高原普遍存在著高原面、盆地面兩級夷平面.盆地面是青藏高原形成于上新世紀初至上新世紀末的第二級夷平面,盆地面內(nèi)平坦開闊,切割微弱.位于裂點以上的金沙江上游谷地是溯源侵蝕尚未達到的地方,這一地區(qū)的盆地面保存完整,較少切割.裂點以下的金沙江河谷位于我國地形最為陡峻的橫斷山區(qū),河流切割成深邃的峽谷.裂點上、下游兩個地貌單元差異明顯[14].又因玉樹地區(qū)處于該裂點的上方,可滿足上文提出的選擇河道控制點的兩條原則，故決定選擇該裂點為計算流域的出口點,繪制斯特拉勒曲線,計算斯特拉勒積分S=0.32,如圖8所示.圖8表明,震區(qū)內(nèi)流域?qū)儆趬涯晖砥诎l(fā)育階段,寬谷緩丘是主要地貌特征.在玉樹地震時,地震觸發(fā)的崩塌滑坡無論從規(guī)模還是數(shù)量上,都遠小于汶川地震和蘆山地震,認為地貌發(fā)育階段不同是玉樹地震崩塌滑坡災害小于汶川、蘆山地震的主要因素(汶川和蘆山地震區(qū)均處于壯年期).

圖7 研究河流縱剖面圖Fig.7 Longitudinal profile of the study river

圖8 計算區(qū)域的斯特拉勒曲線Fig.8 Strahler curve of the computation basin

反之,若按圖6所示的全河段計算相對應流域的斯特拉勒積分,得S′=0.53,據(jù)此推斷整個流域的發(fā)育期處于壯年期,應該是坡陡谷深的地貌,顯然與震區(qū)的自然景觀不符.至此,解決了斯特拉勒分析方法只能適應于小流域的問題,建立了山地系統(tǒng)發(fā)育階段的定量判別方法.

2 不同烈度區(qū)地震觸發(fā)崩塌滑坡分布規(guī)律

發(fā)生在四川龍門山地區(qū)的“5·12”汶川地震、“4·20”蘆山地震,是人類有現(xiàn)代觀測儀器以來,地震觸發(fā)崩塌滑坡數(shù)量最多、資料最翔實的兩次大地震,震后有關(guān)地震觸發(fā)崩塌滑坡的研究成為熱點課題,但國內(nèi)外大多數(shù)研究集中在對該場地震崩塌滑坡災害的統(tǒng)計、地震觸發(fā)崩塌滑坡形成的力學機理等方面,缺少從物理學層面對其總體特征進行描述的整體理論.

龍門山是青藏高原邊緣山脈中的陡度變化最大的山脈,在30多公里范圍內(nèi)海拔從700 m升高到5 000 m.晚新生代中新世以來,龍門山至少有5～10 km的底層被剝蝕掉,上升速度約達0.6 mm/a.近年來的地形變化資料表明,該構(gòu)造帶的九頂山地區(qū)正以0.3～0.4 mm/a速度持續(xù)上升[15]，這種隆升和夷平的持續(xù)作用,造成龍門山河谷深切、地勢陡峻的地貌景觀.龍門山整體處于地貌演化的壯年早期,斜坡系統(tǒng)已經(jīng)演化到了臨界狀態(tài),這就是能在SOC的概念框架下,研究汶川地震、蘆山地震觸發(fā)崩塌滑坡整體分布規(guī)律的前提條件.為此,選擇地震觸發(fā)崩塌滑坡作為研究山地系統(tǒng)災變行為SOC的切入點,在原型問題的物理代表性、資料擁有量等方面,均具有獨特的優(yōu)勢.

2.1 地震觸發(fā)崩塌滑坡實震資料分析

利用震后遙感影像資料進行人工目視解譯是大面積獲取地震觸發(fā)崩塌滑坡信息的主要方法,由于使用的遙感影像資料精度不同、判識人員的判識標準和經(jīng)驗不同等,對同一區(qū)域的判識可能會出現(xiàn)較大差異,因此,現(xiàn)場調(diào)查工作不可忽視.但以上細節(jié)并不是關(guān)鍵問題,影響地震觸發(fā)崩塌滑坡分布規(guī)律最具控制性的因素是地震烈度,在研究崩塌滑坡分布規(guī)律時,應該按地震烈度區(qū)分別統(tǒng)計.

汶川地震后,獲取2008年6月4日的ALOS衛(wèi)星影像(分辨率為10 m),覆蓋范圍約為9 750 km2(30°58.78′N～32°2.3′N,103°33.97′E～104°36.17′E),包括北川、汶川、茂縣、都江堰等區(qū)域.另外,還通過下載獲得汶川震前TM衛(wèi)星圖片.由于這套ALOS遙感數(shù)據(jù)資料Ⅸ度以下烈度區(qū)云層覆蓋率高,因此,僅對Ⅸ度及以上烈度區(qū)進行崩塌滑坡的遙感解譯工作.在上述區(qū)域內(nèi),采用人機交互的目視解譯技術(shù),輔以野外實地調(diào)查,判譯出地震觸發(fā)崩塌滑坡9 341處,解譯結(jié)果如圖9所示.

蘆山地震后,通過多種途徑收集和購買了震區(qū)航空、航天遙感影像資料,主要包括中國科學院遙感與數(shù)字地球研究所提供的三批航片圖(拍攝時間為2013年4月20日和21日,其中,第一批影像分辨率為0.6 m,第二、三批影像有0.4 m和2 m兩種分辨率);購買的資源三號衛(wèi)星影像(拍攝時間為2013年5月13日,分辨率為2.1 m).對以上遙感影像資料進行幾何糾正、融合、拼接、圖像增強處理后,獲得蘆山震區(qū)遙感數(shù)據(jù)5 655 km2(102°27.68′E～103°30.84′E,29°44.54′N～30°26.91′N),包括邛崍市、名山縣、雨城區(qū)、蘆山縣、天全縣、寶興縣、滎經(jīng)縣等區(qū)域.在上述區(qū)域,對地震觸發(fā)的崩塌滑坡災害點進行逐一解譯,并結(jié)合實地考察檢驗解譯結(jié)果,最終確定地震觸發(fā)崩塌滑坡1 608處,解譯結(jié)果如圖10所示.

圖9 汶川地震觸發(fā)的崩塌滑坡分布示意圖Fig.9 The distribution map of landslides induced by Wenchuan earthquake

按不同烈度區(qū),分別對汶川Ⅸ度、Ⅹ度、Ⅺ度,蘆山Ⅶ度、Ⅷ度、Ⅸ度區(qū)的崩塌滑坡進行統(tǒng)計分析,以崩塌滑坡面積(A)作為規(guī)模的度量,分析崩塌滑坡規(guī)模與發(fā)生頻率之間的關(guān)系.由于蘆山地震各烈度區(qū)的遙感影像資料更為完備(Ⅸ度區(qū)遙感覆蓋面積為100%，Ⅷ度區(qū)遙感覆蓋面積為94.5%，Ⅶ度區(qū)遙感覆蓋面積為59.4%),故蘆山震區(qū)增加了崩塌滑坡點密度這一指標,統(tǒng)計結(jié)果見表1.

由表1可看出：汶川Ⅸ度區(qū),地震觸發(fā)的崩塌滑坡面積與出現(xiàn)頻率之間呈現(xiàn)良好的冪律關(guān)系(R2=0.916);汶川Ⅹ度區(qū),崩塌滑坡面積與出現(xiàn)頻率之間冪律關(guān)系式的相關(guān)系數(shù)下降到0.906,可以認為基本服從冪律關(guān)系;汶川Ⅺ度區(qū),崩塌滑坡面積與出現(xiàn)頻率之間的關(guān)系符合對數(shù)正態(tài)分布;蘆山Ⅶ度、Ⅷ度、Ⅸ度烈度區(qū),崩塌滑坡規(guī)模與發(fā)生頻率之間均服從冪律分布,且冪指數(shù)b值相近.蘆山地震不同烈度區(qū)震害的主要區(qū)別是:隨地震烈度的減小,崩塌滑坡點的密度單調(diào)減小.

汶川地震、蘆山地震是我國迄今為止獲得地震山地災害資料最全、精度最高的大地震,據(jù)此得到不同烈度區(qū)地震觸發(fā)崩塌滑坡的統(tǒng)計結(jié)論,作者更加關(guān)注該結(jié)論是汶川地震、蘆山地震的獨特現(xiàn)象,還是具有普適性意義的規(guī)律.因此,在SOC的概念框架下,利用元胞自動機模擬方法,結(jié)合振動臺沙堆模型物理實驗,對上述實震資料顯示的統(tǒng)計規(guī)律,從物理層面詮釋其機理.

圖10 蘆山地震觸發(fā)的崩塌滑坡分布示意圖Fig.10 The distribution map of landslides induced by Lushan earthquake

2.2 地震觸發(fā)崩塌滑坡元胞自動機模擬

元胞自動機(cellular automata或cellular automaton, CA)是空間和時間都離散,物理參量只取有限數(shù)值集的物理系統(tǒng)的理想化模型[16].該模型以規(guī)則網(wǎng)格形式分布、空間離散的元胞個體為基本單元,元胞遵循一定的演化規(guī)則作同步更新來模擬真實的物理系統(tǒng),特別適合復雜系統(tǒng)時空演化過程的動態(tài)模擬研究.

模型方法在SOC的研究中具有非常重要的地位,目前所有對沙堆生長與坍塌機理的模擬算法都是以元胞自動機為數(shù)學基礎(chǔ),對SOC的理解也大部分來源于元胞自動機沙堆模型的數(shù)值模擬.處于青壯年期的山地系統(tǒng)具有SOC的內(nèi)稟屬性,因此,可以構(gòu)建地震觸發(fā)崩塌滑坡沙堆模型,用元胞自動機來模擬.

根據(jù)原型問題的物理特征,與傳統(tǒng)沙堆模型相比,地震觸發(fā)崩塌滑坡的沙堆模型應具有以下特點[17]:

(1) 傳統(tǒng)沙堆模型考察的是一個沙堆在受到多次擾動下坍塌規(guī)模隨時間的分布規(guī)律.地震發(fā)生時,同時引發(fā)多處崩塌滑坡,需要考察多個沙堆在同時受到一次擾動時坍塌規(guī)模的分布規(guī)律.

表1 不同地震烈度區(qū)地震觸發(fā)崩塌滑坡統(tǒng)計Tab.1 Statistics of landslides in zones with different seismic intensity

(2) 在傳統(tǒng)沙堆模型中,沙堆是局部受到擾動，在地震觸發(fā)崩塌滑坡問題中,坡體是整體受到擾動,因此,模型需要采用系統(tǒng)整體受擾的方式.

(3) 傳統(tǒng)沙堆模型均是在同一微擾條件下觀察系統(tǒng)的動力學行為.而地震觸發(fā)崩塌滑坡,首先不同烈度區(qū),地震對坡體的震動力是不同的,其次地震力相對坡體已超過微擾量級,因此,設(shè)計不同擾動強度,來模擬這一物理現(xiàn)象.

(4) 傳統(tǒng)沙堆模型,外界對系統(tǒng)的輸入是物質(zhì)(添加沙粒),發(fā)生連鎖反應時,擾動的傳播是物質(zhì)的傳播,遵守物質(zhì)守恒原則.在地震觸發(fā)崩塌滑坡的問題中,外界對坡體系統(tǒng)輸入地震力,坡體失穩(wěn)需要克服其自穩(wěn)能力,會消耗一部分輸入的能量,擾動的傳播是能量的傳播,并且存在著能量耗散.

根據(jù)以上特點,構(gòu)造地震觸發(fā)崩塌滑坡元胞自動機沙堆模型,基本算法步驟如下:

步驟1 生成沙堆.地震觸發(fā)崩塌滑坡相當于考察大量沙堆,同時受到一次地震作用時坍塌規(guī)模的整體分布規(guī)律,因此,一次性生成N個沙堆,每個沙堆規(guī)模相等,但初始狀態(tài)不同.對于每個沙堆而言,考慮一個L×L的二維系統(tǒng),用(i,j)代表元胞所處的位置(其中, 1≤i,j≤L),每一個元胞有上、下、左、右4個鄰居,Fi,j為反映元胞(i,j)穩(wěn)定性的狀態(tài)值(相當于元胞的能量).所有元胞賦予一個0到閾值Fth之間的初始值,且取隨機數(shù).

步驟2 沙堆演化到臨界態(tài).每個沙堆按照下列規(guī)則連續(xù)反應,直到所有沙堆演化到臨界狀態(tài).

① 找到最大狀態(tài)值Fmax的元胞,把每個元胞的狀態(tài)值都增加Δ=Fth-Fmax(相當于對整個系統(tǒng)的一個擾動),即Fi,j→Fi,j+Δ;

② 若新的Fi,j大于或等于設(shè)定的閾值Fth,即Fi,j≥Fth,則該元胞倒塌,并向鄰居傳播擾動,重新分配Fi,j給它的4個最近鄰:

Fi±1,j→Fi±1,j+αFi,j,

Fi,j±1→Fi,j±1+αFi,j,

Fi,j→0;

(5)

③ 這種擾動傳播可能會導致鏈式反應,如果由于元胞(i,j)的倒塌導致它的鄰居變得不穩(wěn)定,重復②直至所有的元胞狀態(tài)值小于閾值(Fi,j

通過步驟2,產(chǎn)生N個均處于臨界狀態(tài),但各元胞狀態(tài)值不同的沙堆.

步驟3 對N個沙堆同時施加一次擾動.令F′為擾動強度,每個沙堆的元胞狀態(tài)值都統(tǒng)一增加F′,即Fi,j→Fi,j+F′.元胞之間的相互作用仍按照規(guī)則②、③執(zhí)行.

步驟4 改變擾動強度F′,重復步驟3,獲得沙堆在不同擾動強度下坍塌規(guī)模與發(fā)生頻率的關(guān)系.

需要說明,按照Δ=Fth-Fmax施加擾動時,一般只會直接觸發(fā)一兩個元胞,但擾動強度加大后,可能會觸發(fā)一批元胞,它們各自都可能引發(fā)鏈式反應,這些鏈式反應在空間上也會有交叉現(xiàn)象.因此,在算法上采用并行處理,所有受擾元胞在同一時步內(nèi)按照平行更新的方式反應,通過記錄倒塌的元胞數(shù)目,作為沙堆坍塌規(guī)模的度量.

模型參數(shù)取為:α=0.2,Fth=1,L=50.

實驗生成100萬個沙堆(N=106),先連續(xù)反應106次(均取F′=1-Fmax),以確保沙堆演化到臨界狀態(tài).然后取

F′=0.000 01,0.000 02,0.000 04,0.000 08,

1-Fmax,0.001,0.005,0.01,

開展了8組模擬實驗.令沙堆坍塌規(guī)模為S,分析坍塌規(guī)模與發(fā)生頻率的關(guān)系,雪崩密度等于發(fā)生坍塌事件的沙堆數(shù)除以總沙堆數(shù),實驗結(jié)果見表2.

表2 元胞自動機模擬實驗結(jié)果Tab.2 Results of cellular automata simulation

由表2可看出,以擾動強度1-Fmax為界可以劃分為兩個區(qū)間,沙堆模型的動力學特性呈現(xiàn)出不同的性質(zhì):在擾動強度小于1-Fmax的區(qū)間,當F′從0.000 08逐漸減小到0.000 01時,隨擾動強度的降低,坍塌規(guī)模與發(fā)生頻率基本服從同一冪律分布(表2),但隨擾動強度的減小,雪崩密度單調(diào)減小(圖11);在擾動強度大于1-Fmax的區(qū)間,當F′從0.001逐漸增加到0.01時,隨擾動強度的增加,沙堆模型的動力學特性將經(jīng)歷冪律—冪律弱化—對數(shù)正態(tài)分布的漸進式的演變過程(圖12).

對照兩次地震的實震資料統(tǒng)計結(jié)論,汶川地震Ⅸ度區(qū)崩塌滑坡規(guī)模與出現(xiàn)頻率之間的關(guān)系可用冪律描述,Ⅹ度區(qū)這一關(guān)系基本服從冪律分布,Ⅺ度區(qū)更偏向于對數(shù)正態(tài).蘆山地震Ⅶ-Ⅸ度烈度區(qū)地震觸發(fā)崩塌滑坡規(guī)模與頻率均服從冪律分布,但隨地震烈度的減小,崩塌滑坡點密度單調(diào)減小.通過以上元胞自動機模擬,證明了不同烈度區(qū)地震觸發(fā)崩塌滑坡的分布規(guī)律存在必然的聯(lián)系性,這些現(xiàn)象均來源于山地斜坡系統(tǒng)具有SOC內(nèi)稟屬性的物理機制,從而實現(xiàn)了從物理角度對“5·12”汶川地震、“4·20”蘆山地震不同烈度區(qū),地震觸發(fā)崩塌滑坡分布規(guī)律演變機理的詮釋,初步建立了不同烈度區(qū)地震觸發(fā)崩塌滑坡總體特征的完整表征體系.

2.3 地震觸發(fā)崩塌滑坡振動臺沙堆模型實驗

元胞自動機數(shù)值模擬是獲得SOC系統(tǒng)性質(zhì)的主要途徑,但相對地震觸發(fā)崩塌滑坡的原型問題,無疑是一種高度概化的處理,因此有必要配合沙堆模型物理實驗相互印證.

利用大型地震模擬振動臺,輸入精確可控的地震波,開展動力擾動下的沙堆模型實驗,使之更加逼近物理原型.實驗設(shè)備為西南交通大學高速鐵路線路工程教育部重點實驗室的電液伺服驅(qū)動地震模擬振動臺.在振動臺上放置一長2.58 m、寬1.50 m、高1.95 m的箱體,選用粒徑為0.6～50.0 mm經(jīng)過篩分的干燥天然沙石在箱體里堆成一個單面坡沙堆,總重量達6.8 t(圖13).

圖11 雪崩密度ρ與擾動強度F′Fig.11 Avalanche density ρ vs disturbance intensity F′

圖12 不同擾動強度下坍塌規(guī)模的概率密度曲線Fig.12 Probability density curves for different disturbance intensities

圖13 沙堆模型Fig.13 Sandpile model

向振動臺輸入汶川地震臥龍臺站記錄的修正波(圖14),模擬地震對沙堆的擾動.取峰值加速度為0.075g、0.100g、0.125g、0.150g、0.200g、0.350g、0.450g,開展了7組實驗.

圖14 汶川地震臥龍臺記錄修正波Fig.14 Modified Wenchuan acceleration wave recorded by Wolong station

每組實驗步驟如下:

步驟1 在箱內(nèi)堆成單面坡沙堆,坡面靠沙粒的重力下滑自然形成,使坡腳達到開口端邊緣,沙堆達到臨界坡度;

步驟2 輸入指定峰值加速度,稱量從坡腳處落入落沙收集槽的沙粒重量,作為該組的一次實驗值;

步驟3 每次實驗后,將落沙收集槽內(nèi)的沙粒補充回沙堆上,確保沙堆始終處于臨界狀態(tài);

步驟4 重復步驟2、3,每組實驗不少于60次.

以落沙收集槽內(nèi)沙粒的重量作為規(guī)模的度量,令沙粒重量為x,分析坍塌規(guī)模與發(fā)生頻率之間的關(guān)系,坍塌密度等于發(fā)生坍塌事件的實驗次數(shù)除以總實驗次數(shù),實驗結(jié)果見表3.

由表3可看出,當輸入地震波峰值加速度為0.075g、0.100g、0.125g時,落沙量與發(fā)生頻率均服從冪律分布,但隨峰值加速度的減小,坍塌密度單調(diào)減小;當峰值加速度增加到0.15g～0.25g時,落沙量服從對數(shù)正態(tài)分布;當峰值加速度增加到0.35g～0.45g時,該階段樣本表現(xiàn)為具有正態(tài)分布的曲線特征.

振動臺沙堆模型實驗得到了與地震實震現(xiàn)象規(guī)律和元胞自動機模擬類似的結(jié)論,更進一步,雖然目前尚未獲得Ⅻ度區(qū)的實震資料,但不妨預測應具有向正態(tài)分布發(fā)展的趨勢.

振動臺沙堆模型實驗的物理過程更接近地震觸發(fā)崩塌滑坡的原型問題,雖然實驗次數(shù)有限,但仍然對作者提出的不同烈度區(qū)地震觸發(fā)崩塌滑坡分布規(guī)律的普適性提供了有力支持.

表3 振動臺沙堆模型實驗結(jié)果Tab.3 Results of the shaking table test of sandpile model

3 應用案例

具有SOC性質(zhì)的系統(tǒng),在臨界狀態(tài)下受到一系列微小而均勻的擾動,每次擾動下表征反應規(guī)模的物理量可用冪律描述,這種冪律關(guān)系將災害學的研究與物理學聯(lián)系在一起.規(guī)模與頻率的冪次定律描述的是某種量級災害發(fā)生的數(shù)目,這種災害演變總體規(guī)律在災勢分析、風險評估和危險性區(qū)化中的應用,是SOC較為明確的應用領(lǐng)域,對此舉一例進行說明.

除此之外,還將介紹兩個屬于應用技巧層面的案例,旨在說明SOC存在更廣闊的應用前景.

3.1 巖土體地震擾動深度評估

汶川大地震后,在震區(qū)調(diào)查時,觀察到在無深大結(jié)構(gòu)面控制的情況下,道路邊坡地震失穩(wěn)一般表現(xiàn)為淺表的崩塌滑坡,具有“剝皮型”災害特點.與日本OYO公司合作,采用面波儀對國道213線都江堰至映秀段(含水磨支線)震后邊坡進行了勘察,發(fā)現(xiàn)自然坡體0～5 m范圍內(nèi)的表層剪切波速為250 m/s,明顯低于下部土體的波速,判定為汶川地震的強擾動范圍[18].但僅以這些表觀認識和有限的測量數(shù)據(jù)為依據(jù),提出具有普遍性和包容性的結(jié)論是困難的.為此對照1∶2 000的地形圖,對沿線公路邊坡進行詳查,內(nèi)容涵蓋沿線公路邊坡崩塌滑坡的規(guī)模(崩塌滑坡方量)、崩塌滑坡體深度、工點經(jīng)緯度、高程、巖性以及原有工程措施.災害類型主要包括砂巖邊坡滑坡,巖質(zhì)邊坡崩塌,崩坡積體、殘坡積體、沖洪積體等淺層崩塌滑坡.

調(diào)查中發(fā)現(xiàn)不同崩塌體的崩塌深度存在很大區(qū)別.在Ⅸ度裂度區(qū)調(diào)查的61個崩塌體中,崩塌深度小于0.5 m的27個、0.5～1 m的15個、1～2 m的11個、2～5 m的4個、5～10 m的3個、10 m以上的1個,平均崩塌深度為1.2 m.

令崩塌滑坡體深度為h,深度大于h的工點數(shù)為N(h),通過回歸分析,擬合曲線如圖15所示.由圖15可以看出,h與N(h)在雙對數(shù)坐標圖上近乎一直線,關(guān)系式為

N(h)=20.27h-0.994 5.

(6)

由式(6)可知,相關(guān)系數(shù)R2=0.973,表明Ⅸ度裂度區(qū)崩塌滑坡深度與出現(xiàn)頻率之間存在良好的冪律關(guān)系,說明在Ⅸ度區(qū)地震觸發(fā)崩塌滑坡深度也存在SOC特性.這樣若仍以平均崩塌深度作為巖土體地震擾動深度的表征已不合理.

式(6)為Ⅸ度裂度區(qū)崩塌滑坡體深度h與大于h的工點數(shù)N(h)之間的冪律回歸關(guān)系,當樣本數(shù)量足夠大時,在統(tǒng)計意義上可得隨機變量h的分布函數(shù).如據(jù)實震調(diào)查資料可初步作出h的概率分布如圖16所示.

由圖16可知,汶川地震Ⅸ度裂度區(qū)93%的崩塌滑坡體深度均小于5 m,這是基于SOC理論開展的工作,其意義在于建立的是以物理理論為基礎(chǔ)的統(tǒng)計模型,跨越了僅據(jù)典型樣本得出統(tǒng)計關(guān)系的層面,其規(guī)律具有普適性,可以作為汶川地震Ⅸ度區(qū)地震對坡體強擾動效應的定量表征.圖16所示的分布規(guī)律可為汶川震后次生山地災害物源條件估算以及確定邊坡災害防范深度提供指導.

圖15 汶川地震Ⅸ度區(qū)崩塌滑坡深度累積分布Fig.15 Log-log plot of landslides number and depth in zones with seismic intensity Ⅸ

圖16 汶川地震Ⅸ度區(qū)崩塌滑坡體深度概率分布Fig.16 Probability distribution curve of landslide depth in the zone with s eismic intensity Ⅸ in Wenchuan earthquake

3.2 泥石流徑流量的極值計算

在工程規(guī)劃與設(shè)計中,常常需要對未來工程服務(wù)期內(nèi)情況進行預測,包括對外荷載最大值的預計.但人們往往僅掌握有限時段的觀察極值,如每年洪水的最大值.為實現(xiàn)從以前觀測到的極值數(shù)據(jù)運用外推法求得工程服務(wù)期內(nèi)的極值統(tǒng)計問題,在數(shù)學統(tǒng)計學中已建立了極值統(tǒng)計這一領(lǐng)域,重要成果之一是來自帶有指數(shù)型衰減尾部的初始分布的極值將漸近地收斂于Ⅰ型極限形式[19].在以下的研究中,基于泥石流徑流量的頻率-規(guī)?？捎脙缏擅枋鲫P(guān)系,通過數(shù)學上的映射,將其轉(zhuǎn)換成帶有指數(shù)型衰減尾部的分布,從而達到了運用Ⅰ型極限形式的條件,據(jù)此建立了不同設(shè)計基準期蔣家溝泥石流最大徑流量極值分布的數(shù)學模型.

一場泥石流的總徑流量是反映泥石流規(guī)模的主要表征值,也是開展泥石流攔擋工程容量設(shè)計、泥石流堵河可能性預測等工作的關(guān)鍵參數(shù).實際工程中,關(guān)注的是在工程設(shè)計基準期內(nèi)可能出現(xiàn)的泥石流徑流量的極大值,這些極大值可用具有一定概率分布的隨機變量來模擬,但在缺少長序列觀測資料的情況下,確定極值概型是極為困難的.

極值理論是處理一定樣本容量極端值分布特性的理論,其基本概念為:對{Xt,t∈T}隨機變量族,將時間域T等分為n個時段,每個時間段τ=T/n,則得n個獨立的隨機變量ηi(i=1,2,…,n),每個隨機變量都有相同的分布函數(shù)F(x).在這個隨機變量族中,極大值為

M=max(x1,x2,…,xn).

(7)

在上述兩個假定(獨立、具有相同的分布函數(shù))的前提下,極大值的分布函數(shù)為

FM(x)=p(M≤x)=

p(η1≤x,η2≤x,…,ηn≤x)=

(8)

式中:F(x)為初始分布.

若F(x)的密度函數(shù)為f(x),則有

fM(x)=n[F(x)]n-1f(x).

(9)

由于FM(x)與初始分布F(x)有著密切的關(guān)系,這樣在隨機變量的初始分布以及樣本容量n已知時,可寫出極值概率的確切分布.但在通常情況下,確切分布的求解十分困難,在工程實際應用中,常用做法就是使用極值的漸近分布.所謂漸近分布,就是當數(shù)目n趨于無窮大時,概率分布的極限形式,即漸近分布將收斂于確切分布,這樣可用漸近分布代替實際難以獲得的確切分布.

根據(jù)Gumbel的研究[20],極值分布的漸近形式并不取決于初始分布的精確形式,主要決定于初始分布在極值方向的尾數(shù)表現(xiàn),初始分布的中央部分對極值分布的漸近形式影響很小,但是,極值參數(shù)則決定于初始分布的形式.工程上常用的漸近分布有3種類型,其中最大值的Ⅰ型漸近分布多用于荷載的極值統(tǒng)計方面,如橋梁使用荷載、風速、地震震級、洪水流量等.

在泥石流研究領(lǐng)域,P A Johnson等研究表明[21],泥石流規(guī)模與頻率之間存在冪律關(guān)系,是自組織臨界狀態(tài)系統(tǒng)的行為標志.艾南山等認為[22]在泥石流的形成區(qū)松散堆積物組元間的非線性作用,使系統(tǒng)自然地朝著臨界狀態(tài)演化,在臨界狀態(tài)時,諸如降水等外界的細微擾動將被放大而導致規(guī)模不等的泥石流暴發(fā),這種耗散動力學系統(tǒng)的行為特征,可以用自組織臨界狀態(tài)的概念加以解釋.王裕宜[23]從泥石流土體的應力應變自組織臨界特性觸發(fā),認為粘性泥石流具有自組織臨界性.按照SOC系統(tǒng)的定義,泥石流徑流量規(guī)模W與頻率N之間應服從冪律分布,即N=aW-b.假設(shè)在所統(tǒng)計資料中,規(guī)模W的最大值和最小值分別為Wmax和Wmin.設(shè)X=lnW,則

N=a0e-bX,

(10)

式中:a0為常數(shù);b為冪指數(shù);X≥Xmin,Xmin=lnWmin,一般可取為0.

由頻度代替概率的思想,得出X的分布函數(shù)為

1-e-b(x-Xmin).

(11)

利用極值分布的最大吸引場的定理[24],可以判定上式極值分布的極限形式為收斂于Ⅰ型的漸近分布,進而推出X極大值的極限分布函數(shù)為

FM(x)=exp(-exp(-α(x-μ))),

(12)

相應的概率密度函數(shù)為

FM(x)=αe-α(x-μ)exp(-exp(-α(x-μ))),

(13)

α=c2/σQ1=1.28255/σQ1,

(14)

F1(x)=exp(-exp(-3.042×10-6×

(x-538 680))).

(15)

當設(shè)計基準期T=50 a時,已知年最大徑流量的分布,統(tǒng)計時段為1 a,τ=1.n=T/τ=50.根據(jù)極值I型分布的特征可知,設(shè)計基準期最大泥石流徑流量分布的標準差與初始分布相同

σQ50=σQ1=421 620 m3.

均值平移距離為

均值為

這樣可得設(shè)計基準期為50 a時,蔣家溝泥石流最大徑流量的概率分布為

F50(x)=exp(-exp(-3.042×10-6×

(x-1 824 684))).

(16)

同理可推出設(shè)計基準期為100、300 a時,蔣家溝泥石流最大徑流量的概率分布為

F100(x)=exp(-exp(-3.042×10-6×

(x-2 052 542))),

(17)

F300(x)=exp(-exp(-3.042×10-6×

(x-2 413 690))).

(18)

概率密度曲線如圖17所示.

圖17 蔣家溝不同設(shè)計基準期的徑流量極值概率密度Fig.17 Probability density curve of extreme runoff values in different design base periods of Jiangjia ravine

3.3 基于地震活動性參數(shù)b值的地應力評估模型

系統(tǒng)動力學服從冪律是人們對SOC系統(tǒng)的基本認識.但冪指數(shù)b不僅僅是一個統(tǒng)計分析參數(shù),還是反映系統(tǒng)動力學特征與活躍度的主要指標,不同系統(tǒng)b值的差異或者同一系統(tǒng)b值的演化,有著重要的物理意義.以下基于地震活動性參數(shù)b值與地應力的反相關(guān)關(guān)系,推求區(qū)域地應力的研究,雖然是一種統(tǒng)計層面的工作,但借此說明,從冪律到b值的討論,有可能開拓出SOC更廣闊的應用空間.

地應力是長大深埋隧道工程設(shè)計所需的重要參數(shù),現(xiàn)場鉆探是獲得準確初始地應力的可靠手段,但是在鐵路選線的可行性研究階段,由于線路方案尚未確定,不宜開展大規(guī)模的地應力現(xiàn)場鉆探工作,因此,地應力評估模型的發(fā)展就成為具有明確應用背景的課題.

目前,在地應力場分布特征的研究中多采用統(tǒng)計實測地應力、建立回歸模型的方法來分析主應力隨埋深的分布規(guī)律.主要代表性的工作有:中國科學院武漢巖土力學研究所對中國大陸地區(qū)實測地應力沿埋深分布規(guī)律做的相關(guān)研究[29];中國地震局地殼應力研究所根據(jù)統(tǒng)計的地應力實測數(shù)據(jù)對中國陸域八大地塊水平主應力隨埋深分布規(guī)律進行的回歸分析[30].但是上述研究工作提出的回歸模型,由于公式僅有一個自變參量,即埋深H,當模型的統(tǒng)計樣本未涵蓋評估區(qū)域時,就可能造成較大誤差.若能在回歸模型中引入能反映當?shù)氐貞顟B(tài)的參數(shù),則可望通過增加信息量的手段,提高計算精度.

3.3.1b值與地應力的關(guān)系

在前述的古登堡-里克特定律中,a反映平均地震活動水平,b反映大小地震的比例關(guān)系.隨后,在地震預報領(lǐng)域,分析發(fā)現(xiàn)，大震前震源及附近區(qū)域經(jīng)常會出現(xiàn)某些震級檔內(nèi)的地震增多或減小,導致出現(xiàn)大小地震比例失調(diào)、b值減小的異?，F(xiàn)象.由于區(qū)域應力積累水平升高是大地震發(fā)生的必要條件,因此,b值也反映了地應力狀態(tài),且二者呈反比關(guān)系.

之后有學者[31-33]根據(jù)單軸壓縮試驗過程中產(chǎn)生的聲發(fā)射(AE)現(xiàn)象對上述假說進行了驗證.巖石的聲發(fā)射現(xiàn)象描述為:巖石受力變形時,在巖體內(nèi)原先存在或新產(chǎn)生的裂縫周圍形成應力集中,這些局部應力集中區(qū)不均勻發(fā)展,并發(fā)生突然的破裂,從而向四周輻射彈性波.聲發(fā)射活動與地震活動的機制最為接近,在統(tǒng)計參數(shù)上與地震活動性的可對比性也最強,具體表現(xiàn)為:實驗記錄到的聲發(fā)射次數(shù)與聲壓幅值關(guān)系服從冪律分布,其中,冪律關(guān)系式中的冪指數(shù)b值不僅為一常數(shù),而且隨著應力的增加明顯下降,即巖石處于高應力狀態(tài)時,b值低.例如從中國國家地震局地球物理研究所公布的實驗資料(圖18)中可看出,巖石破裂實驗聲發(fā)射b值隨應力的增加而降低,在達到破裂應力之前b值下降較快[34].

基于b值與地應力相關(guān)的原理,由于b值是可以通過歷史地震資料統(tǒng)計獲得的參數(shù),因此,在水平主應力隨埋深分布規(guī)律統(tǒng)計回歸方程的基礎(chǔ)上,可以增加b值作為反映各地地應力狀態(tài)的信息參量,對現(xiàn)有地應力評估模型進行改進.

圖18 5種巖石聲發(fā)射b值隨應力變化圖[34]Fig. 18 Change with stress ofthe b-values of acoustic emission from five types of rock [34]

3.3.2 基于b值的川藏交通廊道地應力評估模型

川藏交通廊道穿越印度洋板塊俯沖歐亞板塊形成的青藏高原,板塊邊界的作用力是構(gòu)造變形的主要動力源,同時造成板塊內(nèi)部地震活躍、大地變形、高地應力等現(xiàn)象[35].如圖19所示,從青藏高原東緣進藏,需穿越三大縫合帶(金沙江縫合帶、班公怒江蛇綠巖縫合帶、雅魯藏布江蛇綠巖縫合帶)和七大斷裂帶(龍門山斷裂帶、鮮水河斷裂帶、甘孜-理塘斷裂帶、金沙江-紅河斷裂帶、甲?？?赤布張錯斷裂帶、納木錯-仲沙斷裂帶、嘉黎至然烏斷裂帶),構(gòu)造活動較強烈.此外,進藏線路靠近龍門山地震帶,穿越甘孜爐霍地震帶、雅魯藏布江地震帶,地震烈度均在Ⅷ度以上,地震活動強烈.

以北緯28°～32.1°、東經(jīng)90.1°～104.1°范圍內(nèi),1970年1月～2013年12月共45 a的地震資料作為基礎(chǔ)數(shù)據(jù),震中分布見圖20.將地震記錄以0.2震級檔為間隔進行分級整理.利用Arcgis軟件將研究區(qū)平面圖以0.1°×0.1°的間距進行網(wǎng)格化,挑選出以每個網(wǎng)格節(jié)點為圓心、半徑為r的圓形統(tǒng)計單元內(nèi)的地震資料[36],對于多數(shù)節(jié)點,統(tǒng)計單元的半徑r=30 km,但由于地震活動程度存在較大差異,對于地震分布較稀疏的局部區(qū)域,統(tǒng)計單元r值可適當擴大,以滿足統(tǒng)計所需的地震樣本數(shù).對于每一個統(tǒng)計單元,采用非線性擬合方法計算出擬合公式以及擬合優(yōu)度R2,確定能滿足整個研究時段的最小完整性震級Mc[37],不同統(tǒng)計單元的Mc是有差異的,計算中選擇擬合優(yōu)度R2取最大值時對應的震級為Mc.利用最佳擬合優(yōu)度法,得到最小完整性震級取Mc時的地震活動性參數(shù)b值,作為該單元中心點(即網(wǎng)格節(jié)點)的計算b值.利用Arcgis軟件將得到的整個研究區(qū)域掃描計算網(wǎng)格點的b值對離散點進行插值計算,經(jīng)數(shù)值范圍和區(qū)域顏色的調(diào)整后得到b值的空間分布等值線圖(圖20).由圖20可知,b值空間分布沿構(gòu)造斷裂帶存在著明顯的空間差異,反映出不同斷裂帶以及同一斷裂帶不同斷裂段現(xiàn)今應力積累水平的差異.川藏交通廊道中，雅安、八美、道孚、甘孜、巴塘至二九六工班、邦達至八宿、通麥、墨竹工卡、曲松等地為異常低b值區(qū),地應力積累水平較高.

圖19 川藏交通廊道線路、斷裂帶、縫合帶、中小地震震中分布圖Fig. 19 Distribution of lines, faults, sutures, and earthquake epicenters on thetransportation corridor from Sichuan to Tibet

圖20 川藏交通廊道b值空間分布Fig. 20 Distribution of b-values of the transportation corridor from Sichuan to Tibet

本文中收集了10個測區(qū),102組地應力原始實測數(shù)據(jù).所收集的原始地應力數(shù)據(jù)測量方法多為鉆孔應力解除法和水壓致裂法,少數(shù)為凱瑟爾效應法.地應力測試工程主要包括:南水北調(diào)西線一期工程、二郎山隧道、嘎隆拉隧道、大崗山水電站、兩河口水電站等.各測點所在區(qū)域的b值見表4.

利用上述資料進行統(tǒng)計回歸分析,建立了實測點最大主應力σ1、埋深H、地震b值三者的關(guān)系式,即川藏交通廊道地應力評估模型為

σ1=2.270 76-1.417b+0.500 4H.

(19)

利用建模時未使用的磽磧測點、康定測點和墨竹工卡測點的33組實測數(shù)據(jù),將式(19)與現(xiàn)有中國陸域八大地塊水平主應力評估模型[30]進行對比分析,有24組數(shù)據(jù)計算值的相對誤差比現(xiàn)有模型的相對誤差小，表明基于地震活動性參數(shù)b值的地應力評估模型精度明顯提高.

4 結(jié)束語

(1) 基于戴維斯地貌發(fā)育理論,按照山地系統(tǒng)的發(fā)育階段,界定了山地系統(tǒng)具備SOC的條件,并提出了基于斯特拉勒積分的山地系統(tǒng)發(fā)育階段判別方法;以地震觸發(fā)崩塌滑坡作為研究山地系統(tǒng)災變行為SOC的切入點,對不同地震烈度區(qū)地震觸發(fā)崩塌滑坡的分布規(guī)律開展了一項專題研究;介紹了3個SOC在山地系統(tǒng)的應用案例;形成了運用SOC認識和描述山地系統(tǒng)災變行為的理論架構(gòu).

(2) SOC是1987年作為非平衡態(tài)統(tǒng)計力學的一個分支建立起來的,對它的現(xiàn)象學研究和對它進行嚴謹?shù)亩x研究仍在進行.SOC已經(jīng)使人們意識到閾值、亞穩(wěn)定性、還有大規(guī)模波動在一大類多體系統(tǒng)的時空行為中起了決定性的作用,但擾動強度的變化對SOC系統(tǒng)動力學行為特性的影響是尚未有人關(guān)注的領(lǐng)域,而自然界,災變事件的擾動強度變化范圍可能會達到幾個數(shù)量級(如Ⅵ級地震到Ⅸ級地震能量相差達32 768倍).在以汶川地震、蘆山地震為原型的研究中,發(fā)現(xiàn)隨擾動強度增加, SOC系統(tǒng)的動力學特性將經(jīng)歷冪律—冪律弱化—對數(shù)正態(tài)分布的演變過程.這一有關(guān)SOC系統(tǒng)行為演化模式的認識,走出了SOC的傳統(tǒng)研究領(lǐng)域,對SOC基礎(chǔ)理論的發(fā)展也是有力的推動.

(3) 在SOC應用領(lǐng)域的研究尚屬初步的工作.通過對SOC系統(tǒng)演化行為模式的深化認識,并與山地災害實踐建立聯(lián)系,推動山地學從唯象學向精確科學跨越,無疑是具有深遠意義的工作.鑒此,本文希望通過范例性的研究產(chǎn)生拋磚引玉之效.

(4) 從科學的實證性來說,相對混沌理論而言, SOC不過分強調(diào)系統(tǒng)的演化依賴于系統(tǒng)的初始條件和細節(jié)部分,從而更容易分析和實驗,沙堆模型的概念似乎可以用來解釋從山脈形成到股市波動的幾乎所有的事.但若言山地系統(tǒng)的許多獨特的細節(jié)均可以借助簡單的元胞自動機數(shù)值模擬來理解,對于多數(shù)地理學家而言,可能是不大現(xiàn)實的,同樣也受到我們的質(zhì)疑.例如,自然界中的坡體結(jié)構(gòu)(節(jié)理、裂隙面等)都具有非均勻的分布,天然散粒體往往也具有寬級配的特征,因此,組元的非均勻性是山地系統(tǒng)的重要特征之一.但在目前的研究中,未考慮組元非均勻性對其動力學的影響,究其原因主要是尚未開發(fā)出能夠反映組元非均勻性的元胞自動機模型.因此,考慮通過元胞幾何尺寸的非均勻性、排列的隨機性、相互作用的各向異性,來構(gòu)建非均勻元胞自動機模型,探索非均勻性對系統(tǒng)動力學影響的獨特細節(jié),揭示山地系統(tǒng)在擾動作用下的深層次災變規(guī)律,將是本團隊下一步的主攻方向.

致謝:從散粒體到山地系統(tǒng)災變行為的SOC,本團隊已經(jīng)過了二十余年的研究歷程.先后有十多位教師和研究生參加過該領(lǐng)域的研究工作,對本文直接做出貢獻的有郭海強、段書蘇、郭沉穩(wěn)等博士研究生和韓俊碩士研究生等.在此一并深致謝意.

[1] BAK P, TANG C, WIESENFELD K. Self-organized criticality: an explanation of the 1/f noise[J]. Physical Review Letters, 1987, 59(4): 381-384.

[2] HELD G A. Experimental study of critical mass fluctuations in an evolving sand-pile[J]. Physical Review Letters, 1990, 65(9): 1120-1123.

[3] BAK P, TANG C. Earthquakes as a self-organized criticality phenomenon[J]. Geophysical Research, 1989, 97(B11): 15635-15637.

[4] 於崇文. 地質(zhì)系統(tǒng)的復雜性[M]. 北京:地質(zhì)出版社,2003: 93-140.

[5] 許強,黃潤秋. 巖石破裂過程的自組織臨界特征初探[J]. 地質(zhì)災害與環(huán)境保護,1996,7(1): 25-30.

XU Qiang, HUANG Runqiu. Disscussion on self-organized critical characters in the course of rock failure[J]. Journal of Geological Hazards and Enviroment Preservation, 1996, 7(1): 25-30.

[6] 姚令侃,李仕雄,蔣良濰. 自組織臨界性及其在散粒體研究中的應用[J]. 四川大學學報:工程科學版,2003,35(1): 8-14.

YAO Lingkan, LI Shixiong, JIANG Liangwei. Self-organized criticality and its application in granular mixtures[J]. Journal of Sichuan University: Engineering Science Edition, 2003, 35(1): 8-14.

[7] 姚令侃,黃藝丹,楊慶華. 地震觸發(fā)崩塌滑坡自組織臨界性研究[J]. 四川大學學報:工程科學版,2010,42(5): 33-43.

YAO Lingkan, HUANG Yidan, YANG Qinghua. The self-organized criticality of landslids triggered by earthquake[J]. Journal of Sichuan University: Engineering Science Edition, 2010, 42(5): 33-43.

[8] DAVIS W M. The geographical cycle[J]. The Geographical Journal, 1899, 14(5): 481-504.

[9] CHORLEY R J, SCHUMM S A, SUGDEN D E. Geomorphology[M]. London: Methuen & Co. Ltd., 1984: 58.

[10] STRAHLER A N. Hypsometric (area-altitude) analysis of erosional topography[J]. Geological Society of America Bulletin, 1952, 63(11): 1117-1142.

[11] 艾南山. 侵蝕流域系統(tǒng)的信息熵[J]. 水土保持學報,1987(2): 1-8.

AI Nanshan. Comentropy in erosion drainage-system[J]. Journal of Soil and Water Conservation, 1987(2): 1-8.

[12] 吳立新,史文中. 地理信息系統(tǒng)原理與算法[M]. 北京:科學出版社,2003: 441.

[13] 楊景春,李有利. 地貌學原理[M]. 北京:北京大學出版社,2005: 24-26.

[14] 中國科學院青藏高原綜合科學考察隊. 西藏地貌[M]. 北京:科學出版社,1983: 30-31,200-210.

[15] YAO Lingkan, FANG Duo. On the self-organized criticality of non-uniform sands[J]. International Journal of Sediment Research, 1998, 13(3): 19-24.

[16] BASTIEN C, MICHEL D. 物理系統(tǒng)的元胞自動機模擬[M]. 祝玉學,趙學龍，譯. 北京:清華大學出版社,2003: 1-2.

[17] 黃藝丹,姚令侃,郭沉穩(wěn). 基于元胞自動機的地震觸發(fā)崩塌滑坡分布規(guī)律[J]. 西南交通大學學報,2013,48(4): 609-615.

HUANG Yidan, YAO Lingkan, GUO Chenwen. Distribution law of landslides triggered by earthquake based on cellular automata[J]. Journal of Southwest Jiaotong University, 2013, 48(4): 609-615.

[18] 胡志旭,姚令侃,王建,等. 面波測試在汶川強震區(qū)土體損傷調(diào)查中的應用[J]. 防災減災工程學報,2010,30(4): 466-470.

HU Zhixu, YAO Lingkan, WANG Jian, et al. Application of surface wave testing in survey of mountain soil damage in Wenchuan strong earthquake area[J]. Journal of Disaster Prevention and Mitigation Engineering, 2010, 30(4): 466-470.

[19] ANG A H-S, TANG W H. 工程規(guī)劃與設(shè)計中的概率概念[M]. 孫芳垂,陳星燾,顧子聰,譯. 北京:冶金工業(yè)出版社,1991: 268-269.

[20] GUMBEL E J. The statistical theory of extreme values and some practical applications[J]. National Bureau of Standards, 1954, 118(1)： 33.

[21] JOHNSON P A, MCCUEN R H, HROMADKA T V. Magnitude and frequency of debris flows[J]. Journal of Hydrology, 1991, 123: 69-82.

[22] 羅德軍,艾南山, 李后強. 泥石流暴發(fā)的自組織臨界現(xiàn)象[J]. 山地研究,1995,13(4): 213-218.

LUO Dejun, AI Nanshan, LI Houqiang. Breakout of debris flow as a self-organized critical phenomenon[J]. Mountain Research, 1995,13(4): 213-218.

[23] 王裕宜,詹錢登,嚴璧玉,等. 泥石流體的流變特性與運移特征[M]. 長沙:湖南科學技術(shù)出版社,2014: 1-3.

[24] 史道濟. 實用極值統(tǒng)計方法[M]. 天津:天津科學技術(shù)技出版社,2006: 81-83.

[25] GUMBEL E.J. Statistics of extremes[M]. New York: Columbia University Press, 1958: 269-272.

[26] 張軍, 熊剛. 云南蔣家溝泥石流運動觀測資料集(1987～1994)[M]. 北京:科學出版社,1997: 1-3.

[27] 吳積善,康志成,田連權(quán),等. 云南蔣家溝泥石流觀測研究[M]. 北京:科學出版社,1990: 1-4.

[28] 康志成,崔鵬,韋方強,等. 東川泥石流觀測研究站觀測實驗資料集(1961～1984)[M]. 北京:科學出版社,2006: 1-3.

[29] 景鋒,盛謙,張勇惠,等.中國大陸淺層地殼實測地應力分布規(guī)律研究[J]. 巖石力學與工程學報,2007,26(10): 2056-2063.

JING Feng, SHENG Qian, ZHANG Yonghui, et al. Research on distribution rule of shallow crustal geostress in China mainland[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(10): 2056-2063.

[30] 楊樹新,姚瑞,崔效峰,等. 中國大陸與各活動地塊、南北地震帶實測應力特征分析[J]. 地球物理學報,2012,55(12): 4207-4217.

YANG Shuxin, YAO Rui, CUI Xiaofeng, et al. Analysis of the characteristics of measured stress in Chinese mainland and its active blocks and North-South seismic belt[J]. Chinese Journal of Geophysics, 2012, 55(12): 4207-4217.

[31] SCHOLZ C H. The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[J]. Bulletin of the Seismolo-gical Society of America, 1968, 58: 399-415.

[32] WIEMER S, WYSS M. Mapping the frequency-magnitude distribution in asperities: An improved technique to calculate recurrence times? [J]. Journal of Geophysical Research, 1997, 102(B7): 15115-15128.

[33] WYSS M, SCHORLEMMER D, WIEMER S. Mapping asperities by minima of local recurrence time: San Jacinto-Elsinore fault zones[J]. Journal of Geophysical Research, 2000, 105(B4): 7829-7844.

[34] 李紀漢.b值影響因素的初步研究[J]. 地震學刊,1987(2): 49-53.

LI Jihan. A preliminary study onb-value[J]. Jouranl of Seismology, 1987(2): 49-53.

[35] 姚令侃,邱燕玲,魏永幸. 青藏高原東緣進藏高等級道路面臨的挑戰(zhàn)[J]. 西南交通大學學報,2012,47(5): 719-734.

YAO Lingkan, QIU Yanling, WEI Yongxing. Challenges in construction of railway and highway from Sichuan to Tibet through eastern margin of Tibetan Plateau[J]. Journal of Southwest Jiaotong University, 2012, 47(5): 719-734.

[36] 易桂喜,聞學澤,范軍,等. 由地震活動參數(shù)分析安寧河-則木河斷裂帶的現(xiàn)今活動習性及地震危險性[J]. 地震學報,2004,26(3): 294-303.

YI Guixi, WEN Xueze, FAN Jun, et al. Assessing current faulting behaviors and seismic risk of the Anninghe-Zemuhe fault zone from seismicity parameters[J]. Acta Seismologica Sinica, 2004, 26(3): 294-303.

[37] 劉麗芳,李志海,蔣長勝. 云南地區(qū)地震目錄最小完整性震級研究[J]. 地震研究,2012,35(4): 491-499.

LIU Lifang, LI Zhihai, JIANG Changsheng. Research on minimum magnitude of completeness for earthquake catalogue in Yunnan region[J]. Journal of Seismological Research, 2012, 35(4): 491-499.

姚令侃(1953—),博士,1996年起至今任職于西南交通大學,現(xiàn)為土木工程學院教授,博士生導師,防災減災研究所所長.主要研究方向為鐵路公路工程災害防治及安全技術(shù).獲中國科學院科技進步二等獎(1995)、四川省科技進步二等獎(2000)、國家科技進步二等獎(2009)、第十屆詹天佑成就獎(2010)、四川省科學技術(shù)進步獎一等獎(2013).四川大學兼職教授、中國科學院成都山地災害與環(huán)境研究所所外專家、中國水土保持學會泥石流滑坡專業(yè)委員會副主任、鐵道部有突出貢獻的中青年專家、第三批、第九批四川省學術(shù)和技術(shù)帶頭人、享受國務(wù)院政府特殊津貼專家.

注:黃藝丹為本文的共同第一作者.

E-mail: yaolk@swjtu.edu.cn

黃藝丹(1982—),博士,2003年起至今任職于西南交通大學,現(xiàn)為土木工程學院講師.主要研究方向為鐵路公路工程災害防治及安全技術(shù).

E-mail: huangyidan@swjtu.edu.cn

(中文編輯:秦 瑜 英文編輯:蘭俊思)

Self-Organized Criticality of Mountain System Catastrophic Behaviors

YAOLingkan1,2,3,HUANGYidan1,2,3

(1. School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China; 2. MOE Key Laboratory of High-Speed Railway Engineering, Chengdu 610031, China; 3. Road and Railway Engineering Research Institute, Sichuan Key Laboratory of Seismic Engineering and Technology, Chengdu 610031, China)

In order to study the macro-dynamic behavior of mountain systems, a view point was proposed by Davis’s theory of the erosion cycle that valleys in late childhood or mountains in middle adulthood have intrinsic properties of self-organized criticality(SOC). Based on Strahler’s integral, a method for distinguishing landform development stages in a large basin is established. For landslides triggered by earthquake, the data of landslides collected in earthquake zones were analyzed, a cellular automata model was built, and shaking table experiments were performed. As a result, an evolution law of landslides distribution in different seismic intensity zones was revealed. Specifically, there exists a strong negative power-law relationship between the sizes and frequencies of landslides in zones with seismic intensity Ⅶ, Ⅷ and Ⅸ; the relationship becomes a weak power law in zones with seismic intensity Ⅹ, and changes into a lognormal distribution in zones with seismic intensity Ⅺ. Finally, applications of SOC to soil depth assessment under earthquake disturbance, extreme runoff calculation in debris flow control design, and crustal stress evaluation based on the seismic activity parameterb-value, are introduced through three case studies.

mountain system; self-organized criticality(SOC); landslides triggered by earthquake; cellular automata

2015-12-01

國家自然科學基金資助項目(41172321,41571004); 國家重點基礎(chǔ)研究發(fā)展計劃課題(2008CB425802)

姚令侃,黃藝丹. 山地系統(tǒng)災變行為自組織臨界性研究[J]. 西南交通大學學報,2016,51(2): 313-330.

0258-2724(2016)02-0313-18

10.3969/j.issn.0258-2724.2016.02.011

P642

A