Python語言在物探數(shù)據(jù)清洗中的探索

謝明宏， 邱崇濤， 祁 程， 張 偉， 何昕欣

(核工業(yè)航測(cè)遙感中心 鈾資源地球物理勘查技術(shù)中心重點(diǎn)實(shí)驗(yàn)室，石家莊 050002)

0 前言

隨著物探技術(shù)、儀器的不斷改進(jìn)和計(jì)算技術(shù)高速發(fā)展，多方法多信息的結(jié)合將成為今后物探工作的一個(gè)趨勢(shì)，地質(zhì)、物探工作者將面臨著大數(shù)據(jù)時(shí)代的到來。而數(shù)據(jù)處理是物探工作中最重要的環(huán)節(jié)之一，主要包括預(yù)處理、數(shù)據(jù)修正、反演計(jì)算、網(wǎng)格化以及數(shù)據(jù)可視化等一系列環(huán)節(jié)，處理成果的好壞、執(zhí)行效率直接決定著后續(xù)的資料解釋工作。數(shù)據(jù)準(zhǔn)備是數(shù)據(jù)預(yù)處理的第一步，在物探工作中，常常要面臨新數(shù)據(jù)結(jié)構(gòu)或要解決新問題，此項(xiàng)工作往往將會(huì)反復(fù)多次。依照以往經(jīng)驗(yàn)，在數(shù)據(jù)分析中，80%的時(shí)間將花費(fèi)在數(shù)據(jù)清理或數(shù)據(jù)準(zhǔn)備等工作中[1]。雖然數(shù)據(jù)準(zhǔn)備(數(shù)據(jù)清理)將會(huì)花費(fèi)了大量的時(shí)間和精力，但是對(duì)于此方面的相關(guān)研究并不多。Hadely Wickham[2]提出了潔凈數(shù)據(jù)(Tidy Data)的概念，筆者將以此為基礎(chǔ)上，結(jié)合常見的物探數(shù)據(jù)格式、數(shù)據(jù)預(yù)處理方法，總結(jié)了不潔凈數(shù)據(jù)類型，并提出了相應(yīng)的補(bǔ)救措施，利用當(dāng)今最為流行的科學(xué)計(jì)算語言Python，嘗試提出數(shù)據(jù)預(yù)處理中數(shù)據(jù)清洗、假值剔除和數(shù)據(jù)內(nèi)插的一些常規(guī)方法，構(gòu)建易于數(shù)據(jù)分析、參量提取和數(shù)據(jù)可視化等數(shù)據(jù)結(jié)構(gòu)或數(shù)據(jù)集，使工作更為輕松、高效。

1 物探數(shù)據(jù)格式分析

1.1 潔凈數(shù)據(jù)的屬性

事實(shí)上，潔凈數(shù)據(jù)為方便計(jì)算機(jī)程序而提出的一種思想，未有十分明確的定義，但需包括如下三個(gè)屬性[2]：①每個(gè)變量構(gòu)成一列，包含有變量名稱和實(shí)測(cè)值(如“電阻率”變量名為“Res”，測(cè)量值為213.26 Ω·m)；②每個(gè)觀測(cè)值構(gòu)成一行(如測(cè)線統(tǒng)計(jì)時(shí)，每一行觀測(cè)值都應(yīng)包含測(cè)線編號(hào)，而不能進(jìn)行將其合并)；③同類型觀測(cè)值構(gòu)成一個(gè)數(shù)據(jù)集或一個(gè)數(shù)據(jù)單元塊。此處的同類型數(shù)據(jù)不同測(cè)區(qū)和測(cè)線，或者是不同步的測(cè)量值。例如航磁測(cè)量數(shù)據(jù)采集間隔是10 次/s，而航放是1 次/s(在多方法測(cè)量時(shí)，因異常追索，加密測(cè)點(diǎn))。

表1 測(cè)線飛行高度質(zhì)量統(tǒng)計(jì)報(bào)表

1.2 常見不潔凈數(shù)據(jù)

物探數(shù)據(jù)大多由行、列構(gòu)成的二維矩形數(shù)據(jù)表，但實(shí)測(cè)或整理的數(shù)據(jù)(集)常出現(xiàn)不滿足潔凈數(shù)據(jù)三個(gè)準(zhǔn)則的情況，此處可稱作不潔凈數(shù)據(jù)集(Messy Dataset)。處理工作中，數(shù)據(jù)(統(tǒng)計(jì))表格設(shè)計(jì)常以閱讀角度呈現(xiàn)，原始數(shù)據(jù)(集)依據(jù)采集時(shí)間順序排列，或者是按照參量分類給出；正因物探數(shù)據(jù)格式多樣性和特殊性，并不完全適合計(jì)算機(jī)程序讀取的標(biāo)準(zhǔn)數(shù)據(jù)結(jié)構(gòu)，所以在數(shù)據(jù)讀取、統(tǒng)計(jì)、建模、可視化時(shí)，則需要更多的并缺乏復(fù)用性的代碼，導(dǎo)致工作效率降低。在物探測(cè)量中，常見不潔凈數(shù)據(jù)包括以下幾種類型：

1.2.1 數(shù)值為列標(biāo)頭而非變量名稱

表1為某一架次航空物探飛行高度質(zhì)量統(tǒng)計(jì)報(bào)表[3]。此表格從閱讀角度上，設(shè)計(jì)較為合理。但飛行高度(范圍)作為表頭(變量)給出，程序?qū)⑹褂酶啻a進(jìn)行讀取。

此時(shí)，將表1中的飛行高度作為變量，表2的數(shù)據(jù)結(jié)構(gòu)更適合于程序閱讀。

1.2.2 一列中存儲(chǔ)多個(gè)變量

物探數(shù)據(jù)中常出現(xiàn)諸如經(jīng)緯度(XY坐標(biāo))、測(cè)量值的變化范圍等復(fù)合參數(shù)。表3中的電阻率范圍[4]存儲(chǔ)在同一列中。此時(shí)應(yīng)將“變化范圍”分割為“最小值”和“最大值”兩列(表4)。

表2 測(cè)線飛行高度統(tǒng)計(jì)報(bào)表

表3 相山地區(qū)主要巖石電阻率統(tǒng)計(jì)一覽表

表4 相山地區(qū)主要巖石電阻率統(tǒng)計(jì)一覽表

1.2.3 變量存儲(chǔ)在行、列中

表5顯示了某地區(qū)航空測(cè)量設(shè)計(jì)測(cè)線。坐標(biāo)值(端點(diǎn))和線號(hào)(端點(diǎn))變量分別處于行、列中。為了計(jì)算方便，應(yīng)將同條測(cè)線的所有變量融合在同一行中(表6)。

表5 黑龍江省某地區(qū)航空物探測(cè)量設(shè)計(jì)測(cè)線一覽表(部分)

表6 黑龍江省某地區(qū)航空物探測(cè)量設(shè)計(jì)測(cè)線一覽表(部分)

1.2.4 一個(gè)數(shù)據(jù)集(表)中存儲(chǔ)多個(gè)類型

此類型常在電(磁)法測(cè)量數(shù)據(jù)中出現(xiàn)(如地球物理勘探協(xié)會(huì)(SEG)推薦的EDI格式、Zone 公司的AVG等)。數(shù)據(jù)結(jié)構(gòu)是一個(gè)測(cè)站(排列)為記錄單元，由多個(gè)數(shù)據(jù)塊(block)構(gòu)成，而各數(shù)據(jù)塊以關(guān)鍵字開頭，包括站點(diǎn)位置、電、磁道信息、接收頻率、視電阻率、相位、離差等一系列參數(shù)。此時(shí)，可根據(jù)數(shù)據(jù)塊的關(guān)鍵字(變量)進(jìn)行分離、組合，形成易于讀取、分析、可視化等的數(shù)據(jù)結(jié)構(gòu)。

1.2.5 同類型數(shù)據(jù)處于多個(gè)數(shù)據(jù)集中

物探測(cè)量工作一般分為剖面測(cè)量(如電法)和面積性測(cè)量(如磁法)。測(cè)量數(shù)據(jù)以測(cè)點(diǎn)為基本單元，多個(gè)測(cè)點(diǎn)構(gòu)成測(cè)線(剖面線)數(shù)據(jù)，再由多個(gè)測(cè)線(剖面線)形成一定網(wǎng)度的測(cè)量數(shù)據(jù)。所以數(shù)據(jù)清洗應(yīng)根據(jù)測(cè)量性質(zhì)來決定數(shù)據(jù)結(jié)構(gòu)，單剖面抑或全區(qū)數(shù)據(jù)。

另外，由于測(cè)量手段限制，觀測(cè)值有時(shí)會(huì)存儲(chǔ)在兩個(gè)或多個(gè)數(shù)據(jù)集中，最為典型的是GPS定位數(shù)據(jù)與地面物探實(shí)測(cè)數(shù)據(jù)處在不同數(shù)據(jù)文件中。

上述問題常用解決辦法包括：①新建一個(gè)數(shù)據(jù)列表，關(guān)聯(lián)各個(gè)數(shù)據(jù)集；②增加新變量(數(shù)據(jù))列；③將所有數(shù)據(jù)集合并在一個(gè)數(shù)據(jù)集中。

1.2.6 數(shù)據(jù)集變量存有差異

物探測(cè)量中，同種儀器因采用不同的測(cè)量模式、測(cè)量裝置、采樣時(shí)間等，獲取的結(jié)果數(shù)據(jù)也不同(如大地電磁測(cè)量模式可分為標(biāo)量和張量?jī)煞N)；激電測(cè)量(IPR-12多功能電法儀)會(huì)因采樣時(shí)間長(zhǎng)短產(chǎn)生不同的數(shù)據(jù)格式。此時(shí)，視具體測(cè)量模式，采用枚舉法讀取數(shù)據(jù)。

1.2.7 丟碼或誤碼現(xiàn)象

由于設(shè)備和數(shù)據(jù)處理方法的原因，物探數(shù)據(jù)中會(huì)出現(xiàn)丟碼或誤碼現(xiàn)象(如航空物探測(cè)量和磁日變采集等數(shù)據(jù)連續(xù)采集中產(chǎn)生的誤碼或漏碼)，數(shù)據(jù)網(wǎng)格化中的不合理值或空值，數(shù)據(jù)采集中的跳變值、畸變值等。誤碼、漏碼雖在數(shù)量或出現(xiàn)頻率并不多，但也絕不容忽視。須確定誤碼或漏碼具體位置和問題緣由后，視具體情況，進(jìn)行數(shù)據(jù)標(biāo)示、剔除或內(nèi)插。

從上述7種現(xiàn)象可看出，不潔凈數(shù)據(jù)①源于數(shù)據(jù)質(zhì)量;②數(shù)據(jù)結(jié)構(gòu)。但隨物探儀器的性能、處理程序的智能化不斷提升，數(shù)據(jù)質(zhì)量隨之提高、數(shù)據(jù)結(jié)構(gòu)也趨于規(guī)范化，所以物探數(shù)據(jù)中的不潔凈主要源于數(shù)據(jù)結(jié)構(gòu)。

從某種意義上講，數(shù)據(jù)清理是一個(gè)數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換過程。在數(shù)據(jù)預(yù)處理中，應(yīng)首先清洗數(shù)據(jù)本身的假值、空值或畸變值等，然后再進(jìn)行數(shù)據(jù)結(jié)構(gòu)的轉(zhuǎn)換。最終形成適合于計(jì)算機(jī)程序運(yùn)算的數(shù)據(jù)結(jié)構(gòu)，即潔凈數(shù)據(jù)。

對(duì)于EDI、AVG格式等大地電磁法中的通用數(shù)據(jù)結(jié)構(gòu)，事實(shí)上在指定程序執(zhí)行反演計(jì)算時(shí)，也需要一個(gè)格式轉(zhuǎn)換的過程；同樣在進(jìn)行其他統(tǒng)計(jì)、分析，仍要依照潔凈數(shù)據(jù)結(jié)構(gòu)提取、整理相關(guān)參量。

2 物探數(shù)據(jù)清洗方法

2.1 數(shù)據(jù)清洗中常用的Python第三方庫(kù)

Python是一種面向?qū)ο蟮慕忉屝缘目缙脚_(tái)高級(jí)計(jì)算機(jī)語言，具有語法簡(jiǎn)單、可移植性強(qiáng)的特點(diǎn)，同時(shí)，具有極為豐富的類庫(kù)，廣泛應(yīng)用于數(shù)據(jù)庫(kù)、多媒體、科學(xué)計(jì)算、網(wǎng)絡(luò)等諸多領(lǐng)域。正因Python具有強(qiáng)大的科學(xué)計(jì)算、繪圖功能，以及在地圖繪制、地理空間數(shù)據(jù)的處理與轉(zhuǎn)換方面具有豐富的類庫(kù)[5]，國(guó)外許多著名軟件(如arcGIS地理信息系統(tǒng)、Montaj Oasis地球物理處理系統(tǒng)等)以它作為主要腳本語言。

在數(shù)據(jù)清洗中，常用的Python第三方庫(kù)主要包括以下四個(gè)：

1)Pandas：基于NumPy構(gòu)建的數(shù)據(jù)分析庫(kù)[6,7]，數(shù)據(jù)格式有Series(一維)和DataFrame(二維)兩種。它兼具了NumPy高性能的數(shù)組計(jì)算及電子表格和關(guān)系型數(shù)據(jù)庫(kù)靈活的數(shù)據(jù)處理功能，提供了復(fù)雜精細(xì)的索引功能更便捷地完成重塑、切片和切塊、聚合及選取數(shù)據(jù)子集等操作[6]。數(shù)據(jù)清理中，Pandas將承擔(dān)著最為重要的角色，包括數(shù)據(jù)讀取、保存，數(shù)據(jù)分離與合并、排序、統(tǒng)計(jì)等。

2)NumPy：運(yùn)算速度非?？斓臄U(kuò)展數(shù)學(xué)庫(kù)，可支持N維數(shù)組和矩陣計(jì)算并提供了大量的數(shù)學(xué)函數(shù)庫(kù)，具有強(qiáng)大的廣播功能[8-9]。與Pandas交互使用，可有效地提高運(yùn)算速度。

3)Scipy：功能全面的算法庫(kù)和數(shù)學(xué)工具包，包括線性代數(shù)、積分、插值、快速傅里葉變換、信號(hào)處理和圖像處理等[8]?；贜umpy上，充分利用Scipy提供的函數(shù)，可優(yōu)化程序設(shè)計(jì)，縮短程序開發(fā)周期。

4)Matplotlib：基于NumPy數(shù)組功能最全、最為經(jīng)典的繪圖庫(kù)[9]，與NumPy和Scipy模塊組合完全可與Matlib語言媲美。

在地質(zhì)、物探、GIS等方面，常用的Python第三方庫(kù)還包括有Shapely(二維圖形分析)、GDAL(地理空間數(shù)據(jù)轉(zhuǎn)換庫(kù))、pyProj(地圖投影換算)、cartopy(地圖繪制)、geoPandas(地圖繪制)、Fiona(地理空間數(shù)據(jù)的讀寫)，pyshp(ESRI Shapefile 文件讀寫)、pyqtGraph(可融入界面的科學(xué)圖形繪制庫(kù))和Bokeh(交互式網(wǎng)頁(yè)圖形繪制)等。

總之，Python在數(shù)據(jù)清理或預(yù)處理中，優(yōu)勢(shì)明顯，可大大縮短開發(fā)周期。

2.2 常規(guī)數(shù)據(jù)清洗操作

在物探數(shù)據(jù)清洗中，常用數(shù)據(jù)操作一般包括變量轉(zhuǎn)換、篩選、聚合、排序、插值等。

1)變量轉(zhuǎn)換：增加或修改變量。變量修改可劃分為單變量修改和多變量修改兩種。單變量修改如電法測(cè)量中的電磁波頻率轉(zhuǎn)換為對(duì)數(shù)格式，多變量修改如利用實(shí)測(cè)電阻和截面積計(jì)算電阻率。對(duì)于此類常規(guī)計(jì)算，Numpy中的通用函數(shù)(ufun)完全可以滿足,再配合廣播(Broadcasting)功能[8-9]將大大減少代碼輸入量，提高運(yùn)算速度。

2)篩選：基于某個(gè)或多個(gè)條件，對(duì)數(shù)據(jù)集進(jìn)行分組和移除變量，如測(cè)線分離，數(shù)據(jù)的分類統(tǒng)計(jì)等,主要由Pandas完成。

3)聚合：多個(gè)測(cè)量值合并為一個(gè)值。如地面磁、放測(cè)量，最終值取多次測(cè)量值的均值；放射性能譜測(cè)量時(shí)，各道計(jì)數(shù)率取某時(shí)間段累加值，主要由Pandas完成。

4)排序：將觀測(cè)值依照某一個(gè)條件或多個(gè)條件排序，Pandas或Numpy均提供相應(yīng)函數(shù)，視具體情況選擇。

5)插值：補(bǔ)充丟失或未獲取的測(cè)量值，如日變修正時(shí)日變數(shù)據(jù)的內(nèi)插；Scipy提供線性(linear)、二次(quadratic)三次(cubic)、樣條(spline)等等諸多插值方法[8]。此處多指一維插值。對(duì)于二維插值(數(shù)據(jù)網(wǎng)格化)，建議使用專業(yè)Geosoft或Surfer軟件提供的網(wǎng)格化方法。

3 物探標(biāo)準(zhǔn)數(shù)據(jù)格式的建立

在潔凈數(shù)據(jù)和不潔凈數(shù)據(jù)屬性的基礎(chǔ)上，結(jié)合物探數(shù)據(jù)特點(diǎn)，對(duì)標(biāo)準(zhǔn)數(shù)據(jù)格式的建立提出如下幾點(diǎn)建議。

3.1 數(shù)據(jù)結(jié)構(gòu)的構(gòu)建

3.1.1 參據(jù)名稱

數(shù)據(jù)集首行應(yīng)標(biāo)注參量名稱，增加可讀性，避免參量混淆。參量名稱標(biāo)明富有明確的含義，例如，經(jīng)緯度可用Lon和Lat注明，坐標(biāo)值用X和Y，必要時(shí)使用X84和Y84 或X54和Y54，表明不同坐標(biāo)系投影坐標(biāo)。再如，原始磁數(shù)據(jù)為rawMag，日變修正后數(shù)據(jù)為MagDiurnalCorr。參量名稱中不能使用空格、中文或特殊字符，避免程序讀取錯(cuò)誤。

3.1.2 CSV格式

CSV格式是指以逗號(hào)為分隔符的ASCII文本文件，尤其適合數(shù)據(jù)量小的數(shù)據(jù)集。其優(yōu)點(diǎn)可借助Excel軟件，易于數(shù)據(jù)閱讀，便于簡(jiǎn)單運(yùn)算、數(shù)據(jù)列的增減。大多數(shù)知名軟件可直接讀取。Python(Pandas)讀取此類型文件也極為方便，僅需一條語句，df = pd.read.csv (r'./mydata.csv')，并能方便地進(jìn)行行、列的操作(分類、統(tǒng)計(jì)、計(jì)算等)。

3.1.3 HDF5格式

HDF5(層次性數(shù)據(jù)格式)是用于存儲(chǔ)大數(shù)據(jù)量、單一數(shù)字型數(shù)組的機(jī)制，可利用層次化組織結(jié)構(gòu)和標(biāo)記含有任意元數(shù)據(jù)的數(shù)據(jù)集建立數(shù)據(jù)模型[10]。它有兩種基本對(duì)象：組(Group)和數(shù)據(jù)集(Dataset)。數(shù)據(jù)集的數(shù)據(jù)結(jié)構(gòu)屬于典型的潔凈數(shù)據(jù)，而組可視為一容器，包含組(次一級(jí))和各數(shù)據(jù)集?？梢哉f，HDF5的層次性機(jī)制非常適合于數(shù)據(jù)量大、結(jié)構(gòu)復(fù)雜物探數(shù)據(jù)的存儲(chǔ)和讀寫。在Python中，h5py(第三方庫(kù))提供HDF5格式文件的接口，并提供了數(shù)據(jù)集建立、讀取、并行計(jì)算等一系列功能。例如：表7顯示了256道航空伽瑪能譜測(cè)量中用于數(shù)據(jù)質(zhì)量控制的早、晚基線測(cè)試數(shù)據(jù)。

表7中的256道能譜數(shù)據(jù)共包含534列，轉(zhuǎn)換成HDF5結(jié)構(gòu)如圖1所示。從圖1中可看出，數(shù)據(jù)結(jié)構(gòu)層次清晰，所以計(jì)算飛行前、后NaI晶體上測(cè)、下測(cè)分辨率和峰漂時(shí)，參量讀取(包括切片)、運(yùn)算等也將極為方便。若利用256道能譜數(shù)據(jù)(測(cè)線)進(jìn)行全譜分析，速度優(yōu)勢(shì)將更為明顯。

3.2 數(shù)據(jù)集的管理

一般情況下，一個(gè)物探測(cè)區(qū)由多條測(cè)線(剖面)或測(cè)網(wǎng)構(gòu)成。對(duì)于面積型測(cè)量，應(yīng)將所有測(cè)線合并為一個(gè)文件，易于數(shù)據(jù)網(wǎng)格化；單測(cè)線分析或繪制剖面平面圖時(shí)，再進(jìn)行測(cè)線分離。對(duì)于剖面測(cè)量，剖面數(shù)據(jù)應(yīng)分別存儲(chǔ)，確保格式的統(tǒng)一，也可通過文件列表進(jìn)行批量處理。

表7 某地區(qū)航放早晚基線測(cè)試256道譜數(shù)據(jù)(部分)

圖1 HDF5結(jié)構(gòu)示例(256道能譜數(shù)據(jù))Fig.1 1 Structure of HDF5 (256 channel spectral data)

4 Python在數(shù)據(jù)清洗中的應(yīng)用

4.1 數(shù)據(jù)結(jié)構(gòu)的清洗

此例是將不潔凈數(shù)據(jù)(表1)轉(zhuǎn)換為潔凈數(shù)據(jù)(表2)，展示Pandas庫(kù)強(qiáng)大的文件讀寫和數(shù)據(jù)結(jié)構(gòu)合并等功能。Python程序代碼示例如下：

import pandas as pd

def meltData(df, cols, myValue):

df.columns = cols

dfTrans = pd.melt(df, ['LineNo'], var_name='FlightHeight', value_name=myValue)

return dfTrans

cols = ['LineNo', 'X<80', '80<=X<100', '100<=X<120', '120<=X<150', '>=150']

#讀取6列數(shù)據(jù)

df1 = pd.read_csv('stat.csv', usecols=[0, 1, 3, 5, 7, 9], skiprows=2, header=None)

dfTrans1 = meltData(df1, cols, myValue='Amount') #轉(zhuǎn)換格式

df2 = pd.read_csv(inFile, usecols=[0, 2, 4, 6, 8, 10], skiprows=2, header=None)

dfTrans2 = meltData(df2, cols, myValue='Proportion') #轉(zhuǎn)換格式

dfTrans1['Proportion'] = dfTrans2['Proportion']#將dfTrans2的'Proportion'加入dfTrans1中

dfTrans1.to_csv('stat_1.csv', index=False)#保存數(shù)據(jù)(CSV格式)

4.2 航測(cè)數(shù)據(jù)清洗和質(zhì)量評(píng)價(jià)

對(duì)于航空物探，偏航距或離地高度是衡量航放測(cè)量飛行質(zhì)量的兩個(gè)重要指標(biāo)。下面以黑龍江省完達(dá)山—太平嶺地區(qū)航測(cè)項(xiàng)目中某一個(gè)架次的航測(cè)數(shù)據(jù)為例，展示使用Pandas庫(kù)進(jìn)行測(cè)線分離后的分類統(tǒng)計(jì)。

1)由Pico公司生產(chǎn)的綜合航空物探測(cè)量系統(tǒng)采集的航測(cè)能譜數(shù)據(jù)(部分)：

------能譜窗寬------

D_TC: 34 239

D_K: 116 133

D_U: 141 158

D_Th: 205 239

D_U1: 88 101

UP_TC: 34 239

UP_U: 141 158

USE1: 117 134

USE2: 142 159

USE3: 206 240

GPS坐標(biāo) WGS84坐標(biāo)系 GAUSS 22N投影 CM=129度

雷達(dá)高度RALT：米，氣壓高度BARO：米，溫度TP：度 Spc1: ISPD Spc2: ISPU

ACFT: Y12(B-3833) PN: wandashan FLTN: 41 OpBy:Yangjinzheng Task: Sec Survey Date: 09/01/14

LINE FN Gps_L Gps_B X Y GPS_Alt RALT BARO TP Time Date

6190 001110 131.2345770 46.9158350 670216.87 5200315.25 165.727 103.36 154.99 22.4 ... 06:05:53.000 09/01/14

6190 001120 131.2351360 46.9153710 670260.94 5200264.87 163.330 96.26 154.15 22.4 ... 06:05:54.000 09/01/14

6190 001130 131.2356740 46.9148930 670303.47 5200212.99 161.216 90.17 151.92 22.4 ... 06:05:55.000 09/01/14

6190 001140 131.2362020 46.9144090 670345.23 5200160.29 159.698 90.89 148.30 22.4 ... 06:05:56.000 09/01/14

..... ...

6200 120600 131.2429890 46.9151710 670859.79 5200259.71 151.597 80.96 137.72 25.4 ... 09:27:59.000 09/01/14

6200 120610 131.2423530 46.9156170 670809.86 5200307.93 150.288 78.65 134.94 25.4 ... 09:28:00.000 09/01/14

6200 120620 131.2417150 46.9160630 670759.90 5200356.09 149.569 76.57 133.83 25.4 ... 09:28:01.000 09/01/14

6200 120630 131.2410730 46.9165030 670709.56 5200403.71 150.372 76.84 131.89 25.4 ... 09:28:02.000 09/01/14

2)Python程序代碼示例

import numpy as np

import pandas as pd

heightClass = [0, 80, 100, 120, 150] #飛行高度分類級(jí)別

data_file = r' cq_B4090106_C_cut_Spc.dat'

#跳過15行，讀取數(shù)據(jù)中0, 4,5,7共4列數(shù)據(jù)，

df = pd.read_table(ffile, delimiter='s+', skiprows=15, usecols=(0, 4, 5, 7))

print (df.head(5)) #顯示前5行

#斜字體為輸出結(jié)果，下同。

LINEXYRALT

06190670216.875200315.25103.36

16190670260.945200264.8796.26

26190670303.475200212.9990.17

36190670345.235200160.2990.89

46190670386.835200107.4186.48

#通過設(shè)計(jì)測(cè)線，計(jì)算飛行偏離度

disErrSet = statPathDev (df, designLines)

#df數(shù)據(jù)集中增加飛行偏離度列(disErr)

df['disErr'] = disErrSet

LINEXYRALTdisErr

06190670216.875200315.25103.36-4.23

16190670260.945200264.8796.26-5.72

26190670303.475200212.9990.17-6.83

36190670345.235200160.2990.89-5.89

46190670386.835200107.4186.48-6.73

56190670428.745200054.6583.71-8.89

for i in range(len(heightClass) - 1):

#根據(jù)飛行高度分類區(qū)間，獲取數(shù)據(jù)

h = df.loc[(df ['RALT'] > heightClass [i]) & (df['RALT'] <= heightClass [i + 1]), ['LINE', 'RALT']]

#獲取各測(cè)線在某個(gè)高度范圍的測(cè)點(diǎn)數(shù)量

indivualCounts =h.groupby('LINE')['RALT'].count()

lineIndex = df['LINE'].unique ( ) #依照測(cè)線飛行順序，提取此架次全部測(cè)線號(hào)

[6190, 6180, 6210, 6200]

for i, ee in enumerate(lineIndex):#迭代每條測(cè)線

#逐一獲取各架次數(shù)據(jù)

df_singleLn = df.loc[(df ['LINE'] == ee), ['LINE', 'X', 'Y', 'RALT', 'disErr']]

#保存數(shù)據(jù)(CSV格式)

df_singleLn.to_csv(myPath + str(ee) + '.csv', index=False)

….

4.3 磁日變數(shù)據(jù)的插值

磁法測(cè)量中，首先要進(jìn)行日變修正。但實(shí)際測(cè)量中，磁日變站數(shù)據(jù)的采集頻率往往小于測(cè)區(qū)數(shù)據(jù)采集頻率。

4.3.1 G858SX銫光泵磁力儀采集的磁日變數(shù)據(jù)

0 54834.178 0.000 16:26:04.60 09/01/14 0

0 54834.204 0.000 16:26:03.60 09/01/14 0

0 54834.166 0.000 16:26:02.60 09/01/14 0

0 54834.245 0.000 16:26:01.60 09/01/14 0

0 54834.268 0.000 16:26:00.60 09/01/14 0

0 54834.245 0.000 16:25:59.60 09/01/14 0

0 54834.292 0.000 16:25:58.60 09/01/14 0

…….

4.3.2 Python程序代碼示例

此段程序分為三部分，①檢查日變數(shù)據(jù)是否丟失；②磁日變數(shù)據(jù)插值；③磁日變數(shù)據(jù)的可視化。

1)數(shù)據(jù)檢查，主要利用Pandas完成

import pandas as pd

import numpy as np

from scipy.interpolate import interp1d

import matplotlib.pyplot as plt

import matplotlib as mpl

#讀取磁日變數(shù)據(jù)(其中3列)，含時(shí)間格式

df_rb = pd.read_table(dirualFile, date_parser=True, header=None, delimiter='s+',

na_values=['NULL'], usecols=(1, 3, 4))

df_rb.columns = ['mag', 'myTime', 'myDate']# 修改數(shù)據(jù)列標(biāo)題

#獲取測(cè)量起、止時(shí)間

myTime = df_rb['myTime']

timeFirst = myTime.values[-1]

timeEnd = myTime.values[0]

#獲取測(cè)量起、止日期

myDay = df_rb['myDate']

firstDay = myDay.values[-1]

endDay = myDay.values[0]

#獲取測(cè)量時(shí)間段，固定頻率(1秒)日期-時(shí)間的索引值

timeRange = pd.date_range(endDay + ' ' + timeEnd, firstDay + ' ' + timeFirst, freq='-1s')

(['2014-09-0116:26:04.600000', '2014-09-0116:26:03.600000',

'2014-09-0116:26:02.600000', '2014-09-0116:26:01.600000',

'2014-09-0116:26:00.600000', '2014-09-0116:25:59.600000',

'2014-09-0116:25:58.600000', … )]

#比較計(jì)算時(shí)間和實(shí)際數(shù)據(jù)采集個(gè)數(shù)，判定是否有數(shù)據(jù)丟失

if len(timeRange) == len(myTime):

myFlag = True

……

2)數(shù)據(jù)插值，主要利用Scipy和Numpy完成。

mag1 = df_rb['mag'] #獲取磁場(chǎng)值

mag1 = np.array(mag1)#轉(zhuǎn)換為Numpy數(shù)組

mylen, = mag1.shape

x1 = np.linspace(0, mylen - 1, mylen)#原始步長(zhǎng)

x1new = np.linspace(0, mylen - 1, (mylen - 1) * 10 + 1) #插值步長(zhǎng)

#選擇cubic方法，一維插值

func = interp1d(x1, mag1, kind='cubic')

y1new = func(x1new)#插值結(jié)果

…..

3)數(shù)據(jù)可視化，利用Matplotlib完成[11]。

fig = plt.figure(figsize=(10, 3))

ax1 = fig.add_subplot(1, 1, 1)

ax1.set_title('Charts of Magnetic Diurnal Variation')

#繪制含有日期格式的圖件(圖2)

ax1.plot_date(timeRange, y1new, linestyle='-', color='b')

#定義X、Y軸屬性

date_format = mpl.dates.DateFormatter('%H:%M')

ax1.xaxis.set_major_formatter(date_format)

ax1.set_xlabel('Time')

ax1.set_ylabel('nT')

fig.autofmt_xdate()

plt.grid(True) #繪制網(wǎng)格線

plt.tight_layout()#緊湊圖框

plt.show()#顯示結(jié)果

從上面四段代碼中可看出，Python程序利用豐富的第三方庫(kù)，代碼簡(jiǎn)潔、易讀，為快速數(shù)據(jù)清洗、分類統(tǒng)計(jì)、數(shù)據(jù)可視化提供強(qiáng)大的技術(shù)支持。

5 結(jié)語

不同的物探方法，所需求的數(shù)據(jù)格式不同。在某測(cè)點(diǎn)上，磁、放測(cè)量獲得一個(gè)或多個(gè)變量，可形成似表格式的數(shù)據(jù)集，而電磁法往往以一組或多組數(shù)據(jù)塊給出；在物探綜合測(cè)量中，數(shù)據(jù)由多種儀器分別采集，需進(jìn)行必要的數(shù)據(jù)合并。這就造成數(shù)據(jù)格式的多樣化。在數(shù)據(jù)處理和資料解釋過程中，相對(duì)于數(shù)據(jù)修正、數(shù)據(jù)網(wǎng)格化、正反演、可視化成圖等工作，數(shù)據(jù)清洗似乎顯得微不足道，但多樣化的物探數(shù)據(jù)格式，導(dǎo)致大部分時(shí)間和精力都花費(fèi)在此環(huán)節(jié)中。為了使數(shù)據(jù)處理人員和物探工作者從繁雜、瑣碎的工作中解脫出來，可借助更為高效的Python編程語言，規(guī)范相應(yīng)的數(shù)據(jù)結(jié)構(gòu)，避免重復(fù)工作，提高工作效率，將更多時(shí)間和精力致力于數(shù)據(jù)分析、地質(zhì)推斷解釋之中，從而提升項(xiàng)目的總體進(jìn)度和成果。

圖2 磁日變曲線(Matplotlib生成)Fig.2 Diurnal variation profile by Matplotlib