東北大學(xué) 2017 年大學(xué)生數(shù)學(xué)建模競賽 A 題一等獎成果。
在 GitHub 上：https://github.com/lalxyy/NEU-2017-MCM-A
由于不能顯示行內(nèi)公式，部分 LaTeX 公式以原形式展現(xiàn)。
摘要
東京時間 2011 年 3 ? 11 ? 14 時 46 分，?本東北地?太平洋海域發(fā)??? 9.0 級地震，并伴隨有 7.0 到 7.4 級的多次強余震。9.0 級地震也造成了歷史記錄以來罕見的?達 39 ?的巨?海嘯。截?到 2011 年 5 ?，超過 2.4 萬?被官?報告為死亡或失蹤；另外，福島第?核電站因此發(fā)?的嚴(yán)重核事故也導(dǎo)致了千余?的?直接死亡，這些?并發(fā)?的重?事故也給?本及周邊國家?guī)砹藝?yán)重的損失。本?通過評估磐城、仙臺、宮古與橫須賀的地震破壞及損失情況，利?多元線性回歸的?式，建?了災(zāi)害有關(guān)因素（地震強度、海嘯等）與災(zāi)害造成的后果（?員傷亡、經(jīng)濟損失）的模型。
本?根據(jù)互聯(lián)?中能夠查得的數(shù)據(jù)進?提取，并獲得了沿岸海浪最??度、與震中距離、地震??烈度、經(jīng)濟?平（GDP）等數(shù)據(jù)，以及?員傷亡、財產(chǎn)損失的數(shù)據(jù)。根據(jù)繪制的散點圖確定多元擬合中每個分量對應(yīng)的多項式次數(shù)，之后將?次單項式換元轉(zhuǎn)化為?次單項式，使?多元線性擬合來得出多項式系數(shù)。利?所得出的模型，本?又簡要分析了仙臺市的災(zāi)情情況，并且指出數(shù)據(jù)點較少帶來的模型的不?之處。
關(guān)鍵詞： 多項式擬合; 地震破壞性評估; 東?本?地震
1 問題背景與分析
1.1 背景
2011 年 3 ? 11 ??本福島發(fā)??? 9.0 級?地震，強震引發(fā)海嘯，導(dǎo)致當(dāng)時全球最?的在役核電站——福島核電站放射性物質(zhì)外泄。該事故根據(jù)國際核事件分級被評為最?級 7 級，與切爾諾貝利核電站泄漏事故等級相同。截? 2014 年 2 ?，“3.11”東?本?地震和?海嘯死亡的?數(shù)達到 15884 ?。 截? 2014 年 5 ?，受東?本?地震和福島第?核電站事故影響，?直接死亡的?數(shù)已經(jīng)達到 1699 ?，超過了在事故中直接死亡的 1603 ?。
建?數(shù)學(xué)模型，在以下?本沿海城市：磐城、仙臺、宮古和橫須賀 中。?較各種規(guī)模??的地震破壞以及所造成的后果，并為當(dāng)?shù)貓蠹垳?zhǔn)備?篇?章，解釋根據(jù)你的模型關(guān)于其中?個城市的發(fā)現(xiàn)。
1.2 問題分析
該問題本質(zhì)為求解?個地震因?向量 x = (x_1 , x_2 , ..., x_m) 到影響程度向量 y = (y_1 , y_2 , ..., y_n) 的函數(shù) F（可能包含多個線性變換），表示為
，怎樣得出 F(x) 也是本?討論的核?內(nèi)容。
由于災(zāi)害產(chǎn)?的影響與?類?產(chǎn)??平的發(fā)展情況和???的主觀能動性有很?的關(guān)聯(lián)，所以單純憑借歷史以來的地震傷亡?數(shù)、財產(chǎn)損失情況不能作為衡量當(dāng)今災(zāi)害情況與財產(chǎn)損失的直接依據(jù)，但是?類??因素又是衡量地震帶來的破壞性的必要考慮因素，且關(guān)系極其復(fù)雜。為了簡化模型、提?模型準(zhǔn)確性與說服?，我們將磐城、仙臺、宮古和橫須賀四個地區(qū)在同?時間橫向?qū)?，不需控制時間的變量，?是通過與震中的距離來衡量地震的影響規(guī)模。
題?所給出的四個地區(qū)的建設(shè)、經(jīng)濟發(fā)展?平不盡相同。圖 1 是 Google Maps 的數(shù)據(jù)，給出了四城的地理位置與震中坐標(biāo)，最右側(cè)海域中的是本次地震的震中，其余點從上到下依次為宮古、仙臺、磐城、橫須賀；
以及四城 2010 年當(dāng)時的國民?產(chǎn)總值（GDP）[1][2][3][4]：
| 城市 | 國民生產(chǎn)總值 |
|---|---|
| 仙臺 | 61.7 |
| 橫須賀 | 12.2882 |
| 磐城 | 16.5905 |
| 宮古 | 5.19 |
表 1: 2010 年四城國民?產(chǎn)總值（單位：億美元）
由此我們將經(jīng)濟因素（代表抗災(zāi)、救災(zāi)?平）納?到考慮因素中。除此之外，根據(jù)常識，地震破壞性本?又與震中位置、震源深度以及城市與震中的距離相關(guān)，也是模型中必要的考慮因素。Google Maps 分別提供了四城市中?到震中的直線距離（表 2）；
| 城市 | 與震中距離 / km |
|---|---|
| 宮古 | 186.45 |
| 仙臺 | 176.75 |
| 磐城 | 207.31 |
| 橫須賀 | 423.28 |
表 2: 震中直線距離
此外，連帶因素還包括由于?本部分區(qū)域近海?產(chǎn)?的由地震引發(fā)的海嘯影響、以及福島核電站放射性物質(zhì)外泄造成的事故。四個城市的海嘯浪?度最?值可以在?象數(shù)據(jù)中查得 [5][6]（表 3）； 以及從各城市與所屬縣政府查得的?員傷亡與財產(chǎn)損失情況（表 4）[7][8][9][10]。
| 城市 | 海嘯沿岸最??度 / ? |
|---|---|
| 宮古 | 8.5 |
| 仙臺 | 8.6 |
| 磐城 | 6.51 |
| 橫須賀 | 2.06 |
表 3: 觀測到的海嘯?度
注：部分觀測點不以城市命名，選取對應(yīng)最近的觀測點。
| 城市 | 傷亡?數(shù) | 經(jīng)濟損失（億美元） |
|---|---|---|
| 宮古 | 611 | 22.4775 |
| 仙臺 | 904 | 119.0022 |
| 橫須賀 | 0 | 0.92 |
| 磐城 | 466 | 0.1478 |
表 4: 四城傷亡?數(shù)與經(jīng)濟損失
注：橫須賀市政府沒有說明有任何?員傷亡情況（參考?獻 10），擬合過程中暫計為 0。（附錄 B.2）
2 模型假設(shè)與說明
2.1 符號說明
| 符號 | 代表意義與說明 |
|---|---|
| A_1 | 仙臺市 |
| A_2 | 橫須賀市 |
| A_3 | 宮古市 |
| A_4 | 磐城市 |
| X | 表示地震整體因素的向量 |
| x_1, ..., x_m | 影響地震破壞性的各個因素 |
| Y | 表示地震破壞性結(jié)果的向量 |
| y_1, ..., y_n | 地震破壞程度的不同衡量特征（財產(chǎn)損失等） |
| F | 最終求解的地震破壞性衡量函數(shù) |
| d_1, ..., d_4 | 對應(yīng) A_1 , ..., A_4 到震中的直線距離，單位為千? |
表 5: 本?所使?的符號說明
2.2 模型假設(shè)
假設(shè) 1：?本本?海岸線在地震中沒有?規(guī)模（超過 10 ?的顯著）偏移；[11]
假設(shè) 2：?類的能動性（主觀性與?產(chǎn)?）與發(fā)展時間相關(guān)且滿?映射關(guān)系；
假設(shè) 3：?類在地震中的抗災(zāi)防災(zāi)?平與經(jīng)濟發(fā)展程度正相關(guān)；
假設(shè) 4：?個城市區(qū)域內(nèi)的建筑物強度可以??個平均值來表示，且可以?建筑物的強度、破壞程度作為衡量城市破壞程度的主要因素。
3 模型建立
3.1 知識背景
3.1.1 地震影響場
采?的地震烈度影響場衰減公式：
沿長軸?向，
沿短軸?向，
其中 d 表示城市與震中的直線距離。[12]
3.1.2 海嘯等級
采?太平洋地區(qū)普遍適?的今村－飯?強度分級來評定海嘯的強度等級，其中 H_av 表示沿最近海岸的平均海嘯?度，I 為輸出強度值。[13]
3.2 具體問題參數(shù)引入
由問題分析所述，已經(jīng)基本可以確定 X、Y 各分量的表?含義。將 x_1 , x_2 , y_1 , ... 等分量指派為表 6、7 中所示的含義；四城中每個城市都與且僅與?組向量對應(yīng)。
| 維度 | 表示含義與單位說明 |
|---|---|
| x_1 | 地震震級（??） |
| x_2 | 國民?產(chǎn)總值（億美元） |
| x_3 | 輸?值在 d_1 , ..., d_4 范圍內(nèi)，與震中直線距離（千?） |
| x_4 | 海嘯沿岸最??度（?） |
表 6: 地震影響因?向量 x 各分量表示含義
| 維度 | 表示含義與單位說明 |
|---|---|
| y_1 | 傷亡?數(shù) |
| y_2 | 經(jīng)濟損失（億美元） |
表 7: 影響結(jié)果 y 各分量表示含義
并通過 Python 語?與 Matplotlib 庫繪制 y 各分量與 x 各分量散點圖：
3.3 基本模型
根據(jù)現(xiàn)有散點圖及統(tǒng)計回歸模型，認(rèn)為 y_1 , y_2 為被解釋變量，x_1 , ..., x_4 為解釋變量。觀察散點圖，認(rèn)為 y_1 與 x_3 ，y_2 與 x_3 滿??次多項式擬合的圖像特征，且函數(shù)對應(yīng)單調(diào)遞減：
其中 a_{ij}, b_{ij} 為需要求出的多項式系數(shù)，i = 1, 2, 3, 4，j = 0, 1, 2，而且擬合曲線應(yīng)當(dāng)單調(diào)遞減。類似方法可得 y_1 與 x_2、y_2 與 x_2 線性關(guān)系，
以及 y_1 與 x_4, y_2 與 x_4 的圖像特征也可用二次多項式來擬合：
y_1 與 y_2 的取值受到以上因素的影響。所以可以得出關(guān)于 x_2, x_3, x_4 的兩個三元多項式擬合函數(shù)
，其中 c_1，c_2 為該多項式擬合模型中的常數(shù)，也就是
4 模型求解
引?擬合?法中的最??乘法，使?多元線性回歸的?法求解該模型。令
，此時原函數(shù)變?yōu)槲逶€性擬合。
對于任意一個取值點 (x_{ij}, y_{kj})，其中 j = 0, ..., 3，i = 1, ..., 6，k = 1, 2（已知四個城市對應(yīng)的數(shù)據(jù)點），有
，對于 y_1，也就是
此時令
，就有
兩側(cè)同時左乘 A^T，得
，整理得
，求出 c_1, a_{21}, a_{31}, a_{41}, a_{32}, a_{42}。
使用計算機求出的最后結(jié)果為
5 模型評價與改進
該模型可以根據(jù)已有的?然災(zāi)害數(shù)據(jù)較為準(zhǔn)確地求得財產(chǎn)損失與?員傷亡的程度，能夠給出范圍內(nèi)的合理值，可?于防災(zāi)減災(zāi)與災(zāi)害造成后果的快速估計與判斷。然?本模型利?的數(shù)據(jù)點較少，如果使??量的數(shù)據(jù)來進?擬合則會使模型更有說服?。
6 推廣與應(yīng)用：仙臺市
仙臺市是?本東北地?的區(qū)域中?，也是本題?已知條件中離震中最為接近的城市。從之前繪制的圖像和模型預(yù)估的情況來看，仙臺的情況并不容樂觀：常識上講，?個發(fā)達區(qū)域的經(jīng)濟發(fā)展?平越?，防災(zāi)減災(zāi)能?也應(yīng)該越強，這也是本?在模型建?時考慮 GDP（國民?產(chǎn)總值）作為防災(zāi)減災(zāi)因素的原因；但是由此看來，仙臺及?本東北地?城市仍然需要加強對于防災(zāi)減災(zāi)?采取的措施，例如削弱海嘯因素帶來的影響，可以體現(xiàn)在本模型中的擬合曲線更為平滑，?且海嘯是導(dǎo)致此次地震及連帶事故中?員傷亡與財產(chǎn)損失的?個主要的原因。
參考文獻
[1] 宮古市役所. 平成 23 年度 財政狀況資料集 [EB/OL]. 2011. http://www.city.miyako.iwate.jp/data/open/cnt/3/5677/1/H23zaiseijoukyou.pdf.
[2] 橫須賀市. 財政?書（平成 23 年度決算）[EB/OL]. 2012. http://www.city.yokosuka.kanagawa.jp/1610/finas/documents/hakusyo23.pdf.
[3] KANEMOTO Y. Metropolitan Employment Area (MEA) Data[EB/OL]. 2011. http://www.csis.u-tokyo.ac.jp/UEA/uea_data_e.htm.
[4] いわき市役所. 平成 23 年度決算 [EB/OL]. 2012. http://www.city.iwaki.lg.jp/www/contents/1001000003486/simple/H23_kessann.pdf.
[5] ?本気象庁. 國內(nèi)の津波観測施設(shè)で観測された津波の観測値 [EB/OL]. 2011. http://www.data.jma.go.jp/svd/eqev/data/2011_03_11_tohoku/tsunami_jp.pdf.
[6] ?本気象庁. 気象庁｜津波観測點（東北地?）[EB/OL]. 2017. http://www.data.jma.go.jp/svd/eqev/data/tsunamimap/.
[7] 宮古市役所. The Great East Japan Earthquake and Tsunami Records of Miyako City[EB/OL]. 2015. http://www.city.miyako.iwate.jp/data/open/cnt/3/5971/1/records_of_miyako_city. pdf.
[8] 仙臺市役所. 東?本?震災(zāi)における本市の被害狀況等 [EB/OL]. 2017. https://www.city.sendai.jp/okyutaisaku/shise/daishinsai/higai.html.
[9] いわき市役所. いわき市災(zāi)害対策本部週報 [EB/OL]. 2015. http://www.city.iwaki.lg.jp/www/contents/1449132951986/simple/higai0524.pdf.
[10] 橫須賀市. 東?本?震災(zāi)関連情報 [EB/OL]. 2011. http://www.city.yokosuka.kanagawa.jp/shinsai311/index.html.
[11] NASA. Japan’s Coastline Before and After the Tsunami[EB/OL]. 2011. https://www.nasa.gov/multimedia/imagegallery/image_feature_1893.html.
[12] 傅再揚, 危福泉, 林巖釗. 福建地震災(zāi)害損失計算分析與思考 [J]. 災(zāi)害學(xué), 2007, 22(1) : 90 – 93.
[13] PAPADOPOULOS G A, IMAMURA F. A proposal for a new tsunami intensity scale[C] // ITS 2001 proceedings. 2001 : 569 – 577.
附錄
A 程序源碼
帶有運?結(jié)果的源碼與部分互動界?可以通過 Jupyter Notebook 打開根?件夾中的 /j-note ?件 夾內(nèi)的 ipynb ?件來查看。
以下代碼均為 Python 語?（Python 3）。
import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model
import sklearn.preprocessing
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimSun'] # 在其他機器上可能需要修改字體名
mpl.rcParams['axes.unicode_minus'] = False
mpl.rcParams['savefig.dpi'] = 108
A.1 數(shù)據(jù)引入
# GDP、震中直線距離、沿岸海浪高度、震級
# x2, x3, x4, x1
miyako_x = np.array([5.19, 186.45, 8.5, 9.0]) # 宮古
iwaki_x = np.array([16.5905, 207.31, 6.51, 9.0]) # 磐城
yokosuka_x = np.array([12.2882, 423.28, 2.06, 9.0]) # 橫須賀
sendai_x = np.array([61.7, 176.75, 8.6, 9.0]) # 仙臺
matrix_x = np.array([miyako_x, iwaki_x, yokosuka_x, sendai_x])
# 人員傷亡、財產(chǎn)損失
miyako_y = np.array([611, 22,4775])
iwaki_y = np.array([466, 0.1478])
yokosuka_y = np.array([0, 0.92])
sendai_y = np.array([904, 119.9022])
matrix_y = np.array([miyako_y, iwaki_y, yokosuka_y, sendai_y])
A.2 散點圖繪制
def draw_scatter_pic(x, y, description, x_label, y_label):
"""
繪制散點圖的（統(tǒng)一方式）函數(shù)
"""
fig = plt.figure(dpi=600)
ax1 = fig.add_subplot(111)
ax1.set_title(description)
plt.xlabel(x_label)
plt.ylabel(y_label)
ax1.scatter(x, y, c = 'r', marker = 'o')
# plt.legend('x1')
plt.show()
# x_gdp = np.array([miyako_x[0], iwaki_x[0], yokosuka_x[0], sendai_x[0]])
x_gdp = matrix_x[:, 0:1]
y_person = np.array([miyako_y[0], iwaki_y[0], yokosuka_y[0], sendai_y[0]])
draw_scatter_pic(x_gdp, y_person, '散點圖：人員傷亡與 GDP', 'GDP', '人員傷亡')
x_distance = np.array([miyako_x[1], iwaki_x[1], yokosuka_x[1], sendai_x[1]])
# fig = plt.figure(dpi=600)
# ax1 = fig.add_subplot(111)
# ax1.set_title('散點圖：震中直線距離與傷亡人數(shù)')
# plt.xlabel('直線距離')
# plt.ylabel('人員傷亡')
# ax1.scatter(x_distance, y_person, c='r', marker='o')
# plt.show()
draw_scatter_pic(x_distance, y_person, '散點圖：震中直線距離與傷亡人數(shù)', '直線距離', '人員傷亡')
x_tsunami_height = np.array([miyako_x[2], iwaki_x[2], yokosuka_x[2], sendai_x[2]])
fig = plt.figure(dpi=600)
ax1 = fig.add_subplot(111)
ax1.set_title('散點圖：沿岸海浪高度與傷亡人數(shù)')
plt.xlabel('海浪高度')
plt.ylabel('人員傷亡')
ax1.scatter(x_tsunami_height, y_person, c='r', marker='o')
plt.show()
y_loss_economy = y_person = np.array([miyako_y[1], iwaki_y[1], yokosuka_y[1], sendai_y[1]])
draw_scatter_pic(x_gdp, y_loss_economy, '散點圖：GDP 與經(jīng)濟損失', 'GDP', '經(jīng)濟損失')
draw_scatter_pic(x_distance, y_loss_economy, '散點圖：震中直線距離與經(jīng)濟損失', '直線距離', '經(jīng)濟損失')
draw_scatter_pic(x_tsunami_height, y_loss_economy, '散點圖：沿岸海浪高度與經(jīng)濟損失', '海浪高度', '經(jīng)濟損失')
miyako_x_extended = np.array([1, miyako_x[0], miyako_x[1], miyako_x[2], miyako_x[1] * miyako_x[1], miyako_x[2] * miyako_x[2]])
iwaki_x_extended = np.array([1, iwaki_x[0], iwaki_x[1], iwaki_x[2], iwaki_x[1] * iwaki_x[1], iwaki_x[2] * iwaki_x[2]])
yokosuka_x_extended = np.array([1, yokosuka_x[0], yokosuka_x[1], yokosuka_x[2], yokosuka_x[1] * yokosuka_x[1], yokosuka_x[2] * yokosuka_x[2]])
sendai_x_extended = np.array([1, sendai_x[0], sendai_x[1], sendai_x[2], sendai_x[1] * sendai_x[1], sendai_x[2] * sendai_x[2]])
np_x = np.array([miyako_x_extended, iwaki_x_extended, yokosuka_x_extended, sendai_x_extended])
np_y = np.array([miyako_y[0], iwaki_y[0], yokosuka_y[0], sendai_y[0]])
A = np.matmul(np.linalg.inv(np.matmul(np.transpose(np_x), np_x)), np.matmul(np.transpose(np_x), np_y))
print(A)
np_y = np.array([miyako_y[1], iwaki_y[1], yokosuka_y[1], sendai_y[1]])
B = np.matmul(np.linalg.inv(np.matmul(np.transpose(np_x), np_x)), np.matmul(np.transpose(np_x), np_y))
print(B)
B 關(guān)于數(shù)據(jù)及其來源的說明
本?使?的所有數(shù)據(jù)均來源于各城市政府官?或其他學(xué)者分析報告，在參考?獻中都已經(jīng)注明。然?數(shù)據(jù)中還有?些不夠準(zhǔn)確之處：
B.1 概念理解
很多數(shù)據(jù)載體（?頁、報表等）只提供??版，故依據(jù)機器翻譯，數(shù)據(jù)理解上可能存在部分出?。
有些城市不提供 GDP 的精確數(shù)據(jù)，只提供了總收?、總?出。模型中的這部分?jǐn)?shù)據(jù)根據(jù)國際通?的 GDP 計算?式計算得出。
B.2 橫須賀市
橫須賀市政府沒有說明有任何?員傷亡情況（參考?獻 10），擬合過程中暫計為 0。
B.3 仙臺市
仙臺市的?些災(zāi)難結(jié)果數(shù)據(jù)（?員傷亡、財產(chǎn)損失）?較違背常理，在第六節(jié)針對仙臺的分析中已提及。
