基于Mapping Toolbox的地震圖件產(chǎn)出

歐陽小波，湯曜偉

（廣東省地震局，廣東 廣州 510000）

0 引言

Matlab 是由MathWorks 公司推出的一款強(qiáng)大的數(shù)值計算軟件[1]，是數(shù)值計算、符號運(yùn)算和圖形處理等多種功能的強(qiáng)有力實現(xiàn)工具，近幾年來Matlab 這個強(qiáng)大的科學(xué)計算軟件包已經(jīng)得到業(yè)界的廣泛認(rèn)可，并已深入到了各行各業(yè)的眾多學(xué)科、在各大公司、科研機(jī)構(gòu)、大學(xué)校園得到了廣泛應(yīng)用。其具有以下優(yōu)點(diǎn)：（1）代碼簡潔直觀，以矩陣為運(yùn)算基礎(chǔ)，內(nèi)含有豐富的函數(shù)庫，在計算方面避免使用者花費(fèi)大量時間去編寫不必要的小程序；（2）數(shù)值計算方面，科學(xué)計算功能強(qiáng)大，涵蓋了從解微分方程，數(shù)字信號處理FFT，到最新的最優(yōu)化技術(shù)諸如神經(jīng)網(wǎng)絡(luò)、遺傳算法等；（3）數(shù)據(jù)可視化方面，Matlab提供了多種基本圖件處理函數(shù)，如plot、scatter、contour 等，可以很容易的繪制出高質(zhì)量的工程圖譜；（4）Matlab 提供了強(qiáng)大的編程接口，支持Matlab 與其他編程語言進(jìn)行數(shù)據(jù)交換。在地理信息數(shù)據(jù)可視化應(yīng)用方面，Matlab 推出了專用于地理信息圖件繪制的工具箱Mapping toolbox，在Matlab 官方的Mapping toolbox user guide 中是這樣介紹Mapping toolbox 的，“地圖工具箱包含了用于在Matlab 環(huán)境下創(chuàng)建地圖、顯示、分析或操作地理空間數(shù)據(jù)的一系列擴(kuò)展函數(shù)庫和圖形用戶接口（GUI），用戶能夠輕松的綜合多種類型的不同地理數(shù)據(jù)源并根據(jù)它們的空間位置關(guān)系完美的展示出來。總的來說，Mapping toolbox 函數(shù)庫提供了以下主要功能：（1）地理空間數(shù)據(jù)的導(dǎo)入和接入；（2）地理圖像與數(shù)據(jù)的網(wǎng)格化處理；（3）地圖在不同地理坐標(biāo)的投射；（4）地圖的展示和互動函數(shù)（如抓取要素、點(diǎn)擊位置）；（5）支持向量及柵格數(shù)據(jù)的地理信息計算等。由此可見Mapping toolbox 在與地理信息有關(guān)的數(shù)據(jù)的可視化功能方面十分的強(qiáng)大。

在防震減災(zāi)工作領(lǐng)域，無論是地震監(jiān)測預(yù)報還是災(zāi)害防御、應(yīng)急救援。數(shù)據(jù)的可視化工作有著不可替代的作用，在地震監(jiān)測預(yù)報科研方面，震中分布圖、震源深度分布圖、結(jié)合研究地區(qū)的地質(zhì)構(gòu)造狀況，地震研究人員可以分析與研究地震與對應(yīng)地區(qū)位置地質(zhì)構(gòu)造之間的關(guān)系，地震熱點(diǎn)地區(qū)地震的時空遷移變化趨勢，探查并了解其中蘊(yùn)含的規(guī)律，為今后的地震預(yù)測研究積累經(jīng)驗以及區(qū)域震害防御提供指導(dǎo)意見。而在災(zāi)害防御、應(yīng)急救援方面，地震烈度區(qū)劃圖又為抗震設(shè)防提供重要依據(jù)。鑒于Matlab 自帶的強(qiáng)大的數(shù)據(jù)科學(xué)計算能力以及Mapping toolbox 提供的多種地圖繪制功能，本文嘗試以Mapping toolbox 作為出圖工具，針對幾種常見地震圖件，根據(jù)其產(chǎn)出規(guī)范要求，利用Matlab 繪制圖件，并給出對應(yīng)實現(xiàn)代碼。作為Matlab 這種強(qiáng)大的第四代計算機(jī)語言在地震圖件產(chǎn)出中的研究探索。

1 Mapping toolbox 繪制地震圖件流程

從地震圖件的角度來看，有文獻(xiàn)將地震圖件大體分成了兩類，第一類是反映某個或者某幾個地震專題要素的空間分布狀態(tài)，例如歷史地震分布圖、臺站分布圖、重點(diǎn)目標(biāo)及危險源分布圖等。第二類是以地震專題統(tǒng)計數(shù)據(jù)為基礎(chǔ)的統(tǒng)計地圖，用于反映專題信息隨著地理位置變化而分布的情況。第一類地震圖件一般采用散點(diǎn)圖的表現(xiàn)形式，第二類圖件一般采用等高線、分級色彩圖等表現(xiàn)形式[2]。對于上面的二類地震圖件，不管如何變化，在Maltab 的Mapping toolbox 繪制過程中，均需要遵循以下幾個流程來繪制圖件。步驟1：基礎(chǔ)繪圖數(shù)據(jù)導(dǎo)入；步驟2：地理底圖的導(dǎo)入；步驟3：使用各種繪圖函數(shù)（點(diǎn)、線、面）等繪制圖形；步驟4：圖形其他要素添加以及標(biāo)注等。下面就這兩類圖件按照這四個步驟繪圖中出現(xiàn)的技術(shù)難點(diǎn)加以介紹。

1.1 地理底圖導(dǎo)入

Mapping toolbox 對目前主流的地理信息數(shù)據(jù)諸如shapefile、GeoTIFF 和KML 都有著對應(yīng)的數(shù)據(jù)導(dǎo)入和寫函數(shù)，另外，在手頭資源不充足的情況下，Mapping toolbox 也提供了網(wǎng)絡(luò)地圖服務(wù)器數(shù)據(jù)導(dǎo)入功能。本文以美國環(huán)境系統(tǒng)研究所公司（ESRI）開發(fā)的一種空間數(shù)據(jù)開放格式shapefile 為例，說明Mapping toolbox 地理底圖導(dǎo)入操作，shaperead 函數(shù)的基本操作是Map=shaperead（filename），直接讀入名為filename 的shapefile 文件，返回值Map 是一個包含大量數(shù)據(jù)的地理數(shù)據(jù)結(jié)構(gòu)體。shaperead 函數(shù)第二種讀入形式為Map=shaperead（filename，Name，Value，...）Names 是地理數(shù)據(jù)結(jié)構(gòu)的其他屬性部分，Value 是Name 的值，通過指定Name 和Value 值，可以有選擇性導(dǎo)入某區(qū)域地圖或某個地理要素到matlab 運(yùn)算空間。然后再通過函數(shù)mapshow（Map） 就可以將導(dǎo)入到空間中的地理底圖顯示出來。

1.2 基礎(chǔ)繪圖數(shù)據(jù)導(dǎo)入轉(zhuǎn)換

前述說過，Matlab 的Mapping toolbox 能夠輕松的綜合多種類型的不同地理數(shù)據(jù)源并根據(jù)它們的空間位置關(guān)系完美的展示出來，在Mapping toolbox 中，向量數(shù)據(jù)與柵格數(shù)據(jù)是其主要處理對象，對于非地理行業(yè)內(nèi)規(guī)范的地理空間數(shù)據(jù)源如我們地震行業(yè)的地震編目數(shù)據(jù)，按Mapping toolbox 中對于地理空間數(shù)據(jù)的定義來說，地震編目數(shù)據(jù)可以稱之為一種地理空間數(shù)據(jù)，對如不同的非規(guī)范化的數(shù)據(jù)源，Mapping toolbox 提供了一種地理數(shù)據(jù)結(jié)構(gòu)體的結(jié)構(gòu)來統(tǒng)一創(chuàng)建和顯示非規(guī)范的數(shù)據(jù)，這也是一種向量數(shù)據(jù)，在這種Mapping toolbox 規(guī)定的地理數(shù)據(jù)結(jié)構(gòu)體中，有幾個成員要素是必須的，第一個是表征幾何形狀的geometry，其值分別是點(diǎn)Point，多點(diǎn)MultiPoint，線Line 或者多邊形Polygon。第二個是經(jīng)緯Lat 和Lon，最后是補(bǔ)充屬性特征如地名、街道名等，對應(yīng)編目數(shù)據(jù)的屬性就是震級大小。

1.3 重要的繪制函數(shù)

對于第一類地震圖件，使用mapshow 與geoshow 函數(shù)能夠很好的完成創(chuàng)建好的地理數(shù)據(jù)結(jié)構(gòu)體的繪制，matlab 也提供了plotm、linem等函數(shù)完成點(diǎn)、線的繪制，對于等高線可以使用contourm 函數(shù)來完成該類線條的繪制工作。使用contourm 函數(shù)需要的是網(wǎng)格化數(shù)據(jù)，對于給定的經(jīng)緯度范圍，需要用到網(wǎng)格創(chuàng)建函數(shù)meshgrid，其調(diào)用格式如下：［X，Y］=meshgrid（x，y）其涵義為通過數(shù)據(jù)重復(fù)在一維數(shù)組X，Y 的每一個交叉點(diǎn)上創(chuàng)建網(wǎng)格點(diǎn)，當(dāng)X 和Y 都是長度為n 的一維數(shù)組時，產(chǎn)生的X，Y 就是x，y 網(wǎng)格化之后的二維矩陣，在實際的地震烈度圖件繪制時候，輸入與x，y 對應(yīng)的經(jīng)緯度數(shù)組序列，即可以產(chǎn)生經(jīng)緯度網(wǎng)格平面。另外在生成網(wǎng)格之后，如果生成的網(wǎng)格較密，那么就需要對原始繪圖數(shù)據(jù)進(jìn)行插值處理，填補(bǔ)原始數(shù)據(jù)是空白的網(wǎng)格點(diǎn)，特別是對于第二類地震圖件，再利用mesh，surf，contour 等繪圖時，當(dāng)利用meshgrid 函數(shù)完成經(jīng)緯度網(wǎng)格化之后，由于觀測數(shù)據(jù)是散亂分布的，想要獲得平滑的數(shù)據(jù)插值，就需要進(jìn)行插值運(yùn)算。其調(diào)用格式為：Z=girddata（x，y，z，X，Y，method） x，y，z 是原始地震數(shù)據(jù)點(diǎn)，X，Y 是由x，Y 利用meshgrid 網(wǎng)格化得到的網(wǎng)格坐標(biāo)，method 指定了插值方法，Matlab 提供了linear（線性），cubic（立方體），nearest（最近鄰點(diǎn)）和v4（Matlab4 網(wǎng)格點(diǎn)）四種插值算法，其中返回值Z 是插值結(jié)果，用cubic 與v4 得到的插值曲面連續(xù)光滑，而linear 和nearest 則不連續(xù)，默認(rèn)情況下Matlab 會采用linear 線性插值。

1.4 圖件的修飾

Matlab 提供了各類圖件修飾的函數(shù)legend、colorbar 等來給圖件進(jìn)行要素標(biāo)識，在震中分布圖中可以用legend 函數(shù)標(biāo)識不同震級大小的地震。clabel、colorbar 顏色漸變條可以用于等震線的繪制。

2 實例繪制

2.1 繪制震中分布圖

地震震中分布圖是地震圖件中一種重要的基礎(chǔ)地震圖件，用于繪制地震震中分布圖的基礎(chǔ)數(shù)據(jù)來源主要是各省級臺網(wǎng)產(chǎn)出的地震目錄。目前其最新的行業(yè)標(biāo)準(zhǔn)規(guī)范是中華人民共和國地震行業(yè)標(biāo)準(zhǔn)（DB/T66—2016），已于2017 年05 年01 日開始實施，在其規(guī)范性附錄中對CSF 地震目錄和地震觀測報告交換格式有明確詳細(xì)的范例[3]。如圖1 所示。

圖1 CSF 地震目錄標(biāo)準(zhǔn)格式范例Fig.1 The example of CSF seismic catalog standard format

按照前述的繪圖流程，根據(jù)地震目錄格式規(guī)范和Mapping toolbox 的要求，我們要進(jìn)行繪圖工作，首先要將規(guī)范的地震目錄格式轉(zhuǎn)換成符合Mapping toolbox 要求的地理數(shù)據(jù)結(jié)構(gòu)體。地震震中分布圖都是采用散點(diǎn)圖的形式，其結(jié)構(gòu)體geometry 值應(yīng)選擇point，建立一個SeismicInformation 的結(jié)構(gòu)體存儲地震目錄中的關(guān)鍵數(shù)據(jù)，現(xiàn)將代碼展示如下：

［a，b］=size（catalog）；%導(dǎo)入之后目錄長度

［SeismicInformation（1：a）.Geometry］ = deal（′Point′）；%創(chuàng)建了一個SeismicInformation 結(jié)構(gòu)體，其Geometry 成員值為點(diǎn)Point

［SeismicInformation（1：a）.lat］=catalog（a，（27：33））；% 導(dǎo)入目錄中的經(jīng)緯度

［SeismicInformation（1：a）.lon］=catalog（a，（35：42））；

［SeismicInformation（1：a）.Mag］=catalog（a，（47：50）；%結(jié)構(gòu)體補(bǔ)充屬性震級

這樣一個完整的用于繪制震中分布圖的地理數(shù)據(jù)結(jié)構(gòu)體就創(chuàng)建完成。

然后再導(dǎo)入地理底圖，在這里使用shaperead函數(shù)導(dǎo)入整個中國區(qū)域地理底圖，并制定顯示區(qū)域（此處顯示廣東?。?/p>

ChinaMap=shaperead（′china.shp′）；%讀入整個中國地圖

mapshow（ChinaMap）；%軟件初始化顯示整個中國地圖

axis off 消除坐標(biāo)軸背景

size=length（ChinaMap）；%q 確定中國地圖結(jié)構(gòu)體數(shù)組長度

Province=’guangdong’；

cla；

for i=1：size %搜素整個地圖結(jié)構(gòu)體數(shù)組

if strcmp（ChinaMap（i）. Name，Province）==1

mapshow（ChinaMap（i））；%僅顯示整個中國地圖中為guangdong 省的區(qū)域

hold on；

end

end

再按震級大小分別顯示SeismicInformation結(jié)構(gòu)體

for i=1：length（SeismicInformation）

if SeismicInformation（i）. Mag＜=1

mapshow（SeismicInformation，′Marker′，′o′，′MarkerSize′，4，′MarkerFaceColor′，′c′，′MarkerEdgeColor′，′k′）；

........

end

end

本文利用廣東臺網(wǎng)2018 年6 月至12 月數(shù)據(jù)繪制震中分布圖如圖2 所示。

圖2 地震震中分布圖例子Fig.2 The example of epicenter distribution

2.2 繪制地震烈度圖

地震烈度圖是描述地震引起地面震動及其影響程度的圖件，一般用等震線描述，等震線又是指不同地震烈度和地震震動強(qiáng)度的分界線，地震烈度的具體數(shù)據(jù)一般按照地震烈度表由現(xiàn)場宏觀考察得到或者由強(qiáng)震計地面峰值加速度值得到，由于筆者未有相關(guān)地震烈度數(shù)據(jù)來源，利用地震發(fā)生之后各臺站地震計速度峰值數(shù)據(jù)代替作圖，為顯示圖件效果。前述地理底圖的代碼已展示，此處主要展示等值線代碼。

lat= Stationlocatio（：，1）；%臺站緯度

longt= Stationlocatio（：，2）；%臺站經(jīng)度

Value=Stationlocatio（：，3）；%地震峰值加速度值

xi=linspace（min（lat），max（lat），50）；

yi=linspace（min（longt），max（longt），50）；

［xi，yi］=meshgrid（xi，yi）；%按50 個點(diǎn)網(wǎng)格化所有經(jīng)緯度

zi=griddata（lat，longt，Value，xi，yi，′linear′）；%網(wǎng)格化之后線性插值

［c，h］=contour（yi，xi，zi）；%繪制等震線

Clabel；%在等震線上標(biāo)注數(shù)值

Colorbar；%加顏色漸變條

作圖其效果如圖3 所示。

圖3 地震烈度等震線示例圖Fig.3 The isoseismal line of the seismic intensity

3 結(jié)語

本文通過以上的兩類不同類型地震圖件的專題繪制研究，闡述了利用Matlab 繪制圖件的四個基本步驟，展示了Matlab 的Mapping toobox強(qiáng)大繪圖功能以及Matlab 極為簡潔的編程語言,經(jīng)過實例驗證，Matlab 能夠繪制出滿足要求的具有較好表現(xiàn)效果的地震圖件，更好的為防震減災(zāi)工作服務(wù)。最后，需要說明的是在日常的地震圖件繪圖中，地震圖件的質(zhì)量很大程度上依賴于自備的地理底圖庫，即地理信息庫的豐滿程度，因此對于繪圖者來說，建立和保存自身常用的制圖底圖模板非常重要。