Python語言在物探數據清洗中的探索

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

(核工業航測遙感中心 鈾資源地球物理勘查技術中心重點實驗室，石家莊 050002)

0 前言

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

1 物探數據格式分析

1.1 潔凈數據的屬性

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

表1 測線飛行高度質量統計報表

1.2 常見不潔凈數據

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

1.2.1 數值為列標頭而非變量名稱

表1為某一架次航空物探飛行高度質量統計報表[3]。此表格從閱讀角度上，設計較為合理。但飛行高度(范圍)作為表頭(變量)給出，程序將使用更多代碼進行讀取。

此時，將表1中的飛行高度作為變量，表2的數據結構更適合于程序閱讀。

1.2.2 一列中存儲多個變量

物探數據中常出現諸如經緯度(XY坐標)、測量值的變化范圍等復合參數。表3中的電阻率范圍[4]存儲在同一列中。此時應將“變化范圍”分割為“最小值”和“最大值”兩列(表4)。

表2 測線飛行高度統計報表

表3 相山地區主要巖石電阻率統計一覽表

表4 相山地區主要巖石電阻率統計一覽表

1.2.3 變量存儲在行、列中

表5顯示了某地區航空測量設計測線。坐標值(端點)和線號(端點)變量分別處于行、列中。為了計算方便，應將同條測線的所有變量融合在同一行中(表6)。

表5 黑龍江省某地區航空物探測量設計測線一覽表(部分)

表6 黑龍江省某地區航空物探測量設計測線一覽表(部分)

1.2.4 一個數據集(表)中存儲多個類型

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

1.2.5 同類型數據處于多個數據集中

物探測量工作一般分為剖面測量(如電法)和面積性測量(如磁法)。測量數據以測點為基本單元，多個測點構成測線(剖面線)數據，再由多個測線(剖面線)形成一定網度的測量數據。所以數據清洗應根據測量性質來決定數據結構，單剖面抑或全區數據。

另外，由于測量手段限制，觀測值有時會存儲在兩個或多個數據集中，最為典型的是GPS定位數據與地面物探實測數據處在不同數據文件中。

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

1.2.6 數據集變量存有差異

物探測量中，同種儀器因采用不同的測量模式、測量裝置、采樣時間等，獲取的結果數據也不同(如大地電磁測量模式可分為標量和張量兩種)；激電測量(IPR-12多功能電法儀)會因采樣時間長短產生不同的數據格式。此時，視具體測量模式，采用枚舉法讀取數據。

1.2.7 丟碼或誤碼現象

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

從上述7種現象可看出，不潔凈數據①源于數據質量;②數據結構。但隨物探儀器的性能、處理程序的智能化不斷提升，數據質量隨之提高、數據結構也趨于規范化，所以物探數據中的不潔凈主要源于數據結構。

從某種意義上講，數據清理是一個數據結構轉換過程。在數據預處理中，應首先清洗數據本身的假值、空值或畸變值等，然后再進行數據結構的轉換。最終形成適合于計算機程序運算的數據結構，即潔凈數據。

對于EDI、AVG格式等大地電磁法中的通用數據結構，事實上在指定程序執行反演計算時，也需要一個格式轉換的過程；同樣在進行其他統計、分析，仍要依照潔凈數據結構提取、整理相關參量。

2 物探數據清洗方法

2.1 數據清洗中常用的Python第三方庫

Python是一種面向對象的解釋性的跨平臺高級計算機語言，具有語法簡單、可移植性強的特點，同時，具有極為豐富的類庫，廣泛應用于數據庫、多媒體、科學計算、網絡等諸多領域。正因Python具有強大的科學計算、繪圖功能，以及在地圖繪制、地理空間數據的處理與轉換方面具有豐富的類庫[5]，國外許多著名軟件(如arcGIS地理信息系統、Montaj Oasis地球物理處理系統等)以它作為主要腳本語言。

在數據清洗中，常用的Python第三方庫主要包括以下四個：

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

2)NumPy：運算速度非?？斓臄U展數學庫，可支持N維數組和矩陣計算并提供了大量的數學函數庫，具有強大的廣播功能[8-9]。與Pandas交互使用，可有效地提高運算速度。

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

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

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

總之，Python在數據清理或預處理中，優勢明顯，可大大縮短開發周期。

2.2 常規數據清洗操作

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

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

2)篩選：基于某個或多個條件，對數據集進行分組和移除變量，如測線分離，數據的分類統計等,主要由Pandas完成。

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

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

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

3 物探標準數據格式的建立

在潔凈數據和不潔凈數據屬性的基礎上，結合物探數據特點，對標準數據格式的建立提出如下幾點建議。

3.1 數據結構的構建

3.1.1 參據名稱

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

3.1.2 CSV格式

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

3.1.3 HDF5格式

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

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

3.2 數據集的管理

一般情況下，一個物探測區由多條測線(剖面)或測網構成。對于面積型測量，應將所有測線合并為一個文件，易于數據網格化；單測線分析或繪制剖面平面圖時，再進行測線分離。對于剖面測量，剖面數據應分別存儲，確保格式的統一，也可通過文件列表進行批量處理。

表7 某地區航放早晚基線測試256道譜數據(部分)

圖1 HDF5結構示例(256道能譜數據)Fig.1 1 Structure of HDF5 (256 channel spectral data)

4 Python在數據清洗中的應用

4.1 數據結構的清洗

此例是將不潔凈數據(表1)轉換為潔凈數據(表2)，展示Pandas庫強大的文件讀寫和數據結構合并等功能。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列數據

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

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

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

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

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

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

4.2 航測數據清洗和質量評價

對于航空物探，偏航距或離地高度是衡量航放測量飛行質量的兩個重要指標。下面以黑龍江省完達山—太平嶺地區航測項目中某一個架次的航測數據為例，展示使用Pandas庫進行測線分離后的分類統計。

1)由Pico公司生產的綜合航空物探測量系統采集的航測能譜數據(部分)：

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

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坐標 WGS84坐標系 GAUSS 22N投影 CM=129度

雷達高度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] #飛行高度分類級別

data_file = r' cq_B4090106_C_cut_Spc.dat'

#跳過15行，讀取數據中0, 4,5,7共4列數據，

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

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

#斜字體為輸出結果，下同。

LINEXYRALT

06190670216.875200315.25103.36

16190670260.945200264.8796.26

26190670303.475200212.9990.17

36190670345.235200160.2990.89

46190670386.835200107.4186.48

#通過設計測線，計算飛行偏離度

disErrSet = statPathDev (df, designLines)

#df數據集中增加飛行偏離度列(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):

#根據飛行高度分類區間，獲取數據

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

#獲取各測線在某個高度范圍的測點數量

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

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

[6190, 6180, 6210, 6200]

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

#逐一獲取各架次數據

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

#保存數據(CSV格式)

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

….

4.3 磁日變數據的插值

磁法測量中，首先要進行日變修正。但實際測量中，磁日變站數據的采集頻率往往小于測區數據采集頻率。

4.3.1 G858SX銫光泵磁力儀采集的磁日變數據

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程序代碼示例

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

1)數據檢查，主要利用Pandas完成

import pandas as pd

import numpy as np

from scipy.interpolate import interp1d

import matplotlib.pyplot as plt

import matplotlib as mpl

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

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']# 修改數據列標題

#獲取測量起、止時間

myTime = df_rb['myTime']

timeFirst = myTime.values[-1]

timeEnd = myTime.values[0]

#獲取測量起、止日期

myDay = df_rb['myDate']

firstDay = myDay.values[-1]

endDay = myDay.values[0]

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

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', … )]

#比較計算時間和實際數據采集個數，判定是否有數據丟失

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

myFlag = True

……

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

mag1 = df_rb['mag'] #獲取磁場值

mag1 = np.array(mag1)#轉換為Numpy數組

mylen, = mag1.shape

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

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

#選擇cubic方法，一維插值

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

y1new = func(x1new)#插值結果

…..

3)數據可視化，利用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) #繪制網格線

plt.tight_layout()#緊湊圖框

plt.show()#顯示結果

從上面四段代碼中可看出，Python程序利用豐富的第三方庫，代碼簡潔、易讀，為快速數據清洗、分類統計、數據可視化提供強大的技術支持。

5 結語

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

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