應(yīng)超，于普兵，穆錦斌，周華民

（浙江省水利河口研究院，浙江杭州310020）

東海海嘯反問題預(yù)報模式研究——以日本“3·11”海嘯為例

應(yīng)超，于普兵，穆錦斌，周華民

（浙江省水利河口研究院，浙江杭州310020）

摘要：采用COMCOT數(shù)值模式，將日本東海岸劃分為18個單位震源，建立了日本東海岸海嘯數(shù)據(jù)庫并用非負(fù)約束的最小二乘法建立反問題預(yù)報模式。將本模式應(yīng)用于日本“3·11”海嘯，計算所得的海嘯初始水位有10 m的抬高與3 m的降低，與前人研究結(jié)果基本一致，預(yù)報的浮標(biāo)水位與實測資料符合良好。對比浙江省近岸潮位站實測海嘯波高，預(yù)報值與實測值偏差較大。若采用反問題反演的震源，通過COMCOT非線性模式求解近岸水位，可以大幅提高預(yù)報精度。

關(guān)鍵詞：COMCOT；海嘯；反問題；最小二乘法

1　引言

海嘯是一種具有強(qiáng)大破壞力的海浪，具有傳播速度快、破壞力大等特點(diǎn)。準(zhǔn)確及時的海嘯預(yù)報是減小海嘯災(zāi)害的重要舉措。海嘯預(yù)報中，海嘯波的到達(dá)時間、波高與預(yù)報速度是需要重點(diǎn)關(guān)心的問題。

反問題預(yù)報方法是海嘯預(yù)報中的一種重要方法，是指用預(yù)先確定好的地震源參數(shù)把地震源劃分為不同的斷層，用不同斷層單位活動產(chǎn)生的水位過程建立數(shù)據(jù)庫，然后回歸分析記錄的海嘯波形確定斷層滑動分布和地震震級。國內(nèi)外很多學(xué)者對反問題預(yù)報模式進(jìn)行過研究，并取得了豐碩的成果。Titov[1]等利用斷層滑動分布通過反問題算法確定越洋海嘯實時模型的初始條件，Wei[2]等采用反問題方法通過預(yù)先計算好的潮位過程的數(shù)據(jù)庫確定遠(yuǎn)場海嘯波形預(yù)測海嘯。在Johnson[3]地震源參數(shù)基礎(chǔ)上，Wei[2]等用反問題方法準(zhǔn)確確定滑層分布，并準(zhǔn)確計算夏威夷各潮位站的波高。Yamazaki[4]等用此方法很好地模擬了1944年Tonakai和1994年Kuril地震海嘯通過太平洋在日本產(chǎn)生的海嘯過程。李林燕[5]等采用反問題方法對馬尼拉海溝建立海嘯數(shù)據(jù)庫并建立預(yù)報模式，但并未采用實測海嘯資料進(jìn)行驗證，僅采用正反問題的計算結(jié)果進(jìn)行相互驗證。以往的這些研究，均未對中國東海附近海域進(jìn)行研究，COMCOT模型及反問題方法在中國東海海域的適用性均需進(jìn)行驗證與探討，另外反問題方法在近岸預(yù)報的適用性問題也應(yīng)進(jìn)行進(jìn)一步分析。

北京時間2011年3月11日13時46分，日本東北部海域發(fā)生了9.0級特大地震，從而誘發(fā)了災(zāi)害性海嘯，造成了巨大的人員傷亡與財產(chǎn)損失。此次海嘯已經(jīng)對我國浙東沿海造成影響，地震發(fā)生后6—9 h，浙江沈家門、大陳、坎門、石坪、石浦、健跳、福建東山及廣東汕頭、汕尾等潮位站先后監(jiān)測到振幅為10—55 cm的海嘯波[6]。因此研究日本“3·11”海嘯反問題預(yù)報模式對我國東海海嘯預(yù)警機(jī)制的建立與完善具有很大的借鑒意義。

2　反問題預(yù)報方法

反問題預(yù)報方法如下：

（1）設(shè)置海嘯浮標(biāo)。海嘯浮標(biāo)一般設(shè)置在地形平坦水深較大的海域，距震源有一定距離但不宜太遠(yuǎn)，從而保證浮標(biāo)既能較快采集數(shù)據(jù)又不被海嘯破壞；

（2）建立海嘯數(shù)據(jù)庫。將關(guān)心區(qū)域內(nèi)，所有可能發(fā)生地震的俯沖帶均劃分為若干個小單元，稱為單位震源。計算每個單位震源發(fā)生單位錯動時所產(chǎn)生的海嘯數(shù)據(jù)，并存為海嘯數(shù)據(jù)庫；

（3）海嘯數(shù)據(jù)記錄與震源反演。當(dāng)?shù)卣鸢l(fā)生后，浮標(biāo)會監(jiān)測海嘯數(shù)據(jù)。當(dāng)首波數(shù)據(jù)監(jiān)測完成后，即可以利用海嘯數(shù)據(jù)庫與反問題算法，得到計算域內(nèi)任意位置的海嘯波振幅與到達(dá)時間，并反演震源信息。

3　海嘯數(shù)據(jù)庫建立

3.1COMCOT數(shù)值模式簡介

COMCOT（Cornell Multi-grid Coupled Tsunami Model）模式由美國Cornell大學(xué)開發(fā)，可以計算海嘯的整個生命周期，包括海嘯生成、傳播、爬高與淹沒的整個過程。根據(jù)研究尺度不同，模型可以靈活選擇球面與直角坐標(biāo)系下的線性與非線性淺水方程，所有方程均采用顯式蛙跳有限差分法進(jìn)行離散。具體方程與數(shù)植格式詳見文獻(xiàn)[7]。COMCOT數(shù)值模型是目前較為成熟的海嘯模型，已經(jīng)成功模擬了2004印尼蘇門答臘海嘯[8]等著名海嘯。

圖1　計算區(qū)域與水深條件圖

3.2網(wǎng)格布置與計算參數(shù)

模型共布設(shè)3層嵌套網(wǎng)格，計算范圍包括東海、臺灣海峽、南中國海、日本及部分西太平洋海域，網(wǎng)格大小為2 min，首層網(wǎng)格地形數(shù)據(jù)采用ETOPO1數(shù)據(jù)，第二、三層網(wǎng)格地形數(shù)據(jù)采用海圖與實測地形插值得到。計算范圍及地形見圖1，網(wǎng)格布置見表1。計算時間步長應(yīng)滿足CFL穩(wěn)定條件，本文取固定時間步長為2 s。建立海嘯數(shù)據(jù)庫是基于線性假定，所以僅計算首層網(wǎng)格。

3.3初始條件與邊界條件

COMCOT模式利用海床位移量來估算地震引起的初始水面高度，這種做法的前提是假定海水不可壓縮，海床為剛性介質(zhì)，水面變動與地震引起的地層錯動同時發(fā)生[9]。海嘯的初始條件的確定是通過輸入地震斷層參數(shù)，由彈性斷層模型計算得到。

本文采用Okada斷層模型，共需輸入9個參數(shù)，分別為：震中經(jīng)度、緯度、震源深度、斷裂長度、斷裂寬度、滑動量、走向角、傾角、滑移角。斷層模型示意圖見圖2。

邊界條件的設(shè)定中，水邊界設(shè)為開邊界；選擇線性淺水波方程時，海陸邊界設(shè)為垂直反射邊界，選擇非線性淺水波方程時，采用移動邊界方案。

3.4日本東海岸單位滑塊布置

圖2　斷層模型示意圖

NOAA的太平洋海洋環(huán)境實驗室（PMEL）將太平洋潛在地震海嘯由15個板塊804個海嘯源組成，每個單位海嘯源的大小均為100 km×50 km，并預(yù)先計算每個海嘯源在深海傳播的結(jié)果組成太平洋海嘯數(shù)據(jù)庫[10]。本文針對日本“3·11”海嘯進(jìn)行研究，所以僅選取震源附近的18個單位震源建立海嘯數(shù)據(jù)庫，單位震源布置見圖3，震源參數(shù)見表2。

圖3　單位震源分布圖

4　反問題模式驗證

4.1反問題模式的數(shù)學(xué)方法

假設(shè)震源區(qū)域被分割為NS個單位震源。第i個震源在地震發(fā)生后t時刻，在計算域內(nèi)任意一點(diǎn)(x,y )所產(chǎn)生的水面波動可記為Gi(x,y,t )，那么t時刻處的總水面波動可表示為：

表1　網(wǎng)格布置與模式設(shè)置

表2　單位震源參數(shù)表

式（1）中，ci表示第i個單位震源的權(quán)重系數(shù)，需要通過反問題算法來確定。

地震發(fā)生后，一旦海嘯波到達(dá)浮標(biāo)處，浮標(biāo)就會記錄水面波動。假設(shè)浮標(biāo)所處位置為(x0,y0) ，那么浮標(biāo)所記錄的數(shù)據(jù)集可以表示為Zk(x0,y0,tk) , { k=1,Nt}，式中，Nt代表浮標(biāo)記錄的數(shù)據(jù)總數(shù)。從而權(quán)重系數(shù)ci可由非負(fù)約束的最小二乘法求得，即式（2）：

式中，A為預(yù)先建立的海嘯數(shù)據(jù)庫，c為待求的權(quán)重系數(shù)，b為浮標(biāo)記錄數(shù)據(jù)集。

值得注意的是，浮標(biāo)記錄個數(shù)Nt需大于等于參于反演的單位震源個數(shù)Ns，這樣ci才有唯一解。

4.2日本“3·11”海嘯反問題模式驗證

“3·11”海嘯發(fā)生后，浮標(biāo)21401在地震發(fā)生后78 min時記錄了完整的海嘯首波波形。本文采用表2所示的單位震源建立海嘯數(shù)據(jù)庫，并根據(jù)21401浮標(biāo)前78 min的海嘯實測數(shù)據(jù)進(jìn)行反演，求解各單位震源的權(quán)重系數(shù)，計算結(jié)果表明主要發(fā)生錯動的單元為3b—7b，其余單元的權(quán)重系數(shù)均為0，單位震源權(quán)重系數(shù)分布見表3，并根據(jù)OKADA模型計算地震發(fā)生后的初始水位，見圖4?！?·11”地震海嘯發(fā)生后，Tang[11]等利用浮標(biāo)站壓力數(shù)據(jù)，反演出震源可由六塊沿著日本海溝的100 km×50 km的單位震源組合而成，斷層的詳細(xì)參數(shù)如表4所示，利用這些斷層參數(shù)與OKADA模型作為初始條件，在海嘯到岸前5 h計算出了美國海岸線30多個城市的海嘯波爬高和淹沒情況，模擬值與后來的觀測值較為符合。Wei[12]等也利用同樣的斷層設(shè)置，模擬計算了海嘯在日本近岸的爬高和淹沒，結(jié)果表明，模擬計算的茨城縣和青森縣的淹沒范圍準(zhǔn)確度達(dá)到85.5%。利用Tang等反演的震源計算地震發(fā)生后的初始水位，見圖5。對比初始水位圖可以看出，兩者初始水位分布形態(tài)較為相似，且量級上也非常接近，均反演出初始海面位移有10 m的抬升和3 m的下沉。

表3　單位震源權(quán)重系數(shù)表

為進(jìn)一步驗證模型，采用21401，21419與21413 這3個浮標(biāo)的實測數(shù)據(jù)對本文的預(yù)報結(jié)果進(jìn)行驗證（見圖6）。由浮標(biāo)水位驗證圖可以看出，反問題模式的預(yù)報結(jié)果在首波到達(dá)時間與首波波高上的計算上基本與實測值相符。雖然反問題模式在首波波谷處的計算上有所偏差，但是從預(yù)報的角度來說，本模式已經(jīng)可以在極短的時間內(nèi)回答預(yù)報所關(guān)心的問題。

5　反問題模式近岸預(yù)報準(zhǔn)確度分析

表4　Tang等反演的震源信息

圖4　本文反演震源所產(chǎn)生初始水位（單位：m）

圖5　Tang反演震源所產(chǎn)生初始水位（單位：m）

圖6　浮標(biāo)實測數(shù)據(jù)驗證圖

圖7　浙江沿海潮位站驗證站布置圖

表5　近岸潮位站驗證表

圖8　浙江沿海潮位站驗證圖

由于反問題模式的基本假定是海嘯波在大洋中為線性傳播，計算海嘯數(shù)據(jù)庫時也僅采用了線性淺水方程，所以對于近岸點(diǎn)直接采用式（1）進(jìn)行預(yù)報難免會存在誤差，可以考慮首先利用反問題法反演震源信息，再用COMCOT非線性方程求解近岸海嘯波傳播過程。本文采用前述兩種方法分別計算“3·11”海嘯期坎門、石浦、沈家門潮位站的海嘯參數(shù)，COMCOT非線模擬時網(wǎng)格布置同3.2節(jié)，三層網(wǎng)格均參于計算。各潮位站布置見圖7，各方法計算值與實測值的對比見表5與圖8。由圖與表可見，各方法所計算的首波波峰到達(dá)時間與實測資料相差均在20 min以內(nèi)，具有較高的精度。反問題法直接預(yù)報的首波波峰值，較實測值有明顯的偏小，這是因為，反問題模式并無考慮海嘯波的近海非線性效應(yīng)與海嘯波的爬高過程。若采用反問題法反演的震源，采用COMCOT非線性模式進(jìn)行計算，首波波峰有了明顯的增大，與實測值的誤差基本在20%以內(nèi)，僅沈家門站偏差較大，這可能是地形精度的影響所致。本節(jié)計算說明，在近岸淺水區(qū)域，反問題預(yù)報法僅能對海嘯波的到達(dá)時間作初步估計，若要取得較準(zhǔn)確的預(yù)報結(jié)果，需利用反問題法反演的震源，采用非線性模型進(jìn)行計算得到。

6　結(jié)論

本文采用COMCOT數(shù)值模式建立了日本東海岸的海嘯數(shù)據(jù)庫，并采用非負(fù)約束的最小二乘法構(gòu)建海嘯預(yù)報模式。應(yīng)用本模式對日本“3·11”海嘯進(jìn)行反演，計算結(jié)果表明，本模式可在2 min內(nèi)完成反演與預(yù)報，計算的初始水位與前人研究成果基本一致，計算的浮標(biāo)水位波動與實測資料較為符合。在近岸處，反問題模式計算的首波波峰到達(dá)時間仍有一定的精度，但在首波波峰值的計算上有較大偏差。若需獲得較準(zhǔn)確的近岸預(yù)報結(jié)果，應(yīng)采用反問題法反演的震源信息，應(yīng)用非線性方程進(jìn)行計算。

參考文獻(xiàn)

[1] Titov V V, Mofjeld H O, Gonzalez F I, et al. Offshore Forecasting of Alaskan Tsunamis in Hawaii[M]//Tsunami Research at the End of a Critical Decade. Netherlands: Springer, 2001: 75-90.

[2] Wei Y, Cheung K F, Curtis G D, et al. Inverse Algorithm for Tsunami Forecasts[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2003, 129(2): 60-69.

[3] Johnson J M. Heterogeneous Coupling Along Alaska-Aleutians as Inferred From Tsunami, Seismic, and Geodetic Inversions[J]. Advances in Geophysics, 1998, 39: 1-116.

[4] Yamazaki Y, Wei Y, Cheung K F, et al. Forecast of Tsunamis from the Japan-Kuril-Kamchatka Source Region[J]. Natural Hazards, 2006, 38(3): 411-435.

[5]李林燕,毛獻(xiàn)忠.南海海嘯反問題預(yù)報模式[J].水動力學(xué)研究與進(jìn)展:A輯, 2012, 27(1): 62-67.

[6]王培濤,于福江,趙聯(lián)大,等. 2011年3月11日日本地震海嘯越洋傳播及對中國影響的數(shù)值分析[J].地球物理學(xué)報, 2012, 55(9): 3088-3096.

[7] Liu P L F, Woo S B, Cho Y S. Computer Programs for Tsunami Propagation and Inundation[J]. 1998.

[8] Wang X M, Liu P L F. An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami[J]. Journal of Hydraulic Research, 2006, 44(2): 147-154.

[9]潘文亮,王盛安. COMCOT數(shù)值模式的介紹和應(yīng)用[J].海洋預(yù)報, 2009, 26(3): 45-52.

[10] Gica E, Spillane M C, Titov V V, et al. Development of the Forecast Propagation Database for NOAA's Short-Term Inundation Forecast for Tsunamis(SIFT)[R]. 2008.

[11] Tang L J, Titov V V, Bernard E N, et al. Direct Energy Estimation of the 2011 Japan Tsunami Using Deep-Ocean Pressure Measurements[J]. Journal of Geophysical Research: Oceans (1978-2012), 2012, 117(C8): C08008.

[12] Wei Y, Chamberlin C, Titov V V, et al. Modeling of the 2011 Japan Tsunami: Lessons for Near-Field Forecast[J]. Pure and Applied Geophysics, 2013, 170(6-8): 1309-1331.

Study on inversion forecasting model for East China Sea——A case study of Japan“3·11”tsunami

YING Chao，YU Pu-bing，MU Jin-bin，ZHOU Hua-min
（Zhejiang Inst．of Hydraulics＆Estuary，Hangzhou 310020 China）

Abstract：The Japan’s east coast was divided into 18 unit source to set up a tsunami database by using COMCOT numerical model. An inversion forecasting model was established by using least non-negative square method based on the database. The model was applied to 2011 Tohoku tsuami, the initial tsunami water level with 10m increase and 3m decrease calculated by the model was basically the same as previous research, the buoy level of prediction is in good agreement with measured data. Comparing with tsunami heights measured by tidal stations at coastal area of Zhejiang Province, the deviation of forecasted and measured value is large. But the prediction accuracy.can be greatly improved by solving COMCOT nonlinear equations with source parameters inversed by the forecasting model.

Key words：COMCOT；tsunami；inverse problem epicenter；least square method

作者簡介：應(yīng)超（1988-），男，工程師，碩士，主要從事海岸防災(zāi)減災(zāi)及海洋水環(huán)境研究。E-mail:xinqing928@126.com

基金項目：浙江省自然基金（LY13E090001）；浙江省創(chuàng)新人才培養(yǎng)項目（2012F20049）；浙江省公益基金（2014C33057）；浙江省公益技術(shù)研究社會發(fā)展項目（2014C33057）；浙江省科技計劃（2015F50064）

收稿日期：2014-07-31

DOI:10.11737/j.issn.1003-0239.2015.03.005

中圖分類號：P731.36

文獻(xiàn)標(biāo)識碼：A

文章編號：1003-0239（2015）03-0036-07