印度洋區(qū)域地震海嘯數(shù)值預(yù)報(bào)系統(tǒng)的建立

王宗辰，史健宇，原野

（國(guó)家海洋環(huán)境預(yù)報(bào)中心（自然資源部海嘯預(yù)警中心）自然資源部海洋災(zāi)害預(yù)報(bào)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100081）

1 引言

2017 年6 月20 日，中華人民共和國(guó)國(guó)家發(fā)展和改革委員會(huì)與原國(guó)家海洋局聯(lián)合發(fā)布了《“一帶一路”建設(shè)海上合作設(shè)想》，勾勒出3條藍(lán)色經(jīng)濟(jì)通道，其中1 條便是：以中國(guó)沿海經(jīng)濟(jì)帶為支撐，連接中國(guó)-中南半島經(jīng)濟(jì)走廊，經(jīng)南海向西進(jìn)入印度洋，銜接中巴、孟中印緬經(jīng)濟(jì)走廊，直至非洲的“海上絲綢之路”。

從地質(zhì)構(gòu)造的角度來(lái)看，海上絲綢之路依次經(jīng)過(guò)了歐亞板塊、印度洋板塊和非洲板塊。板塊交界處常伴有較高的地震活動(dòng)性，因此也是海嘯頻發(fā)之地。統(tǒng)計(jì)顯示，自1604年以來(lái)該區(qū)域發(fā)生海嘯事件117 次[1]，其中南中國(guó)海及周邊海域79 次、印度洋區(qū)域38次（見(jiàn)圖1）。

蘇門(mén)答臘俯沖帶位于印度洋-澳大利亞板塊與巽他次級(jí)板塊的交匯部位，印度洋-澳大利亞板塊在這里向巽他次級(jí)板塊下面俯沖，是世界上最活躍的板塊構(gòu)造邊緣。自2004年以來(lái)，該區(qū)域Mw8.0級(jí)以上地震引發(fā)的海嘯事件達(dá)到7 次，8.6 級(jí)以上地震海嘯3 次，其中最嚴(yán)重的2004 年蘇門(mén)答臘島Mw9.2 級(jí)地震海嘯共造成24～28萬(wàn)人傷亡和失蹤，以及十億美元數(shù)量級(jí)的經(jīng)濟(jì)損失，成為有記錄以來(lái)死亡人數(shù)最多的自然災(zāi)害之一[2-4]。海嘯的破壞形式不僅包括近岸海水驟升造成的漫灘、漫堤和潰堤淹沒(méi)，還包括海嘯波到達(dá)近岸和港口引發(fā)的強(qiáng)流對(duì)基礎(chǔ)設(shè)施和船只的沖擊和破壞[5-6]。

自主開(kāi)展海上絲綢之路地震海嘯監(jiān)測(cè)和預(yù)警分析工作迫在眉睫。實(shí)際上，在太平洋區(qū)域海嘯預(yù)警與減災(zāi)系統(tǒng)框架下，南中國(guó)海區(qū)域已經(jīng)擁有獨(dú)立的海嘯預(yù)警系統(tǒng)[7]（責(zé)任范圍見(jiàn)圖1），由自然資源部海嘯預(yù)警中心負(fù)責(zé)運(yùn)維。以聯(lián)合國(guó)教科文組織海洋學(xué)委員會(huì)南中國(guó)海區(qū)域海嘯預(yù)警中心（South China Sea Tsunami Advisory Center/United Nations Educational, Scientific and Cultural Organization-Intergovernmental Oceanographic Commission，SCSTAC/UNESCO-IOC）的名義向南中國(guó)海周邊9 個(gè)國(guó)家提供海嘯預(yù)警服務(wù)。

本文將基于圖形處理單元（Graphic Processing Unit，GPU）并行加速技術(shù)建立印度洋海嘯數(shù)值預(yù)報(bào)系統(tǒng)。首先，收集統(tǒng)計(jì)了印度洋和南中國(guó)海區(qū)域內(nèi)的歷史地震海嘯事件，然后利用線性淺水方程建立了覆蓋印度洋的海嘯傳播計(jì)算模型；進(jìn)一步沿環(huán)印度洋域內(nèi)國(guó)家的海岸線選取了540個(gè)預(yù)報(bào)點(diǎn)對(duì)岸段危險(xiǎn)進(jìn)行警示，對(duì)重點(diǎn)城市和港口給出海嘯首波預(yù)計(jì)抵達(dá)時(shí)間和最大波幅；最后，利用2004 年蘇門(mén)答臘地震海嘯及之后的3次事件對(duì)數(shù)值預(yù)報(bào)結(jié)果進(jìn)行釋用和性能分析。

2 預(yù)報(bào)系統(tǒng)建立

海嘯數(shù)值預(yù)報(bào)系統(tǒng)包括海嘯源計(jì)算、海嘯傳播計(jì)算和海嘯預(yù)警分析3個(gè)模塊。海嘯源計(jì)算是為了獲取初始海表面位錯(cuò)分布，作為海嘯初始場(chǎng)驅(qū)動(dòng)海嘯傳播方程；繼而借助GPU 加速技術(shù)進(jìn)行海嘯傳播計(jì)算，產(chǎn)出海嘯傳播時(shí)間和最大波幅場(chǎng)；最后，預(yù)警分析模塊通過(guò)對(duì)波幅的后處理獲得岸段危險(xiǎn)性分析結(jié)果。

2.1 地震監(jiān)測(cè)和海嘯源計(jì)算

實(shí)際業(yè)務(wù)中，海嘯數(shù)值預(yù)報(bào)系統(tǒng)的啟用通常需要兩個(gè)條件：一是發(fā)生在海底的淺源地震且震級(jí)達(dá)到Mw7.1 級(jí)以上；二是依賴震源特征參數(shù)用于快速計(jì)算海嘯源。

條件一要求對(duì)地震事件具備實(shí)時(shí)監(jiān)測(cè)能力，能夠在地震發(fā)生后迅速給出地震發(fā)生的時(shí)間、位置、震級(jí)和深度。自然資源部海嘯預(yù)警中心開(kāi)展了全球地震監(jiān)測(cè)系統(tǒng)建設(shè)，現(xiàn)已具備全球地震快速監(jiān)測(cè)能力[8-9]，本文采用其速報(bào)參數(shù)。

條件二要求能夠快速估計(jì)震源的破裂長(zhǎng)度、寬度和滑移量，以及反演地震的震源機(jī)制。本文選擇了Blaser等[10]推導(dǎo)的大洋俯沖帶逆沖型定標(biāo)律來(lái)估算震源的破裂尺度；然后利用地震矩公式計(jì)算均一滑移量[11]；再采用W-phase 方法準(zhǔn)實(shí)時(shí)反演震源機(jī)制解[12-15]。

滿足上述兩個(gè)條件之后，即可利用OKADA 模型計(jì)算海底位錯(cuò)，即初始海表面位移作為海嘯源[16]。

2.2 數(shù)值預(yù)報(bào)模型

海嘯波在大洋中傳播時(shí)，波幅遠(yuǎn)小于水深，波長(zhǎng)遠(yuǎn)大于水深。因此，可以忽略非線性效應(yīng)，線性長(zhǎng)波方程便可很好地模擬海嘯波，但是應(yīng)考慮科氏力的影響。球坐標(biāo)系下的方程表達(dá)式為：

式中：η為相對(duì)于平均海平面的自由表面位移；φ為緯度；ψ為經(jīng)度；R為地球半徑；P為沿緯度單位寬度的通量；Q為沿經(jīng)度單位寬度的通量；f為科氏力系數(shù)；g為重力加速度。該模型對(duì)海嘯波在大洋的傳播模擬已經(jīng)經(jīng)過(guò)驗(yàn)證[11]。

模型計(jì)算區(qū)域覆蓋整個(gè)印度洋和南中國(guó)海以及蘇祿海和蘇拉威西海，具體范圍40°S ～28 °N、15°～130°E，地形水深采用GEBCO_14 Grid。線性方程的差分方案采用球坐標(biāo)系下蛙跳格式有限中心差分，輻射邊界條件，具體配置詳見(jiàn)表1。硬件計(jì)算環(huán)境為2 個(gè)14 核Intel(R) Xeon(R) CPU E5-2680 v4 芯片，主 頻為2.4 GHz，128 G 內(nèi)存；GPU 采 用Nvidia Tesla K40c 顯卡，計(jì)算耗時(shí)見(jiàn)表1。對(duì)于震級(jí)較小的計(jì)算過(guò)程應(yīng)采用2 arc-min 空間分辨率，時(shí)間成本變成約4 min。

表1 印度洋海嘯預(yù)報(bào)模型配置和計(jì)算效率

2.3 預(yù)報(bào)點(diǎn)選取和后處理

當(dāng)海嘯波傳播到水深較淺的近岸區(qū)域，非線性效應(yīng)變得顯著，線性方程不再適用。如果采用非線性模型直接計(jì)算沿岸波幅，則需要高精度地形資料配合以及更短的時(shí)間步長(zhǎng)，時(shí)間成本將大幅提高。為了兼顧岸段預(yù)報(bào)的準(zhǔn)確性和時(shí)效性，本預(yù)報(bào)系統(tǒng)采用格林公式[17]，將低分辨率模型輸出的離岸點(diǎn)（Offshore）海嘯波幅換算到水深極淺的沿岸預(yù)報(bào)點(diǎn)（Coast）。

式中：Ac和Ao分別為近岸和離岸預(yù)報(bào)點(diǎn)最大海嘯波幅；Ho和Hc分別對(duì)應(yīng)離岸和近岸點(diǎn)水深。根據(jù)模型分辨率的不同，Wang 等[17]對(duì)離岸點(diǎn)的水深做了條件限制。對(duì)于4 arc-min 的模型分辨率，Ho應(yīng)大于992 m；Hc使用固定值1 m。并且，離岸點(diǎn)和預(yù)報(bào)點(diǎn)之間的距離在不超過(guò)100 km 的條件應(yīng)該盡量取小。

按照約60 km 的間隔距離沿著北印度區(qū)域內(nèi)的澳大利亞、印度尼西亞和泰國(guó)等27個(gè)國(guó)家或者主權(quán)歸屬地選取了540 對(duì)岸段預(yù)報(bào)點(diǎn)（見(jiàn)圖3）。離岸點(diǎn)和近岸點(diǎn)平均間隔約50 km，超過(guò)95%的離岸點(diǎn)水深在（100 m，1 000 m）的區(qū)間，平均水深450 m，根據(jù)經(jīng)驗(yàn)，將近岸點(diǎn)水深調(diào)整為5 m（變形格林公式）。這540 個(gè)預(yù)報(bào)點(diǎn)中包括了環(huán)印度洋國(guó)家的78 個(gè)重點(diǎn)城市和港口。

3 預(yù)報(bào)結(jié)果和系統(tǒng)性能

選取2004—2010年發(fā)生在印度洋的4次Mw7.8級(jí)以上地震海嘯事件對(duì)數(shù)值預(yù)報(bào)結(jié)果和性能進(jìn)行說(shuō)明，震源參數(shù)見(jiàn)表2。引發(fā)2004 年蘇門(mén)答臘海嘯的地震震級(jí)超過(guò)了Mw9.0，且地震破裂過(guò)程非常復(fù)雜[2-3]。因此，在計(jì)算海嘯源時(shí)，采用了Grilli 等[18]的多點(diǎn)源策略。

表2 海嘯事件震源參數(shù)

3.1 預(yù)報(bào)結(jié)果釋用

以2004年蘇門(mén)答臘大海嘯為例，解釋說(shuō)明預(yù)報(bào)系統(tǒng)產(chǎn)出的定量海嘯預(yù)報(bào)結(jié)果，包括大洋最大波幅圖和岸段危險(xiǎn)警示圖（見(jiàn)圖2、3），以及重點(diǎn)城市、港口和驗(yàn)潮站的海嘯預(yù)計(jì)抵達(dá)時(shí)間。圖片右側(cè)標(biāo)注了地震發(fā)生的時(shí)間、位置、震級(jí)以及震源深度和快速反演的震源機(jī)制解。從圖2中可以看到海嘯能量的幾個(gè)主要傳播方向，包括震中附近的蘇門(mén)答臘島、泰國(guó)西海岸、緬甸南部、孟加拉灣、斯里蘭卡、印度東南海岸、馬爾代夫、索馬里、馬達(dá)加斯加?xùn)|南島弧以及印度洋中脊西南支。

圖2 印度洋區(qū)域海嘯預(yù)計(jì)傳播時(shí)間和最大波幅預(yù)報(bào)（以2004年蘇門(mén)答臘海嘯為例）

圖3 印度洋區(qū)域海嘯預(yù)計(jì)傳播時(shí)間和岸段波幅預(yù)報(bào)（以2004年蘇門(mén)答臘海嘯為例）

從圖3 我們能夠看到，震源附近的蘇門(mén)答臘第一時(shí)間受到海嘯沖擊，尤其位于其西北的亞齊省以及臨近的韋島產(chǎn)生了超過(guò)3 m 的海嘯波幅，最大波幅超過(guò)10 m。這意味著在沒(méi)有巨型防波堤的情況下，上述區(qū)域?qū)a(chǎn)生大面積海嘯入侵和淹沒(méi)。除此之外，泰國(guó)西海岸、馬來(lái)西亞半島西北海岸、緬甸南部西岸、安達(dá)曼尼科巴島、斯里蘭卡東岸、印度東南海岸以及馬爾代夫群島東側(cè)都也都出現(xiàn)了超過(guò)3 m的海嘯波，大面積淹沒(méi)不可避免，預(yù)報(bào)結(jié)果與災(zāi)后調(diào)查結(jié)果基本一致[19]。

3.2 性能說(shuō)明

W-phase 方法只能快速反演得到基于點(diǎn)源的震源機(jī)制，如果一個(gè)地震的破裂屬性存在顯著的空間變化，那么單點(diǎn)源震源機(jī)制刻畫(huà)的海嘯源代表性則較差。震源附近的海嘯傳播時(shí)間和波幅計(jì)算結(jié)果往往會(huì)產(chǎn)生較大偏差。除此之外，如果驗(yàn)潮站位于有遮擋的海灣、受防波堤保護(hù)的港口、或者驗(yàn)潮站離岸地形水深條件復(fù)雜，以及格林公式換算路徑上存在很寬的陸架或陡峭的陸坡，也都會(huì)對(duì)海嘯預(yù)報(bào)結(jié)果造成偏差[17]。一方面原因是模型分辨率不足，另一方面是不滿足格林公式的適用條件[11]。

基于以上原因，并非所有的驗(yàn)潮站均可用于驗(yàn)證本文的預(yù)報(bào)系統(tǒng)，或者說(shuō)，驗(yàn)潮站的觀測(cè)結(jié)果不完全能夠指示岸段的危險(xiǎn)性。以蘇門(mén)答臘地震海嘯為例，泰國(guó)的RANONG、KURABURI、PHUKET、KRABI 和TRANG 5 個(gè)驗(yàn)潮站記錄的海嘯波高普遍約為0.7 m，最大為1.1 m；而災(zāi)害調(diào)查結(jié)果顯示，上述驗(yàn)潮站所在岸段波高均值超過(guò)5 m；斯里蘭卡的COLOMBO 觀測(cè)到1.5 m 的海嘯波，災(zāi)后調(diào)查顯示其所在岸段海嘯波高普遍超過(guò)4 m。上述原因可能是驗(yàn)潮站在記錄到最大海嘯波高前就已經(jīng)失去功能。

因此，在選擇驗(yàn)潮站比對(duì)預(yù)報(bào)結(jié)果時(shí)，我們將不滿足格林公式適用條件的驗(yàn)潮進(jìn)行了剔除。4 次海嘯事件共篩選了44個(gè)驗(yàn)潮站次的數(shù)據(jù)作為觀測(cè)，其中2007 年的南蘇門(mén)答臘Mw8.4 地震海嘯事件缺少海嘯傳播時(shí)間數(shù)據(jù)，導(dǎo)致只有32 個(gè)站次數(shù)據(jù)可用。下面選擇海嘯傳播時(shí)間和海嘯最大波幅兩個(gè)指標(biāo)進(jìn)行分析。

計(jì)算的海嘯傳播時(shí)間與觀測(cè)相比，平均絕對(duì)誤差約為15 min，最大絕對(duì)誤差約為40 min（見(jiàn)圖4，具體的計(jì)算和觀測(cè)結(jié)果對(duì)比見(jiàn)表3）。對(duì)于國(guó)家尺度的海嘯預(yù)警，我們認(rèn)為這樣的誤差量級(jí)基本可以接受。更準(zhǔn)確的海嘯傳播時(shí)間依賴精細(xì)的海嘯源刻畫(huà)和高精度地形水深資料。

表3 海嘯傳播至驗(yàn)潮站的走時(shí)計(jì)算結(jié)果與觀測(cè)對(duì)比

續(xù)表3

圖4 海嘯傳播至驗(yàn)潮站的走時(shí)計(jì)算結(jié)果與觀測(cè)對(duì)比

按照國(guó)際慣用海嘯危險(xiǎn)等級(jí)和王宗辰等[11]對(duì)波幅預(yù)警的評(píng)價(jià)標(biāo)準(zhǔn)對(duì)計(jì)算結(jié)果進(jìn)行了分析（見(jiàn)圖5）。44 個(gè)站次的海嘯波幅絕對(duì)計(jì)算誤差約為17 cm，相對(duì)誤差為27%；其中危險(xiǎn)性合理判斷比例達(dá)到86%，高估和低估比例分別為9%和5%，對(duì)于大尺度的國(guó)家預(yù)警而言，這樣的波幅預(yù)報(bào)水平完全可以接受。如果再考慮圖3 和災(zāi)害調(diào)查結(jié)果的比較，系統(tǒng)性能分還可提高。

圖5 驗(yàn)潮站海嘯最大波幅計(jì)算結(jié)果與觀測(cè)對(duì)比

海嘯傳播計(jì)算用時(shí)約40 s，GMT 繪制圖件用時(shí)約6 s，整個(gè)數(shù)值預(yù)報(bào)流程耗時(shí)可以控制在50 s 之內(nèi)。如果在震后0.5 h可獲得震源機(jī)制解，數(shù)值預(yù)報(bào)的流程耗時(shí)基本可以忽略。依賴網(wǎng)站、郵件和短信等通訊手段，海嘯數(shù)值預(yù)報(bào)結(jié)果和危險(xiǎn)性信息可以在1 min 之內(nèi)發(fā)出。以蘇門(mén)答臘地震海嘯為例，除了位于震源周?chē)?00 km 之內(nèi)的蘇門(mén)答臘島和安達(dá)曼尼科巴島之外，其他國(guó)家和相關(guān)機(jī)構(gòu)仍有1.5～2.5 h的時(shí)間做出減災(zāi)部署。

4 結(jié)論

本文依托自然資源部海嘯預(yù)警中心的全球地震監(jiān)測(cè)和反演系統(tǒng)實(shí)時(shí)獲取地震速報(bào)和震源機(jī)制解，結(jié)合地震震級(jí)-破裂尺度定標(biāo)律和OKADA 公式計(jì)算海嘯源；再利用基于GPU 加速的線性海嘯傳播計(jì)算模型建立了印度洋區(qū)域的海嘯數(shù)值預(yù)報(bào)系統(tǒng)；沿著印度洋周邊國(guó)家的海岸線，按照50～60 km 的標(biāo)準(zhǔn)選擇了540 對(duì)離岸-近岸預(yù)報(bào)點(diǎn)，利用變形格林公式對(duì)岸段海嘯波幅進(jìn)行定量預(yù)報(bào)。

選取2004 年以來(lái)發(fā)生在蘇門(mén)答臘島周邊4 次Mw8.0級(jí)以上地震海嘯事件對(duì)預(yù)報(bào)系統(tǒng)的定量預(yù)警產(chǎn)品進(jìn)行釋用，并對(duì)系統(tǒng)性能進(jìn)行了測(cè)試。結(jié)果顯示，本系統(tǒng)能夠在得到震源機(jī)制解之后1 min 內(nèi)產(chǎn)出整個(gè)印度洋區(qū)域海嘯預(yù)警分析結(jié)果。4 次海嘯事件的后報(bào)結(jié)果顯示，海嘯傳播時(shí)間的絕對(duì)計(jì)算誤差約為15 min，危險(xiǎn)性合理判斷比例達(dá)到86%。

印度洋地震海嘯數(shù)值預(yù)報(bào)系統(tǒng)的建立不僅能夠?yàn)閲?guó)內(nèi)“海絲路”相關(guān)單位和機(jī)構(gòu)提供準(zhǔn)實(shí)時(shí)的定量海嘯預(yù)警報(bào)服務(wù)，借助網(wǎng)站和社交媒體，該系統(tǒng)還可以為印度洋周邊國(guó)家提供海嘯危險(xiǎn)性參考。