東北大學(xué) 2017 年大學(xué)生數(shù)學(xué)建模競賽 A 題一等獎(jiǎng)成果弥鹦。
在 GitHub 上：https://github.com/lalxyy/NEU-2017-MCM-A
由于不能顯示行內(nèi)公式赠摇，部分 LaTeX 公式以原形式展現(xiàn)芍瑞。
摘要
東京時(shí)間 2011 年 3 ? 11 ? 14 時(shí) 46 分，?本東北地?太平洋海域發(fā)??? 9.0 級地震皮胡，并伴隨有 7.0 到 7.4 級的多次強(qiáng)余震痴颊。9.0 級地震也造成了歷史記錄以來罕見的?達(dá) 39 ?的巨?海嘯。截?到 2011 年 5 ?屡贺，超過 2.4 萬?被官?報(bào)告為死亡或失蹤蠢棱；另外，福島第?核電站因此發(fā)?的嚴(yán)重核事故也導(dǎo)致了千余?的?直接死亡甩栈，這些?并發(fā)?的重?事故也給?本及周邊國家?guī)砹藝?yán)重的損失泻仙。本?通過評估磐城、仙臺量没、宮古與橫須賀的地震破壞及損失情況玉转，利?多元線性回歸的?式，建?了災(zāi)害有關(guān)因素（地震強(qiáng)度殴蹄、海嘯等）與災(zāi)害造成的后果（?員傷亡冤吨、經(jīng)濟(jì)損失）的模型。
本?根據(jù)互聯(lián)?中能夠查得的數(shù)據(jù)進(jìn)?提取饶套，并獲得了沿岸海浪最??度、與震中距離垒探、地震??烈度妓蛮、經(jīng)濟(jì)?平（GDP）等數(shù)據(jù)，以及?員傷亡圾叼、財(cái)產(chǎn)損失的數(shù)據(jù)蛤克。根據(jù)繪制的散點(diǎn)圖確定多元擬合中每個(gè)分量對應(yīng)的多項(xiàng)式次數(shù)，之后將?次單項(xiàng)式換元轉(zhuǎn)化為?次單項(xiàng)式夷蚊，使?多元線性擬合來得出多項(xiàng)式系數(shù)构挤。利?所得出的模型，本?又簡要分析了仙臺市的災(zāi)情情況惕鼓，并且指出數(shù)據(jù)點(diǎn)較少帶來的模型的不?之處筋现。
關(guān)鍵詞： 多項(xiàng)式擬合; 地震破壞性評估; 東?本?地震
1 問題背景與分析
1.1 背景
2011 年 3 ? 11 ??本福島發(fā)??? 9.0 級?地震，強(qiáng)震引發(fā)海嘯，導(dǎo)致當(dāng)時(shí)全球最?的在役核電站——福島核電站放射性物質(zhì)外泄矾飞。該事故根據(jù)國際核事件分級被評為最?級 7 級一膨，與切爾諾貝利核電站泄漏事故等級相同。截? 2014 年 2 ?洒沦，“3.11”東?本?地震和?海嘯死亡的?數(shù)達(dá)到 15884 ?豹绪。 截? 2014 年 5 ?，受東?本?地震和福島第?核電站事故影響申眼，?直接死亡的?數(shù)已經(jīng)達(dá)到 1699 ?瞒津，超過了在事故中直接死亡的 1603 ?。
建?數(shù)學(xué)模型括尸，在以下?本沿海城市：磐城巷蚪、仙臺、宮古和橫須賀 中姻氨。?較各種規(guī)模??的地震破壞以及所造成的后果钓辆，并為當(dāng)?shù)貓?bào)紙準(zhǔn)備?篇?章，解釋根據(jù)你的模型關(guān)于其中?個(gè)城市的發(fā)現(xiàn)肴焊。
1.2 問題分析
該問題本質(zhì)為求解?個(gè)地震因?向量 x = (x_1 , x_2 , ..., x_m) 到影響程度向量 y = (y_1 , y_2 , ..., y_n) 的函數(shù) F（可能包含多個(gè)線性變換）前联，表示為
，怎樣得出 F(x) 也是本?討論的核?內(nèi)容娶眷。
由于災(zāi)害產(chǎn)?的影響與?類?產(chǎn)??平的發(fā)展情況和???的主觀能動(dòng)性有很?的關(guān)聯(lián)似嗤，所以單純憑借歷史以來的地震傷亡?數(shù)、財(cái)產(chǎn)損失情況不能作為衡量當(dāng)今災(zāi)害情況與財(cái)產(chǎn)損失的直接依據(jù)届宠，但是?類??因素又是衡量地震帶來的破壞性的必要考慮因素烁落，且關(guān)系極其復(fù)雜豌注。為了簡化模型伤塌、提?模型準(zhǔn)確性與說服?轧铁，我們將磐城童本、仙臺、宮古和橫須賀四個(gè)地區(qū)在同?時(shí)間橫向?qū)?绑蔫，不需控制時(shí)間的變量，?是通過與震中的距離來衡量地震的影響規(guī)模鄙煤。
題?所給出的四個(gè)地區(qū)的建設(shè)晾匠、經(jīng)濟(jì)發(fā)展?平不盡相同。圖 1 是 Google Maps 的數(shù)據(jù)梯刚，給出了四城的地理位置與震中坐標(biāo)凉馆，最右側(cè)海域中的是本次地震的震中，其余點(diǎn)從上到下依次為宮古亡资、仙臺澜共、磐城、橫須賀锥腻；
以及四城 2010 年當(dāng)時(shí)的國民?產(chǎn)總值（GDP）[1][2][3][4]：
| 城市 | 國民生產(chǎn)總值 |
|---|---|
| 仙臺 | 61.7 |
| 橫須賀 | 12.2882 |
| 磐城 | 16.5905 |
| 宮古 | 5.19 |
表 1: 2010 年四城國民?產(chǎn)總值（單位：億美元）
由此我們將經(jīng)濟(jì)因素（代表抗災(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ì)外泄造成的事故。四個(gè)城市的海嘯浪?度最?值可以在?象數(shù)據(jù)中查得 [5][6]（表 3）经窖； 以及從各城市與所屬縣政府查得的?員傷亡與財(cái)產(chǎn)損失情況（表 4）[7][8][9][10]坡垫。
| 城市 | 海嘯沿岸最??度 / ? |
|---|---|
| 宮古 | 8.5 |
| 仙臺 | 8.6 |
| 磐城 | 6.51 |
| 橫須賀 | 2.06 |
表 3: 觀測到的海嘯?度
注：部分觀測點(diǎn)不以城市命名，選取對應(yīng)最近的觀測點(diǎn)画侣。
| 城市 | 傷亡?數(shù) | 經(jīng)濟(jì)損失（億美元） |
|---|---|---|
| 宮古 | 611 | 22.4775 |
| 仙臺 | 904 | 119.0022 |
| 橫須賀 | 0 | 0.92 |
| 磐城 | 466 | 0.1478 |
表 4: 四城傷亡?數(shù)與經(jīng)濟(jì)損失
注：橫須賀市政府沒有說明有任何?員傷亡情況（參考?獻(xiàn) 10）冰悠，擬合過程中暫計(jì)為 0。（附錄 B.2）
2 模型假設(shè)與說明
2.1 符號說明
| 符號 | 代表意義與說明 |
|---|---|
| A_1 | 仙臺市 |
| A_2 | 橫須賀市 |
| A_3 | 宮古市 |
| A_4 | 磐城市 |
| X | 表示地震整體因素的向量 |
| x_1, ..., x_m | 影響地震破壞性的各個(gè)因素 |
| Y | 表示地震破壞性結(jié)果的向量 |
| y_1, ..., y_n | 地震破壞程度的不同衡量特征（財(cái)產(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：?類的能動(dòng)性（主觀性與?產(chǎn)?）與發(fā)展時(shí)間相關(guān)且滿?映射關(guān)系；
假設(shè) 3：?類在地震中的抗災(zāi)防災(zāi)?平與經(jīng)濟(jì)發(fā)展程度正相關(guān)宪卿；
假設(shè) 4：?個(gè)城市區(qū)域內(nèi)的建筑物強(qiáng)度可以??個(gè)平均值來表示，且可以?建筑物的強(qiáng)度万栅、破壞程度作為衡量城市破壞程度的主要因素佑钾。
3 模型建立
3.1 知識背景
3.1.1 地震影響場
采?的地震烈度影響場衰減公式：
沿長軸?向，
沿短軸?向烦粒，
其中 d 表示城市與震中的直線距離休溶。[12]
3.1.2 海嘯等級
采?太平洋地區(qū)普遍適?的今村－飯?強(qiáng)度分級來評定海嘯的強(qiáng)度等級代赁，其中 H_av 表示沿最近海岸的平均海嘯?度，I 為輸出強(qiáng)度值兽掰。[13]
3.2 具體問題參數(shù)引入
由問題分析所述芭碍，已經(jīng)基本可以確定 X、Y 各分量的表?含義孽尽。將 x_1 , x_2 , y_1 , ... 等分量指派為表 6窖壕、7 中所示的含義；四城中每個(gè)城市都與且僅與?組向量對應(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)濟(jì)損失（億美元） |
表 7: 影響結(jié)果 y 各分量表示含義
并通過 Python 語?與 Matplotlib 庫繪制 y 各分量與 x 各分量散點(diǎn)圖：
3.3 基本模型
根據(jù)現(xiàn)有散點(diǎn)圖及統(tǒng)計(jì)回歸模型，認(rèn)為 y_1 , y_2 為被解釋變量熏挎，x_1 , ..., x_4 為解釋變量速勇。觀察散點(diǎn)圖，認(rèn)為 y_1 與 x_3 坎拐，y_2 與 x_3 滿??次多項(xiàng)式擬合的圖像特征烦磁，且函數(shù)對應(yīng)單調(diào)遞減：
其中 a_{ij}, b_{ij} 為需要求出的多項(xiàng)式系數(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 的圖像特征也可用二次多項(xiàng)式來擬合：
y_1 與 y_2 的取值受到以上因素的影響。所以可以得出關(guān)于 x_2, x_3, x_4 的兩個(gè)三元多項(xiàng)式擬合函數(shù)
磅轻，其中 c_1珍逸，c_2 為該多項(xiàng)式擬合模型中的常數(shù)，也就是
4 模型求解
引?擬合?法中的最??乘法聋溜，使?多元線性回歸的?法求解該模型谆膳。令
，此時(shí)原函數(shù)變?yōu)槲逶€性擬合撮躁。
對于任意一個(gè)取值點(diǎn) (x_{ij}, y_{kj})漱病，其中 j = 0, ..., 3，i = 1, ..., 6把曼，k = 1, 2（已知四個(gè)城市對應(yīng)的數(shù)據(jù)點(diǎn)）杨帽，有
，對于 y_1嗤军，也就是
此時(shí)令
注盈，就有
兩側(cè)同時(shí)左乘 A^T，得
叙赚，整理得
老客，求出 c_1, a_{21}, a_{31}, a_{41}, a_{32}, a_{42}僚饭。
使用計(jì)算機(jī)求出的最后結(jié)果為
5 模型評價(jià)與改進(jìn)
該模型可以根據(jù)已有的?然災(zāi)害數(shù)據(jù)較為準(zhǔn)確地求得財(cái)產(chǎn)損失與?員傷亡的程度，能夠給出范圍內(nèi)的合理值胧砰，可?于防災(zāi)減災(zāi)與災(zāi)害造成后果的快速估計(jì)與判斷鳍鸵。然?本模型利?的數(shù)據(jù)點(diǎn)較少，如果使??量的數(shù)據(jù)來進(jìn)?擬合則會使模型更有說服?尉间。
6 推廣與應(yīng)用：仙臺市
仙臺市是?本東北地?的區(qū)域中?偿乖，也是本題?已知條件中離震中最為接近的城市。從之前繪制的圖像和模型預(yù)估的情況來看乌妒，仙臺的情況并不容樂觀：常識上講汹想，?個(gè)發(fā)達(dá)區(qū)域的經(jīng)濟(jì)發(fā)展?平越?，防災(zāi)減災(zāi)能?也應(yīng)該越強(qiáng)撤蚊，這也是本?在模型建?時(shí)考慮 GDP（國民?產(chǎn)總值）作為防災(zāi)減災(zāi)因素的原因古掏；但是由此看來，仙臺及?本東北地?城市仍然需要加強(qiáng)對于防災(zāi)減災(zāi)?采取的措施侦啸，例如削弱海嘯因素帶來的影響槽唾，可以體現(xiàn)在本模型中的擬合曲線更為平滑，?且海嘯是導(dǎo)致此次地震及連帶事故中?員傷亡與財(cái)產(chǎn)損失的?個(gè)主要的原因光涂。
參考文獻(xiàn)
[1] 宮古市役所. 平成 23 年度 財(cái)政狀況資料集 [EB/OL]. 2011. http://www.city.miyako.iwate.jp/data/open/cnt/3/5677/1/H23zaiseijoukyou.pdf.
[2] 橫須賀市. 財(cái)政?書（平成 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] ?本気象庁. 気象庁｜津波観測點(diǎn)（東北地?）[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)害対策本部週報(bào) [EB/OL]. 2015. http://www.city.iwaki.lg.jp/www/contents/1449132951986/simple/higai0524.pdf.
[10] 橫須賀市. 東?本?震災(zāi)関連情報(bào) [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] 傅再揚(yáng), 危福泉, 林巖釗. 福建地震災(zāi)害損失計(jì)算分析與思考 [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 程序源碼
帶有運(yùn)?結(jié)果的源碼與部分互動(dòng)界?可以通過 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'] # 在其他機(jī)器上可能需要修改字體名
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])
# 人員傷亡齐佳、財(cái)產(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 散點(diǎn)圖繪制
def draw_scatter_pic(x, y, description, x_label, y_label):
"""
繪制散點(diǎn)圖的（統(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, '散點(diǎn)圖：人員傷亡與 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('散點(diǎn)圖：震中直線距離與傷亡人數(shù)')
# plt.xlabel('直線距離')
# plt.ylabel('人員傷亡')
# ax1.scatter(x_distance, y_person, c='r', marker='o')
# plt.show()
draw_scatter_pic(x_distance, y_person, '散點(diǎn)圖：震中直線距離與傷亡人數(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('散點(diǎn)圖：沿岸海浪高度與傷亡人數(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, '散點(diǎn)圖：GDP 與經(jīng)濟(jì)損失', 'GDP', '經(jīng)濟(jì)損失')
draw_scatter_pic(x_distance, y_loss_economy, '散點(diǎn)圖：震中直線距離與經(jīng)濟(jì)損失', '直線距離', '經(jīng)濟(jì)損失')
draw_scatter_pic(x_tsunami_height, y_loss_economy, '散點(diǎn)圖：沿岸海浪高度與經(jīng)濟(jì)損失', '海浪高度', '經(jīng)濟(jì)損失')
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é)者分析報(bào)告私恬，在參考?獻(xiàn)中都已經(jīng)注明。然?數(shù)據(jù)中還有?些不夠準(zhǔn)確之處：
B.1 概念理解
很多數(shù)據(jù)載體（?頁炼吴、報(bào)表等）只提供??版本鸣，故依據(jù)機(jī)器翻譯，數(shù)據(jù)理解上可能存在部分出?硅蹦。
有些城市不提供 GDP 的精確數(shù)據(jù)荣德，只提供了總收?、總?出童芹。模型中的這部分?jǐn)?shù)據(jù)根據(jù)國際通?的 GDP 計(jì)算?式計(jì)算得出涮瞻。
B.2 橫須賀市
橫須賀市政府沒有說明有任何?員傷亡情況（參考?獻(xiàn) 10），擬合過程中暫計(jì)為 0假褪。
B.3 仙臺市
仙臺市的?些災(zāi)難結(jié)果數(shù)據(jù)（?員傷亡署咽、財(cái)產(chǎn)損失）?較違背常理，在第六節(jié)針對仙臺的分析中已提及嗜价。
