中國(guó)地震局地質(zhì)研究所 中國(guó)地震局活動(dòng)構(gòu)造與火山重點(diǎn)實(shí)驗(yàn)室 吳熙彥

基于開(kāi)源GIS的活動(dòng)斷層矢量數(shù)據(jù)融合方案初探

中國(guó)地震局地質(zhì)研究所 中國(guó)地震局活動(dòng)構(gòu)造與火山重點(diǎn)實(shí)驗(yàn)室 吳熙彥

活動(dòng)斷層是確定區(qū)域內(nèi)最大發(fā)震震級(jí)和破壞強(qiáng)度的重要資料。隨著計(jì)算機(jī)技術(shù)的發(fā)展，在國(guó)家支持下，地震活動(dòng)斷層調(diào)查相關(guān)項(xiàng)目也獲得了大量活動(dòng)斷層調(diào)查空間數(shù)據(jù)。將這些數(shù)據(jù)進(jìn)行綜合，提取關(guān)鍵信息，形成一套格式、結(jié)構(gòu)統(tǒng)一的數(shù)據(jù)，有利于提高效率和進(jìn)一步的數(shù)據(jù)分析。

矢量數(shù)據(jù)融合；活動(dòng)斷層；數(shù)據(jù)；pyshp

一、全國(guó)活動(dòng)斷層數(shù)據(jù)融合意義

中國(guó)位于世界兩大地震帶——環(huán)太平洋地震帶與歐亞地震帶之間，地震頻發(fā)，造成的災(zāi)害損失巨大。然而，地震短臨預(yù)報(bào)仍是世界難題，短期內(nèi)無(wú)法攻克。目前，國(guó)際和國(guó)內(nèi)防震減災(zāi)的一項(xiàng)主要措施是通過(guò)確定某區(qū)域可能發(fā)生地震時(shí)造成的破壞強(qiáng)度，設(shè)定抗震設(shè)防烈度，規(guī)定建構(gòu)筑物要達(dá)到的抗震要求，減少地震災(zāi)害造成的損失?；顒?dòng)斷層則是確定區(qū)域內(nèi)最大發(fā)震震級(jí)和破壞強(qiáng)度的重要資料。隨著計(jì)算機(jī)技術(shù)的發(fā)展，在國(guó)家支持下，地震活動(dòng)斷層調(diào)查相關(guān)項(xiàng)目也獲得了大量活動(dòng)斷層調(diào)查空間數(shù)據(jù)。這些數(shù)據(jù)獲取的時(shí)間跨度大，來(lái)源不一。因此，數(shù)據(jù)的格式、結(jié)構(gòu)和制作標(biāo)準(zhǔn)等并不一致。使用時(shí)需要調(diào)用和查看不同的數(shù)據(jù)集，調(diào)用不同的圖層，效率較低。同類數(shù)據(jù)分層分包存儲(chǔ)，也不利于空間查詢、統(tǒng)計(jì)等分析研究。將這些數(shù)據(jù)進(jìn)行綜合，提取關(guān)鍵信息，形成一套格式、結(jié)構(gòu)統(tǒng)一的數(shù)據(jù)，有利于提高效率和進(jìn)一步的數(shù)據(jù)分析。

二、全國(guó)活動(dòng)斷層數(shù)據(jù)融合流程設(shè)計(jì)

不同時(shí)期、不同來(lái)源的活動(dòng)斷層GIS數(shù)據(jù)融合的內(nèi)容包括空間圖形和屬性值[1]。將數(shù)據(jù)融合的思路分解為圖形的融合和屬性值的融合，以先圖形后屬性的方式進(jìn)行流程設(shè)計(jì)。數(shù)據(jù)融合的步驟如下。

（1）轉(zhuǎn)換文件格式：原有數(shù)據(jù)采用了不同的文件格式。數(shù)據(jù)處理的第一步，是將所有數(shù)據(jù)統(tǒng)一轉(zhuǎn)換為shape file文件。Shapefile文件是被廣泛接受的一種GIS數(shù)據(jù)文件。常用的GIS軟件都能將GIS數(shù)據(jù)轉(zhuǎn)換為這種格式。

（2）轉(zhuǎn)換空間坐標(biāo)系：不同來(lái)源的GIS數(shù)據(jù)可能在不同的空間坐標(biāo)系下，應(yīng)在處理之前轉(zhuǎn)換到統(tǒng)一的坐標(biāo)系下。目前，我國(guó)最新的坐標(biāo)系統(tǒng)為2000國(guó)家大地坐標(biāo)系，它是國(guó)家基礎(chǔ)地理信息數(shù)據(jù)（基礎(chǔ)地理底圖）所采用的坐標(biāo)系，由國(guó)務(wù)院批準(zhǔn)自2008年7月1日啟用。為推廣2000國(guó)家大地坐標(biāo)系的使用，國(guó)家測(cè)繪地理信息局公布了現(xiàn)有測(cè)繪成果轉(zhuǎn)換到該坐標(biāo)系的技術(shù)指南。依據(jù)技術(shù)指南，可以將現(xiàn)存數(shù)據(jù)都轉(zhuǎn)換到2000國(guó)家大地坐標(biāo)系。數(shù)據(jù)轉(zhuǎn)換到統(tǒng)一坐標(biāo)系后，可實(shí)現(xiàn)多來(lái)源空間數(shù)據(jù)疊加，進(jìn)行編輯和修改。

（3）分區(qū)處理數(shù)據(jù)：空間數(shù)據(jù)疊加后，可以看到部分區(qū)域內(nèi)只有一層圖層，每個(gè)對(duì)象只有一種圖形表達(dá)，而部分區(qū)域有一層以上的重疊數(shù)據(jù)。對(duì)于只有唯一空間數(shù)據(jù)的區(qū)域，不需進(jìn)行空間操作，直接保留。對(duì)于重疊區(qū)域，須進(jìn)行圖形合并。因此，在圖形合并前，需將將數(shù)據(jù)集切塊，切出重疊區(qū)的部分待下一步處理。

（4）重疊區(qū)圖形數(shù)據(jù)合并：地圖合并技術(shù)一般包括匹配實(shí)體和將匹配的同名實(shí)體合并兩個(gè)過(guò)程[2][3]。實(shí)體自動(dòng)匹配主要依靠幾何相似度、拓?fù)?、?shù)據(jù)模型和屬性等。同名實(shí)體自動(dòng)合并主要根據(jù)預(yù)設(shè)權(quán)重和待合并圖形幾何形態(tài)生成新的圖形。然而，以線要素表達(dá)的活動(dòng)斷層線的幾何形態(tài)并非地質(zhì)對(duì)象的真實(shí)形態(tài)，即地殼斷裂帶的真實(shí)幾何形態(tài)，而是地質(zhì)學(xué)家對(duì)地質(zhì)對(duì)象的一種解釋和抽象。同一實(shí)體的不同解釋和抽象可能使幾何形態(tài)相似程度較低。因此，活動(dòng)斷層的幾何形態(tài)的，并不具備很好地自動(dòng)匹配和自動(dòng)合并的條件。因此，采用以下三步合并數(shù)據(jù)。1）選取重疊區(qū)某一圖層為最終保留層（稱基準(zhǔn)圖層），用人機(jī)交互的方式修改該圖層；2）為每一圖層增加“同名實(shí)體編號(hào)”字段。所有同名的實(shí)體賦予唯一的一個(gè)編號(hào)填入該字段；3）利用“同名實(shí)體編號(hào)”連接其他數(shù)據(jù)屬性表到基準(zhǔn)圖層，將同名實(shí)體的屬性全部合并入基準(zhǔn)圖層。最終重疊區(qū)僅保留基準(zhǔn)圖層一層數(shù)據(jù)。

（5）確定屬性表結(jié)構(gòu)：合并重疊區(qū)和非重疊區(qū)的空間數(shù)據(jù)，形成一層新的圖層。新數(shù)據(jù)的屬性表是原有數(shù)據(jù)所有屬性的集合。新屬性表存在重復(fù)信息，需要重新設(shè)計(jì)與規(guī)劃。根據(jù)屬性內(nèi)容，重新設(shè)計(jì)屬性表結(jié)構(gòu)，在不損失信息條件下精簡(jiǎn)表結(jié)構(gòu)。

（6）屬性融合：根據(jù)重新設(shè)計(jì)的表結(jié)構(gòu)，為屬性表新增字段，準(zhǔn)備填入信息融合后的屬性值。屬性處理有3種情況。1）屬性值直接賦給新表中對(duì)應(yīng)屬性字段；2）多個(gè)信息重復(fù)的屬性值，合并到新表中的一個(gè)屬性字段；3）多個(gè)信息不同的屬性值，合并到新表中的一個(gè)屬性字段。根據(jù)新形成的數(shù)據(jù)和新設(shè)計(jì)的表結(jié)構(gòu)，設(shè)計(jì)屬性遷移和合并規(guī)則。

（7）刪除冗余字段：將上一步形成的新數(shù)據(jù)中的冗余字段刪除，獲得最終成果數(shù)據(jù)。

圖1 流程圖

三、方案可行性分析

1.pyshp介紹

方案的第1-5步都較容易實(shí)現(xiàn)。方案實(shí)現(xiàn)的關(guān)鍵是第6步屬性融合。筆者采用了pyshp作為實(shí)現(xiàn)第6步的試驗(yàn)工具。pyshp全稱Python Shapefile Library，是基于python語(yǔ)言的開(kāi)源函數(shù)庫(kù)，提供了讀寫(xiě)ESRI shapefile格式文件的方法。pyshp函數(shù)庫(kù)具有編譯工具簡(jiǎn)單、免安裝、代碼開(kāi)源、使用簡(jiǎn)單的特點(diǎn)，并且提供了簡(jiǎn)單易懂的使用說(shuō)明和示例代碼，在ESRI shapefile編輯方面便利，非常適合地理數(shù)據(jù)融合工具。pyshp提供三個(gè)基礎(chǔ)類:Reader、Writer、Editor。Reader 提供了ESRI shapefile的文件頭信息、幾何圖形和屬性表的讀取方法；Writer提供了ESRI shapefile的文件頭信息、幾何圖形和屬性表的寫(xiě)入方法；Editor是Writer子類，繼承了Writer的所有方法，主要用于已有文件的編輯。pyshp可以在任何pyshon平臺(tái)直接使用import命令導(dǎo)入模塊文件調(diào)用函數(shù)。本文使用的是pyshp 1.2.3。

2.關(guān)鍵步驟的實(shí)現(xiàn)

這里僅討論關(guān)鍵的兩種屬性處理。

（1）多個(gè)信息重復(fù)的屬性值，合并到新表中的一個(gè)屬性字段。例如，斷層活動(dòng)年代的編碼方式不同，一部分代碼值為{a1,a2,a3,...ai}，另一部分代碼值為{b1,b2,b3,...ai}。實(shí)際上，a1和b1代表相同屬性。該屬性的新代碼值為{c1,c2,c3,...ci}。此時(shí)，設(shè)計(jì)一個(gè)代碼映射表，將代碼一一映射到新代碼集。核心python代碼如下：

iSetField=self.getEsfFieldIndex(esf, SetField)#新字段

iReferField=self.getEsfFieldIndex(esf, ReferField)#舊代碼值所在字段nRecord=0

for xRecord in esf.records:

#如果新代碼字段是空的，就將舊代碼值的對(duì)應(yīng)新代碼填入

if xRecord[iSetField] is None or xRecord[iSetField]==0:

reSValue=str(u(xRecord[iReferField]))

xRecord[iSetField]=int(AgeDic.get(reSValue,0))#從代碼映射表獲得新代碼

nRecord+=1

return nRecord

（2）多個(gè)信息不同的屬性值，合并到新表中的一個(gè)屬性字段。例如，原數(shù)據(jù)有一個(gè)表示斷層是“出露”還是“隱伏”的字段，還有另一個(gè)表示“正斷”還是“逆斷”的字段。兩種信息可以合并，成為新的斷層性質(zhì)字段。新編代碼的前幾位代表前一個(gè)字段值，后幾位代表后一個(gè)字段值。核心python代碼如下：

for xRecord in esf.records:

Sfirvalue=str(u(xRecord[iFirField])).strip() #獲得第一個(gè)字段值

xFir=NewCodeFirDic.get(Sfirvalue,0) #從代碼映射表獲得第一個(gè)字段值的新代碼

SSecvalue=str(u(xRecord[iSecField])).strip() #獲得第二個(gè)字段值

xSec=NewCodeFirDic.get(SSecvalue,0) #從代碼映射表獲得第二個(gè)字段值的新代碼

xRecord[iSetField]= xFir.join(xSec)

nRecord+=1

return nRecord

[1]唐文靜．多源地理空間矢量數(shù)據(jù)融合[M]．清華大學(xué)出版社,2014．

[2]張橋平,李德仁,龔健雅．地圖合并技術(shù)[J]．測(cè)繪通報(bào),2001(1)：6-8．

[3]COBB M,CHUNG M,FOLEY H．A Rule-based Approach for the Conflation of Attributed Vector Data．GeoInformatica,1998,2(1)：7-35．

中國(guó)大陸主要地震構(gòu)造帶活動(dòng)斷層探察專項(xiàng)（IGCEA1621）子課題“郯廬斷裂帶三維地震構(gòu)造分段模型與地震預(yù)測(cè)研究——數(shù)據(jù)庫(kù)建設(shè)”；中國(guó)地震局地震行業(yè)科研專項(xiàng)（20138001、201408023、200908001、201008001）。