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

王宗辰，史健宇，原野

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

2.2 數值預報模型

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

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

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

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

2.3 預報點選取和后處理

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

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

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

3 預報結果和系統(tǒng)性能

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

表2 海嘯事件震源參數

3.1 預報結果釋用

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

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

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

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

3.2 性能說明

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

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

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

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

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

續(xù)表3

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

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

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

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

4 結論

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

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

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