震例研究中前兆數(shù)據(jù)曲線快速矢量化程序

石富強，戴 勇，張 博

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

0 引言

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

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

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

1 程序基本原理

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

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

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

1.2 像素點坐標向時間軸（T 軸）和實數(shù)軸（y軸）的轉換

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

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

圖1 直接提取的像素點坐標及其簡單平均結果

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

為此，我們將前面1.1節(jié)識別出來的時間序列通過滑動平均將曲線置于“0”附近，對于給定的橫軸像素點坐標i計算拾取的n個點縱坐標平均值A＝mean（y（i，1）：y（i，n））。當A＜0 時，我們認為曲線向下脈沖，取i坐標下的最小值min（y（i，1：n））為脈沖峰值；當A＞0時，我們認為曲線向上脈沖，取i坐標下的最大值max（y（i，1：n））為脈沖峰值；當A＝0時，曲線光滑沒有脈沖信號，自動服從前面提到的平均光滑處理，信號值為mean（y（i，1：n））。至此，我們全部完成了數(shù)據(jù)的篩選、轉換以及平滑處理工作。下面，簡單介紹程序的使用流程，以便讀者有直觀深入的理解。

2 程序操作流程簡介

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

圖2 某測項2009年8月觀測值圖片文件

圖3 程序使用流程

3 程序檢驗

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

4 總結

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

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

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

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

另注：該程序盡量在高版本的Matlab下使用，作者是在Matlab R2011b平臺下編寫。

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

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

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

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

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

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

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

［7］ 張桂銘，劉文鋒.基于震例研究的地震預測預報分析［J］.中國地震，2013（4）：528－536.

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

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

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

［11］ 牛安福，張凌空，閆偉，等.中國鉆孔應變觀測能力及在地震預報中的應用［J］.大地測量與地球動力學，2011，31（2）：48－52.

［12］ 邱澤華，張寶紅，池順良，等.汶川地震前姑咱臺觀測的異常應變變化［J］.中國科學：地球科學，2010，40（8）：1031－1039.

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

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