石富強(qiáng)，戴 勇，張 博

（1.陜西省地震局，西安 710068；2.內(nèi)蒙古自治區(qū)地震局，呼和浩特 010010；3.甘肅省地震局，蘭州 730000）

0 引言

歷史震例的回溯性研究是地震分析預(yù)測(cè)的學(xué)習(xí)過(guò)程，也是對(duì)歷史震例重新認(rèn)識(shí)、不斷升華的過(guò)程。廣大地震科技工作者在長(zhǎng)期的監(jiān)測(cè)預(yù)報(bào)中反復(fù)推敲歷史震例，尋找其中的共性問(wèn)題和個(gè)性問(wèn)題，并做了大量的探索、挖掘和深化研究工作［1－8］。在地震預(yù)報(bào)還處在探索階段的今天，這項(xiàng)工作還將持續(xù)下去。但同時(shí)也面臨著另一個(gè)問(wèn)題：過(guò)去很多歷史震例，盡管在當(dāng)時(shí)震例總結(jié)時(shí)給出了很好的曲線(xiàn)形態(tài)，也匯編成冊(cè)（如：《中國(guó)震例》）。但由于時(shí)間久遠(yuǎn)、觀(guān)測(cè)手段革新、觀(guān)測(cè)儀器更新?lián)Q代、觀(guān)測(cè)臺(tái)點(diǎn)更換、數(shù)據(jù)庫(kù)升級(jí)換代等等問(wèn)題，過(guò)去好多臺(tái)站的原始觀(guān)測(cè)數(shù)據(jù)已經(jīng)丟失，這無(wú)疑給震例再認(rèn)識(shí)帶來(lái)巨大困難。如何從有限的《中國(guó)震例》圖片中提取高質(zhì)量的前兆數(shù)據(jù)，是震例再研究的一項(xiàng)重要的基礎(chǔ)性工作。另外，提取數(shù)據(jù)的質(zhì)量直接決定著異常再認(rèn)識(shí)的科學(xué)性和準(zhǔn)確性。

目前互聯(lián)網(wǎng)上關(guān)于曲線(xiàn)矢量化的程序也很多，如：Origin里的Digitize.opk插件，R2V32，Engauge Digitizer，GetdataDraghDigitizer，Graph Digitizer Scout以及Image2data等等。但作者在試用時(shí)發(fā)現(xiàn)存在以下問(wèn)題：①這些軟件針對(duì)的都是橫軸為實(shí)數(shù)時(shí)間序列的曲線(xiàn)，而震例里面的時(shí)間序列，橫軸基本都是年月日（如：20080512）式的時(shí)間序列，用這些軟件提取出來(lái)的數(shù)據(jù)還需要進(jìn)一步轉(zhuǎn)化才能方便使用；②這些程序多需要手動(dòng)逐點(diǎn)雙擊拾取點(diǎn)，非常累人且效率較低（如Digitize.opk和R2V32），同時(shí)取出的數(shù)據(jù)點(diǎn)主觀(guān)依賴(lài)性較大；③在自動(dòng)矢量化方面這些程序僅能針對(duì)相對(duì)單一且干凈的圖片，震例里面的圖片大多相對(duì)復(fù)雜，有多條曲線(xiàn)同圖對(duì)比或不間斷標(biāo)地震事件以及“儀器檢修、調(diào)零”等標(biāo)注字樣，用這些軟件自動(dòng)矢量化時(shí)還需要借助第三方軟件（如：Photoshop）手動(dòng)擦除標(biāo)注文字及多余曲線(xiàn)，增加工作量且效率不高。

為此，本文作者基于Matlab編寫(xiě)了一種能夠快速?gòu)默F(xiàn)有圖片文件中提取歷史震例前兆曲線(xiàn)的程序，并將其命名為Quick＿Pick（程序下載地址：http：／／pan.baidu.com／s／1c0vkmze）。

1 程序基本原理

圖像的亮度值直接存于圖像數(shù)組中。對(duì)于RGB圖像，該數(shù)組為M×N×3階矩陣；對(duì)于灰度圖像，該數(shù)組為M×N 階矩陣。其中，M、N 表示圖像像素的行、列數(shù)。矩陣中的元素代表該點(diǎn)的亮度值，其范圍為［0，255］，分別代表黑色和白色。

1.1 圖片的讀取及曲線(xiàn)數(shù)據(jù)篩選

程序首先讀入需要被矢量化的圖片，并利用matlab命令“imadjust”提高圖像對(duì)比度。一般圖片中，曲線(xiàn)基本為黑色或者彩色，而背景基本為白色，因此我們?cè)诔绦蚶镌O(shè)置一個(gè)閾值（如C＝254），將亮度值小于閾值C 的像素點(diǎn)亮度全部修改為“0”（黑色），將亮度值大于閾值C 的像素點(diǎn)亮度全部修改為“255”（白色）。這樣可以很容易從巨大的像素矩陣中找出疑似曲線(xiàn)像素點(diǎn)的集合R，當(dāng)然其中可能包含其他曲線(xiàn)或者地震標(biāo)注等其他文字信息。為了進(jìn)一步剔除其他曲線(xiàn)以及多余的文字標(biāo)注，我們?cè)陲@示的像素點(diǎn)集R 圖像上手動(dòng)拾取一個(gè)能夠包圍目標(biāo)曲線(xiàn)的多邊形框，利用matlab自帶的“inpolygon”命令自動(dòng)拾取多邊形內(nèi)的像素點(diǎn)坐標(biāo)（圖2d）。至此，基本完成了數(shù)據(jù)的篩選。

1.2 像素點(diǎn)坐標(biāo)向時(shí)間軸（T 軸）和實(shí)數(shù)軸（y軸）的轉(zhuǎn)換

對(duì)于橫軸為實(shí)數(shù)序列的曲線(xiàn)來(lái)說(shuō)，將1.1節(jié)拾取的像素點(diǎn)坐標(biāo)直接做線(xiàn)性等比例轉(zhuǎn)換便可以得到曲線(xiàn)矢量化后的數(shù)據(jù)。但對(duì)于橫軸時(shí)間序列為“YYMMDD”格式的地震前兆觀(guān)測(cè)數(shù)據(jù)，則會(huì)涉及到閏年（366天）和平年（365天）等問(wèn)題。為此，我們利用matlab 自帶命令“datenum”將指定日期“YYMMDD”轉(zhuǎn)換為實(shí)數(shù)時(shí)間序列；然后再將橫軸和縱軸的像素點(diǎn)坐標(biāo)（i，j）分別線(xiàn)性等比例轉(zhuǎn)化為實(shí)數(shù)時(shí)間序列坐標(biāo)下的（t，y）；最后將橫軸數(shù)據(jù)轉(zhuǎn)到“YYMMDD”格式下的（T，y）。至此，基本完成了圖像的初步矢量化工作并提取保存了曲線(xiàn)數(shù)據(jù)。

1.3 數(shù)據(jù)平滑處理

圖1 直接提取的像素點(diǎn)坐標(biāo)及其簡(jiǎn)單平均結(jié)果

上述過(guò)程得到的曲線(xiàn)形態(tài)在整體上與原始圖像非常相像，但放大到細(xì)節(jié)，會(huì)出現(xiàn)如圖1 的鋸齒形態(tài)。這是由于原始圖片曲線(xiàn)的粗度占用多個(gè)像素點(diǎn)造成在一個(gè)橫坐標(biāo)i處拾取了多個(gè)（n）縱坐標(biāo)j值，這將給頻譜分析等添加不可估量的數(shù)據(jù)成分。為避免這種影響，一般在i坐標(biāo)下對(duì)多個(gè)j坐標(biāo)求平均值得到光滑曲線(xiàn)［9］。但這種求平均在遇到“毛刺”、“突跳”等單個(gè)脈沖型信號(hào)時(shí)，會(huì)直接將原來(lái)的脈沖幅度折半，有失信號(hào)的真實(shí)性。在現(xiàn)實(shí)中，“毛刺”、“突跳”等單個(gè)脈沖型信號(hào)往往又被認(rèn)為是可靠的前兆異常［10－13］。因此，有必要盡可能保留并還原這種信號(hào)。

為此，我們將前面1.1節(jié)識(shí)別出來(lái)的時(shí)間序列通過(guò)滑動(dòng)平均將曲線(xiàn)置于“0”附近，對(duì)于給定的橫軸像素點(diǎn)坐標(biāo)i計(jì)算拾取的n個(gè)點(diǎn)縱坐標(biāo)平均值A(chǔ)＝mean（y（i，1）：y（i，n））。當(dāng)A＜0 時(shí)，我們認(rèn)為曲線(xiàn)向下脈沖，取i坐標(biāo)下的最小值min（y（i，1：n））為脈沖峰值；當(dāng)A＞0時(shí)，我們認(rèn)為曲線(xiàn)向上脈沖，取i坐標(biāo)下的最大值max（y（i，1：n））為脈沖峰值；當(dāng)A＝0時(shí)，曲線(xiàn)光滑沒(méi)有脈沖信號(hào)，自動(dòng)服從前面提到的平均光滑處理，信號(hào)值為mean（y（i，1：n））。至此，我們?nèi)客瓿闪藬?shù)據(jù)的篩選、轉(zhuǎn)換以及平滑處理工作。下面，簡(jiǎn)單介紹程序的使用流程，以便讀者有直觀(guān)深入的理解。

2 程序操作流程簡(jiǎn)介

假設(shè)我們?cè)诠ぷ髦杏龅綀D2所示的前兆數(shù)據(jù)曲線(xiàn)。由于各種原因該測(cè)項(xiàng)原始數(shù)據(jù)已經(jīng)丟失，為了研究和深挖掘的需要，我們必須盡可能還原原始數(shù)據(jù)，從圖1中提取較為可靠的數(shù)據(jù)來(lái)替代原始數(shù)據(jù)做深入分析。首先，我們將圖片截取放在程序目錄下面（圖3a），運(yùn)行“Quick＿Pick”讀入圖片并按照?qǐng)D3b“1”、“2”、“3”順序，單擊鼠標(biāo)左鍵選定坐標(biāo)3個(gè)控制點(diǎn)，程序自動(dòng)跳出圖3c輸入框；根據(jù)前面控制點(diǎn)“1”、“2”、“3”分別輸入X 軸起始點(diǎn)坐標(biāo)“20090801”和“20090831”以及Y 軸起始點(diǎn)“58.10”和“58.45”后點(diǎn)“確定”鍵，此時(shí)程序會(huì)將圖片上灰度值為“0”的全部可疑信息篩選出來(lái)顯示為圖3d；再單擊鼠標(biāo)左鍵生成一個(gè)盡可能簡(jiǎn)單的閉合多邊形區(qū)域（順時(shí)針或逆時(shí)針均可），使所有有用信息全部落于多邊形區(qū)域內(nèi)，排除“降雨”等標(biāo)記的多余信息（多邊形閉合的時(shí)候，雙擊鼠標(biāo)左鍵）。這樣，我們便完成了一張無(wú)原始數(shù)據(jù)圖片的矢量化工作。提取的數(shù)據(jù)將自動(dòng)保存于程序目錄下“某臺(tái)某測(cè)項(xiàng)2009年8月觀(guān)測(cè)值.dat”文件中（圖3e），可以直接在Mapsis日常分析軟件下使用。

圖2 某測(cè)項(xiàng)2009年8月觀(guān)測(cè)值圖片文件

圖3 程序使用流程

3 程序檢驗(yàn)

由于是事先假定的曲線(xiàn)，在矢量化完成后，我們可以拿出原始數(shù)據(jù)和提取的數(shù)據(jù)做對(duì)比，將原始數(shù)據(jù)和提取的數(shù)據(jù)在同比例的坐標(biāo)中繪圖，如圖4所示。對(duì)比發(fā)現(xiàn)，提取的數(shù)據(jù)在長(zhǎng)尺度上與原始數(shù)據(jù)形態(tài)幾乎完全一致。但個(gè)別突跳點(diǎn)（如18日晚上19日凌晨的突跳）幅度略有降低，這是由于個(gè)別像素點(diǎn)灰度值不夠被限所致。將18日至20日數(shù)據(jù)單獨(dú)顯示如圖5所示，可以發(fā)現(xiàn)：在原始曲線(xiàn)較光滑的時(shí)候，提取數(shù)據(jù)與原始數(shù)據(jù)完全一致；在原始曲線(xiàn)變化不穩(wěn)的時(shí)候，略有差距。這是因?yàn)椋紨?shù)據(jù)變化不穩(wěn)定、上下波動(dòng)時(shí)，在圖2中像素點(diǎn)已經(jīng)靠在一起或者重疊，而程序很難正確將他們分離。但整體來(lái)說(shuō)，提取的信號(hào)基本符合原始曲線(xiàn)的變化規(guī)律，能夠替代原始曲線(xiàn)做近似化的深入分析。

4 總結(jié)

本文基于Matlab設(shè)計(jì)了一組快速矢量化圖片的程序，該程序在歷史觀(guān)測(cè)數(shù)據(jù)缺失的情況下，研究歷史觀(guān)測(cè)數(shù)據(jù)或者從已發(fā)表論文中提取前人數(shù)據(jù)做深挖掘研究具有很強(qiáng)的實(shí)用性，且操作簡(jiǎn)便。

對(duì)于過(guò)于老的圖片（像素不高，分辨率極低，噪聲成分很高），本文還沒(méi)有考慮，這種情況需要通過(guò)濾波等技術(shù)做圖像再處理。Matlab自身也有很強(qiáng)的圖像處理能力［14］，在遇到這種情況的時(shí)候讀者可以將Matlab圖像處理和本程序一起使用。

圖4 提取數(shù)據(jù)與原始數(shù)據(jù)的對(duì)比（全部）

圖5 提取數(shù)據(jù)與原始數(shù)據(jù)的對(duì)比（局部）

另注：該程序盡量在高版本的Matlab下使用，作者是在Matlab R2011b平臺(tái)下編寫(xiě)。

致謝：非常感謝“2014年青年分析預(yù)報(bào)人員培訓(xùn)班”全體同學(xué)，他們給予了本文作者莫大的支持和鼓勵(lì)，并提出許多優(yōu)秀的建議和意見(jiàn)！

［1］ 蔡仲瓊.新疆烏魯木齊地區(qū)地震水化震例深入剖析［J］.內(nèi)陸地震，1989，3（3）：230－239.

［2］ 顧瑾萍，張曉東，黃輔瓊，等.我國(guó)震例前兆異常統(tǒng)計(jì)特征和應(yīng)用研究［J］.地震，2004，24（2）：59－65.

［3］ 蔣海昆，苗青壯，吳瓊，宋金.基于震例的前兆統(tǒng)計(jì)特征分析［J］.地震學(xué)報(bào)，2009，31（3）：245－259.

［4］ 簡(jiǎn)春林，晏銳，黃輔瓊.中國(guó)大陸23個(gè)震例流體異常時(shí)空分布特征研究［J］.地震，2006，26（1）：99－106.

［5］ 王道.哈薩克斯坦及相鄰地區(qū)地下流體震例資料［J］.內(nèi)陸地震，2004，17（4）：371－384.

［6］ 衛(wèi)定軍，羅國(guó)富，司學(xué)蕓，等.基于支持向量機(jī)回歸的寧夏地震前兆［J］.地震研究，2014，37（2）：186－191.

［7］ 張桂銘，劉文鋒.基于震例研究的地震預(yù)測(cè)預(yù)報(bào)分析［J］.中國(guó)地震，2013（4）：528－536.

［8］ 鄭兆苾，張國(guó)民，何康，等.中國(guó)大陸地震震例異常統(tǒng)計(jì)與分析［J］.地震，2006，26（2）：29－37.

［9］ 付昆昆，鄭百林，李鑫.基于Matlab的圖像曲線(xiàn)數(shù)據(jù)提取方法［J］.汕頭大學(xué)學(xué)報(bào)（自然科學(xué)版），2010，25（2）：51－63.

［10］ 盧雙苓，于慶民，鐘普浴，等.泰安地震臺(tái)SSQ－2數(shù)字水平擺曲線(xiàn)畸變分析［J］.大地測(cè)量與地球動(dòng)力學(xué)，2010，30（增）：53－61.

［11］ 牛安福，張凌空，閆偉，等.中國(guó)鉆孔應(yīng)變觀(guān)測(cè)能力及在地震預(yù)報(bào)中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31（2）：48－52.

［12］ 邱澤華，張寶紅，池順良，等.汶川地震前姑咱臺(tái)觀(guān)測(cè)的異常應(yīng)變變化［J］.中國(guó)科學(xué)：地球科學(xué)，2010，40（8）：1031－1039.

［13］ 田山，汪翠枝，李君英，等.數(shù)字前兆觀(guān)測(cè)地震異常短臨特征的初步研究［J］.地震地磁觀(guān)測(cè)與研究，2009，30（2）：57－67.

［14］ 成波.數(shù)字圖像處理及MATLAB實(shí)現(xiàn)［M］.重慶：重慶大學(xué)出版社，2003.