王宗辰，史健宇，原野

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 預(yù)報點選取和后處理

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

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

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

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

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

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

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

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

圖2 印度洋區(qū)域海嘯預(yù)計傳播時間和最大波幅預(yù)報（以2004年蘇門答臘海嘯為例）

圖3 印度洋區(qū)域海嘯預(yù)計傳播時間和岸段波幅預(yù)報（以2004年蘇門答臘海嘯為例）

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

3.2 性能說明

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

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

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

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

表3 海嘯傳播至驗潮站的走時計算結(jié)果與觀測對比

續(xù)表3

圖4 海嘯傳播至驗潮站的走時計算結(jié)果與觀測對比

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

圖5 驗潮站海嘯最大波幅計算結(jié)果與觀測對比

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

4 結(jié)論

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

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

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