王宗辰,史健宇,原野
(國(guó)家海洋環(huán)境預(yù)報(bào)中心(自然資源部海嘯預(yù)警中心)自然資源部海洋災(zāi)害預(yù)報(bào)技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京100081)
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)。
蘇門答臘俯沖帶位于印度洋-澳大利亞板塊與巽他次級(jí)板塊的交匯部位,印度洋-澳大利亞板塊在這里向巽他次級(jí)板塊下面俯沖,是世界上最活躍的板塊構(gòu)造邊緣。自2004年以來(lái),該區(qū)域Mw8.0級(jí)以上地震引發(fā)的海嘯事件達(dá)到7 次,8.6 級(jí)以上地震海嘯3 次,其中最嚴(yán)重的2004 年蘇門答臘島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 年蘇門答臘地震海嘯及之后的3次事件對(duì)數(shù)值預(yù)報(bào)結(jié)果進(jìn)行釋用和性能分析。
海嘯數(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é)果。
實(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]。
海嘯波在大洋中傳播時(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ì)算效率
當(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)城市和港口。
選取2004—2010年發(fā)生在印度洋的4次Mw7.8級(jí)以上地震海嘯事件對(duì)數(shù)值預(yù)報(bào)結(jié)果和性能進(jìn)行說(shuō)明,震源參數(shù)見(jiàn)表2。引發(fā)2004 年蘇門答臘海嘯的地震震級(jí)超過(guò)了Mw9.0,且地震破裂過(guò)程非常復(fù)雜[2-3]。因此,在計(jì)算海嘯源時(shí),采用了Grilli 等[18]的多點(diǎn)源策略。
表2 海嘯事件震源參數(shù)
以2004年蘇門答臘大海嘯為例,解釋說(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è)主要傳播方向,包括震中附近的蘇門答臘島、泰國(guó)西海岸、緬甸南部、孟加拉灣、斯里蘭卡、印度東南海岸、馬爾代夫、索馬里、馬達(dá)加斯加?xùn)|南島弧以及印度洋中脊西南支。
圖2 印度洋區(qū)域海嘯預(yù)計(jì)傳播時(shí)間和最大波幅預(yù)報(bào)(以2004年蘇門答臘海嘯為例)
圖3 印度洋區(qū)域海嘯預(yù)計(jì)傳播時(shí)間和岸段波幅預(yù)報(bào)(以2004年蘇門答臘海嘯為例)
從圖3 我們能夠看到,震源附近的蘇門答臘第一時(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]。
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)性。以蘇門答臘地震海嘯為例,泰國(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 年的南蘇門答臘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ā)出。以蘇門答臘地震海嘯為例,除了位于震源周圍100 km 之內(nèi)的蘇門答臘島和安達(dá)曼尼科巴島之外,其他國(guó)家和相關(guān)機(jī)構(gòu)仍有1.5~2.5 h的時(shí)間做出減災(zāi)部署。
本文依托自然資源部海嘯預(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ā)生在蘇門答臘島周邊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)性參考。