基于粒子系統(tǒng)的滑坡災(zāi)害過程模擬仿真方法

杜志強(qiáng)，李 靜

(1.武漢大學(xué) 測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079；2. 地球空間信息技術(shù)協(xié)同創(chuàng)新中心，湖北 武漢 430079；3.香港中文大學(xué) 太空與地球信息科學(xué)研究所，香港 999077)

0 引 言

極端天氣頻發(fā)，自然災(zāi)害多有發(fā)生，其危害面廣，破壞性大，是對人民生命財(cái)產(chǎn)安全的最大威脅和對社會經(jīng)濟(jì)發(fā)展的重大制約因素[1-2]?；聦儆谕话l(fā)性自然災(zāi)害，其危害不言而喻[3]。由于突發(fā)性自然災(zāi)害不可預(yù)測，爆發(fā)突然，其拍攝往往比較困難[4]，所以，突發(fā)性自然災(zāi)害模擬仿真尤為重要，其有助于增強(qiáng)人們對自然災(zāi)害的認(rèn)知，并幫助決策者采取正確的防災(zāi)救災(zāi)措施。

當(dāng)前滑坡模擬的方法種類很多，其中將滑坡運(yùn)動假想為單點(diǎn)運(yùn)動的集中質(zhì)量模型[5]和基于連續(xù)介質(zhì)的力學(xué)模型[6]是由Hunger[7]于1995年歸納的滑坡兩大類模擬仿真模型；Heim提出的“雪橇模型”是最著名的集中質(zhì)量模型，它是描述滑坡運(yùn)動的第一個(gè)模型，其將滑坡假想為一個(gè)單一的點(diǎn)，模型的優(yōu)點(diǎn)在于它的簡單明晰，雖然可計(jì)算出滑坡運(yùn)動的重力，但無法預(yù)測滑坡的滑動方向；基于連續(xù)介質(zhì)力學(xué)模型在滑坡滑動距離及相關(guān)參數(shù)預(yù)測方面有重要作用，Sassa[6]利用三維模型來模擬滑坡體的分布特征，在研究過程中將滑坡體視為像泥石流一樣的流體。理論上來講，雖然可以應(yīng)用模型定量地計(jì)算出滑坡的滑動體積、影響范圍，但由于模型自身的精確度和基本假定等因素影響，滑坡的實(shí)際運(yùn)動不可能完全符合模型的計(jì)算結(jié)果，另一方面，滑坡本身的誘發(fā)因素很多，在滑動過程中也有較強(qiáng)的隨機(jī)性和模糊性。而粒子系統(tǒng)是一種可以生成模糊及不確定對象的計(jì)算模型的方法[8]，到目前為止，其被認(rèn)為是用來模擬不規(guī)則物體運(yùn)動的最成功的圖形生成算法。

本文進(jìn)行滑坡模擬的主要思路是利用PhysX物理引擎中的粒子系統(tǒng)中的粒子來代表滑坡體中的泥石等成分，通過粒子在模擬災(zāi)害環(huán)境中的運(yùn)動，來模擬滑坡體在真實(shí)環(huán)境中的分布及運(yùn)動情況。

1 粒子系統(tǒng)基本原理

粒子系統(tǒng)采用大量具有生命周期和特征屬性的微小粒子作為基本單元來描述模糊對象的運(yùn)動特征[9]。在粒子系統(tǒng)中，每個(gè)粒子單元都具有形狀、大小、顏色、位置、速度大小和方向、生命周期等屬性，其中生命周期用于描述粒子在系統(tǒng)中的生存時(shí)間[10]，圖1為粒子在系統(tǒng)中隨時(shí)間變化的流程圖。

圖1 粒子生存周期Fig.1 Life span of particle

其將滑坡假想為一個(gè)單一的點(diǎn)，模型的優(yōu)點(diǎn)在于它的簡單明晰，雖然可計(jì)算出滑坡運(yùn)動的重力，但無法預(yù)測滑坡的滑動方向；基于連續(xù)介質(zhì)力學(xué)模型在滑坡滑動距離及相關(guān)參數(shù)預(yù)測方面有重要作用，Sassa[6]利用三維模型來模擬滑坡體的分布特征，在研究過程粒子在系統(tǒng)中的運(yùn)動以系統(tǒng)中的參數(shù)為基準(zhǔn)，并添加一些隨機(jī)變量進(jìn)行控制，使粒子的運(yùn)動變化帶有很強(qiáng)的隨機(jī)性，因此，可以很逼真地模擬出模糊對象的運(yùn)動特征。

粒子系統(tǒng)是由很多細(xì)小的粒子的集合共同表示一個(gè)模糊對象所構(gòu)成的[11]。用粒子系統(tǒng)來模擬滑坡動態(tài)過程的具體方法如下：

1）產(chǎn)生新的粒子加入系統(tǒng)中；

2）賦予每一新粒子以一定屬性；

3）刪除那些已經(jīng)超過其生命周期的粒子；

4）根據(jù)粒子的動態(tài)屬性對粒子進(jìn)行移動和變換；

5）繪制并顯示由有生命的粒子組成的圖像。

1.1 粒子屬性

粒子的運(yùn)動狀態(tài)源于其屬性的設(shè)置，以下為本文中所用到的粒子屬性[11]。

1）粒子初始位置

粒子的位置信息用三維坐標(biāo)（X,Y,Z）進(jìn)行描述。本文中粒子的初始位置設(shè)定為災(zāi)害發(fā)生的源點(diǎn)。

2）粒子初始速度

粒子系統(tǒng)生成時(shí)的形狀決定了粒子的運(yùn)動方向。如果為球形，粒子將從粒子系統(tǒng)的原點(diǎn)向外運(yùn)動，如果是圓形或者矩形，粒子將向XY平面的上方運(yùn)動，同時(shí)，可以調(diào)節(jié)粒子的發(fā)射角度來控制粒子的運(yùn)動方向。本文中粒子的運(yùn)動方向由滑坡該點(diǎn)處的坡向決定。

3）粒子生命周期

生命周期是用時(shí)間的方式描述了粒子的生命歷程，本文中通過控制時(shí)間間隔來控制粒子生命周期。

1.2 粒子動力學(xué)特征

粒子產(chǎn)生后，在三維空間中開始運(yùn)動，此時(shí)，粒子的屬性在運(yùn)動過程中會不斷改變。粒子的動力學(xué)特征是指每個(gè)粒子的屬性值在每一幀都要更新變化。

粒子在n+1幀的屬性由n幀的屬性決定，在第n+1幀時(shí)粒子的運(yùn)動位置[12]為：

式中，Xn+1表示在第f+1幀時(shí)粒子的X坐標(biāo)，VXn表示在第f幀時(shí)粒子的速度，△t表示第f幀到第f+1幀的時(shí)間差，其他同理。

第n+1幀時(shí)粒子的速度[12]為：

式中，Vn+1表示在第f+1幀時(shí)粒子的速度，a為粒子的加速度。

2 滑坡模擬模型設(shè)計(jì)

本文進(jìn)行滑坡模擬的主要思路是利用PhysX物理引擎中的粒子系統(tǒng)中的粒子來代表滑坡體中的泥石等成分，通過粒子在模擬災(zāi)害環(huán)境中的運(yùn)動，來模擬滑坡體在真實(shí)環(huán)境中的分布及運(yùn)動情況。

本文所用到的數(shù)據(jù)包括30 m全國地形數(shù)據(jù)、高分辨率影像數(shù)據(jù)（主要為高分?jǐn)?shù)據(jù)）以及間接觀測數(shù)據(jù)（如坡度、坡向、水平地震力等）。首先對高分辨率影像數(shù)據(jù)進(jìn)行影像糾正、信息提取等處理之后得到滑坡分析的宏觀信息。利用間接觀測數(shù)據(jù)對滑坡相關(guān)信息進(jìn)行參數(shù)率定，得到有效信息。根據(jù)極限平衡分析方法原理，生成滑坡的有效下滑數(shù)據(jù)。最后，利用有效下滑數(shù)據(jù)對PhysX物理引擎中的粒子系統(tǒng)進(jìn)行初始化，粒子發(fā)射器在滑坡范圍內(nèi)發(fā)射滑坡粒子，滑坡粒子將會受到環(huán)境影響而運(yùn)動。每一時(shí)刻，循環(huán)所有在生存周期內(nèi)的粒子，判斷滑坡粒子是否觸地，若不是，說明滑坡粒子還未落地，則該粒子暫時(shí)不處理，若是，則根據(jù)粒子位置判斷該粒子所在的格網(wǎng)，將對應(yīng)格網(wǎng)處數(shù)組加一，由此可得滑坡粒子數(shù)量數(shù)組，即落入每個(gè)格網(wǎng)中的粒子總數(shù)。每個(gè)格網(wǎng)的模擬時(shí)間即為該格網(wǎng)中的粒子數(shù)與時(shí)間間隔的乘積，若模擬時(shí)間大于0，則說明該格網(wǎng)處于受災(zāi)范圍區(qū)域，由此可得到滑坡災(zāi)害范圍。

本文進(jìn)行的滑坡模擬的流程圖如圖2所示。

圖2 滑坡模擬流程圖Fig.2 Process of landslide simulation

2.1 多參數(shù)率定方法原理

傳統(tǒng)的滑坡穩(wěn)定性分析方法，幾乎沒有考慮數(shù)據(jù)獲取過程中的人為誤差、系統(tǒng)誤差對數(shù)據(jù)的影響。如根據(jù)滑坡變形速率的分析方法，沒有經(jīng)過去差方法或只經(jīng)過簡單的線性擬合方法直接利用變形數(shù)據(jù)特點(diǎn)很容易引起滑坡預(yù)警的誤報(bào)錯(cuò)報(bào)[12]。如圖3紅色變形監(jiān)測區(qū)域，滑坡變形速率突然增大而后又出現(xiàn)減速現(xiàn)象。

圖3 傳統(tǒng)滑坡監(jiān)測方法容易引起誤報(bào)的監(jiān)測點(diǎn)Fig.3 Monitoring of false positives caused by traditional landslide monitoring method

滑坡穩(wěn)定性動態(tài)分析需要基于多種來源的觀測數(shù)據(jù)，包括地面監(jiān)測數(shù)據(jù)、高分辨率SAR數(shù)據(jù)、光學(xué)遙感數(shù)據(jù)、高時(shí)間、高空間分辨率數(shù)據(jù)、歷史災(zāi)情數(shù)據(jù)，在充分分析各種數(shù)據(jù)的特點(diǎn)和使用目的的基礎(chǔ)上利用尺度轉(zhuǎn)換技術(shù)、空間數(shù)據(jù)插值技術(shù)和屬性數(shù)據(jù)空間化技術(shù)實(shí)現(xiàn)多源數(shù)據(jù)的融合。在利用不同時(shí)刻的不能直接測定的土壤參數(shù)數(shù)據(jù)時(shí)，可利用線性回歸或多參數(shù)反演的方法來獲取相對時(shí)間的參數(shù)值，也可采用參數(shù)智能率定的方法，確保參數(shù)的實(shí)時(shí)準(zhǔn)確性。

2.2 極限平衡分析方法原理

極限平衡分析方法是邊坡穩(wěn)定分析中最常用的方法，它是通過分析在臨近破壞狀況下，土體外力與內(nèi)部強(qiáng)度所提供抗力之間的平衡，計(jì)算土體在自身和外荷作用下的土坡穩(wěn)定性程度，通常以邊坡穩(wěn)定系數(shù)表示[14]。

基于極限平衡分析方法[15]的計(jì)算公式（5），原理如圖4所示。

圖4 平衡分析方法原理示意圖Fig.4 Principle diagram of the limit equilibrium method

參數(shù)說明見表1。

表1 參數(shù)說明Tab.1 Parameter specification

2.3 滑坡粒子系統(tǒng)模型設(shè)計(jì)

根據(jù)粒子系統(tǒng)的基本原理和滑坡災(zāi)害體屬性分析，本小節(jié)將對滑坡災(zāi)害體的粒子系統(tǒng)進(jìn)行詳細(xì)設(shè)計(jì)?；铝Ｗ佑闪Ｗ影l(fā)射器發(fā)射，故粒子發(fā)射器的設(shè)計(jì)是極其重要的，其用于產(chǎn)生新粒子，并初始化粒子屬性。由于本文進(jìn)行的是滑坡的有源模擬，故可選用點(diǎn)發(fā)射源?；铝Ｗ拥奶匦园ɑ履M參數(shù)（包括模擬間隔、模擬時(shí)間等）、滑坡DEM數(shù)據(jù)、坡向數(shù)據(jù)、滑坡靜止密度、粘滯系數(shù)，粒子點(diǎn)特性（包括大小、顏色、速度等）以及滑坡邊坡穩(wěn)定系數(shù)等粒子屬性中的參數(shù)可直接或間接地從極限平衡分析算法中得到。

3 滑坡模擬仿真實(shí)驗(yàn)及結(jié)果分析

2015年11月13日22時(shí)50分許，浙江省麗水市蓮都區(qū)雅溪鎮(zhèn)里東村發(fā)生山體滑坡。山體滑坡塌方量30余萬立方米，導(dǎo)致27戶房屋被埋，21戶房屋進(jìn)水。截至2015年11月19日12時(shí)39分，現(xiàn)場找到第39位失聯(lián)人員。此次山體滑坡已有38人遇難，1人經(jīng)搶救后生命體征平穩(wěn)。

圖5為麗水滑坡實(shí)地影像圖：

圖5 麗水滑坡實(shí)地影像圖Fig.5 The image of landslide in Lishui

利用本軟件對從雅溪鎮(zhèn)里東村處發(fā)生的滑坡災(zāi)害進(jìn)行模擬，滑坡源頭設(shè)置在北緯28°38′49.2″，東經(jīng)119°54′36″，模擬開始時(shí)間為2015年11月13日22點(diǎn)50分，持續(xù)時(shí)間0.5 h，模擬時(shí)間間隔設(shè)置為5 s。圖6為模擬結(jié)果。

圖6 滑坡模擬結(jié)果圖Fig.6 The result of the landslide simulation

圖7 實(shí)際災(zāi)害與模擬災(zāi)害疊加圖Fig.7 The overlay chart of the actual disaster with the simulation result

由于模擬采用的DEM數(shù)據(jù)精度以及時(shí)效性等原因?qū)е聻?zāi)害模擬的范圍與實(shí)際災(zāi)害范圍基本吻合，但在細(xì)節(jié)上存在差異。另外，由于模擬采用的是有源模擬，故源點(diǎn)位置的選擇對模擬結(jié)果的影響至關(guān)重要。本例中選擇的泥石流的位置位于里東村。從圖7中看出模擬結(jié)果與實(shí)際災(zāi)害范圍基本吻合。

4 結(jié)束語

在地質(zhì)及物理學(xué)界，已有不少對山體滑坡等自然災(zāi)害現(xiàn)象的研究工作，但它們關(guān)注的重點(diǎn)往往是這些災(zāi)害現(xiàn)象的成因、運(yùn)動機(jī)制及數(shù)值模擬等。這些計(jì)算模型由于太過復(fù)雜并不適合于真實(shí)感模擬。本文在傳統(tǒng)模擬方法的基礎(chǔ)上引入了粒子系統(tǒng)，對滑坡過程進(jìn)行了動態(tài)模擬。模擬結(jié)果表明，本文方法能較好地模擬滑坡發(fā)生過程，展現(xiàn)滑坡危害區(qū)域，提供可靠防災(zāi)減災(zāi)資料。本文的模擬結(jié)果與實(shí)際災(zāi)情范圍具有較強(qiáng)的一致性，但由于DEM精度以及源點(diǎn)位置等因素，導(dǎo)致模擬結(jié)果出現(xiàn)部分偏差。