王振振

(葛洲壩集團試驗檢測有限公司，武漢 430015)

?

Python科學(xué)計算包在實驗數(shù)據(jù)處理中的應(yīng)用

王振振

(葛洲壩集團試驗檢測有限公司，武漢 430015)

Python科學(xué)計算包以其開源免費、面向多用途等優(yōu)點，已成功應(yīng)用于眾多科學(xué)計算領(lǐng)域。本文介紹了通過工程應(yīng)用中一個帶粗差剔除功能的線性回歸程序例子，簡要介紹了Python科學(xué)計算包在計量測試工作中實驗結(jié)果處理的應(yīng)用，展示了其強大的功能和易用性，為其在該領(lǐng)域的普及推廣提供了應(yīng)用經(jīng)驗。

Python；粗差剔除；數(shù)據(jù)處理；拉依達準則

0 引言

Python 是一種開源、跨平臺、通用型高級動態(tài)語言。作為一種廣泛使用的通用型程序設(shè)計語言，Python語言應(yīng)用領(lǐng)域涵蓋GUI編程、網(wǎng)絡(luò)通訊、科學(xué)計算、硬件通訊和多媒體編程等各個層面，同時還可以通過C語言進行擴展。Python弱類型、解釋型語言的特性，使得開發(fā)者可以交互式運行命令，方便及時檢驗、調(diào)試數(shù)據(jù)，非常適合科學(xué)計算編程。近些年來，隨著NumPy、SciPy、matplotlib等眾多函數(shù)庫的完善，Python已經(jīng)具備了足夠的功能，能夠充分滿足常見的科學(xué)計算需求[1]。Python語言形成的科學(xué)計算包，與流行的商業(yè)軟件Matlab相比，具備開源免費、面向多用途等優(yōu)點。本文將簡要介紹Python科學(xué)計算包，并通過一個帶粗差剔除功能的線性回歸程序例子，闡述Python科學(xué)計算包在實驗數(shù)據(jù)處理中的應(yīng)用。

1 Python科學(xué)計算包簡介

Python 科學(xué)計算包泛指Python語言下一系列面向科學(xué)計算的程序庫集合。它沒有確切的定義標準，使用者根據(jù)自己的需要選擇工具組合。當(dāng)前也存在高度整合、便于安裝的Python科學(xué)計算發(fā)行版如Python(x,y)、Anaconda等。這些發(fā)行版提供了基于圖形界面的工程管理、數(shù)據(jù)存取、數(shù)據(jù)可視化、編碼調(diào)試等一體化功能，安裝、使用非常方便，很大程度上推動了Python在數(shù)學(xué)、科學(xué)和工程計算領(lǐng)域的應(yīng)用普及。

計算機數(shù)值運算離不開向量計算，正是由于多數(shù)運算都是基于向量和矩陣而非單值，Matlab才具備了實現(xiàn)編制快速計算程序的能力。在Python語言中，第三方庫NumPy提供了向量計算所必需的數(shù)據(jù)結(jié)構(gòu)與基礎(chǔ)函數(shù)，從而成為了科學(xué)計算包的基礎(chǔ)。

數(shù)據(jù)結(jié)構(gòu)ndarray(n-dimensional array，多維數(shù)組，簡寫為array)及其相關(guān)運算函數(shù)即是NumPy的核心功能。ndarray是一種多維數(shù)組類型，通過下標訪問內(nèi)部元素，并支持切片(slice)，用以實現(xiàn)向量和矩陣。有別于Python 內(nèi)置類型list，它要求內(nèi)部元素具備相同的數(shù)據(jù)類型。NumPy中內(nèi)置了很多向量和矩陣的運算輔助函數(shù)，這些函數(shù)大多數(shù)都由C語言優(yōu)化實現(xiàn)快速計算，使用者無需擔(dān)心效率問題。以下代碼片段演示了ndarray向量基本運算。

>>> x = np.array([1, 2, 3, 4])

>>> x

array([1, 2, 3, 4])

>>> y=x**2 # 向量元素的二次方

>>> y

array([1, 4, 9, 16])

>>> y+1 # 向量元素與常數(shù)的加法

array([2, 5, 10, 17])

>>> x+y # 向量與向量的相加

array([2, 6, 12, 20])

>>> x.max(),x.mean(),x.std() # 向量的最大值、平均值、標準差

(4, 2.5, 1.1180339887498949)

NumPy提供了基本向量元素的三角函數(shù)、反三角函數(shù)、統(tǒng)計函數(shù)、隨機數(shù)等一系列常用功能。得益于Python弱類型動態(tài)語言的特性，正如上例所示，NumPy進行向量元素與常數(shù)運算、向量與向量運算沒有特別的語法機制進行區(qū)別，使用起來非常簡單方便。作為NumPy的擴充，SciPy函數(shù)庫提供了更為廣泛的科學(xué)計算函數(shù)，包括圖像處理、信號分析、線性方程組、插值與擬合、數(shù)理統(tǒng)計等一系列功能，限于篇幅本文不再作深入探討[2]。

此外，科學(xué)計算通常需要多樣的可視化工具來展示計算成果。Python科學(xué)計算常用的繪圖庫Matplotlib提供了與Matlab接口類似的交互式圖表。Matplotlib與NumPy數(shù)據(jù)結(jié)構(gòu)兼容，支持繪制曲線圖、直方圖、散點圖以及3D圖表，圖表中坐標、線形、標簽等元素都可以配置，生成的圖表還支持縮放交互，是一款非常強大的繪圖庫。下面的程序片段演示了Matplotlib的基本用法，運行結(jié)果見圖1。

>>> import matplotlib.pyplot as plt

>>> import numpy as np

>>> data = np.loadtxt(′D:data.csv′) # loadtxt載入csv數(shù)據(jù)文件

>>> plt.plot(data)

>>> plt.show()

圖1 Matplotlib繪制曲線圖形

2 帶有粗差剔除功能程序的實現(xiàn)

下面將結(jié)合某款正在研制中的混凝土性能自動測試系統(tǒng)的程序編制過程，介紹Python科學(xué)計算包在實驗數(shù)據(jù)處理中的應(yīng)用。

2.1 問題描述

該系統(tǒng)通過混凝土試件自動進行激振測試，將獲取的振動信號進行分析，計算特征值，研究特征值與混凝土凝結(jié)程度的關(guān)系。由于受材料不確定性的影響，振動測試結(jié)果一般都具有模糊性，然而經(jīng)實驗發(fā)現(xiàn)，單次特征值計算結(jié)果的模糊性，對混凝土性能指標隨時間變化的整體趨勢沒有影響。從圖2中可以發(fā)現(xiàn)，雖然振動測試特征值單次計算結(jié)果數(shù)值分布較為分散，但是整體上具備一定區(qū)間內(nèi)的線性關(guān)系。經(jīng)多次實驗驗證，通過人工作圖法在散點圖中剔除粗差，繪制穿過最多散點的直線段，就能得到測量值與混凝土凝結(jié)程度的對應(yīng)關(guān)系?，F(xiàn)在需要通過編制程序?qū)崿F(xiàn)這一過程，即自動對帶有粗差的實驗結(jié)果進行處理，尋找線性區(qū)間并對線性曲線直線段進行擬合。

圖2 實驗所得特征值與混凝土凝結(jié)時長的關(guān)系曲線

2.2 解決思路

為了從圖2中離散樣本的實驗結(jié)果求解出線性關(guān)系，需要解決曲線擬合、粗差剔除、尋找最佳擬合區(qū)間三個問題。

針對這種簡單的直線段擬合，本例采用最小二乘法實現(xiàn)?？紤]到用于統(tǒng)計計算的實驗結(jié)果數(shù)據(jù)較多(一般樣本數(shù)量n>50)，可以采取拉依達準則(即三倍標準差準則)對樣本數(shù)據(jù)進行粗差判別[3]?！坝美肋_準則判斷粗大誤差的基本思想是以給定的置信概率99.7%為標準，以三倍測量列的標準偏差限為依據(jù)，凡超過此界限的誤差，就認為它不屬于隨機誤差的范疇，而是粗大誤差。含有粗大誤差的測量值稱為異常值,異常值是不可取的，應(yīng)該從測量數(shù)據(jù)中剔除”[4]。

粗差剔除處理流程為：

1)對樣本數(shù)據(jù)使用最小二乘法進行初步擬合；

2)計算擬合所得離差，采用三倍標準差準則判斷；

3)若標準差超過評判準則，說明存在粗差。剔除異常值后，對剩余的樣本數(shù)據(jù)再次擬合，進入步驟1)；

4)若標準差不超過評判準則，說明粗差剔除完畢，所得擬合曲線滿足要求。

為了尋找最佳擬合區(qū)間，即對應(yīng)作圖法中尋找最長連續(xù)直線段，可采用簡單的遍歷機制，從小樣本容量逐步擴大至所有樣本，逐一進行粗差剔除和擬合運算，以擬合所得相關(guān)系數(shù)為標準判斷，保留最佳擬合結(jié)果。

2.3 程序編制

下面將簡述如何運用Python科學(xué)計算包實現(xiàn)程序編制。在本例的真實應(yīng)用環(huán)境中，由于采集所得樣本數(shù)量足夠大，為了更快速的得到線性區(qū)間的擬合曲線，程序首先對樣本數(shù)據(jù)進行了簡單的極值篩選。數(shù)值過小與過大的樣本不列入粗差判別的范圍。

2.3.1 最小二乘法擬合函數(shù)

類似于Matlab，Numpy函數(shù)庫中已經(jīng)提供了最小二乘法多項式擬合函數(shù)polyfit。調(diào)用polyfit(x, y, deg)函數(shù)可以對x、y兩個數(shù)組進行指定階數(shù)deg 的多項式擬合，返回相應(yīng)階數(shù)的多項式各項系數(shù)。在本例中，由于需要同時獲取多項式擬合計算過程中標準差，從運行效率考慮，沒有使用已有的polyfit函數(shù)，而是單獨編寫了可返回相關(guān)系數(shù)r的最小二乘法擬合函數(shù)leastsq。計算公式如下：

以下代碼演示了實現(xiàn)過程。

defleastsq(x,y):

″″″

一階多項式y(tǒng)=kx+b擬合，返回k、b以及相關(guān)系數(shù)r

″″″

x,y=numpy.array(x),numpy.array(y) # 生成array對象

meanx,meany=x.mean(),y.mean() # 求x、y數(shù)組的平均值

sumxy= (x*y).sum() # 向量相乘然后求和

xsum,ysum= 0.0, 0.0

foriinrange(len(x)):

xsum+= (x[i] -meanx)*(y[i]-meany)

ysum+= (x[i] -meanx)**2

k=xsum/ysum

b=meany-k*meanx

r=sum((x-meanx)*(y-meany)) /(math.sqrt(sum((x-meanx)**2)) *math.sqrt(sum((y-meany)**2)))

returnk,b,r#返回擬合的兩個系數(shù),相關(guān)系數(shù)

從實現(xiàn)過程可以看出，使用Python根據(jù)對照計算公式編制程序直觀明了，非常適合科研與工程技術(shù)人員等進行數(shù)據(jù)計算與處理。

2.3.2 剔除粗差與保留最佳擬合區(qū)間

該系統(tǒng)的實際應(yīng)用要求程序即時對已采集到的數(shù)據(jù)進行分析，即振動測試與數(shù)據(jù)處理自動執(zhí)行。同時考慮三倍標準差準則對樣本數(shù)量的要求，程序中設(shè)定當(dāng)樣本數(shù)量n大于15時才進行數(shù)據(jù)粗差剔除與擬合處理[5]。從樣本中末端15筆數(shù)據(jù)開始，逐一擴大計算的樣本數(shù)量，依次對末端15、16、17……直至全部數(shù)據(jù)進行粗差剔除和擬合運算，保留該過程中最佳擬合結(jié)果作為最終的擬合結(jié)果。

最終實現(xiàn)代碼片段如下：

ANS= {′k′:None, ′b′:None, ′r′:0}

if15 >len(sx):

return

forninrange(15,len(DATA_X)):

done=False

sx=copy.copy(DATA_X[-n:]) # 取末端數(shù)據(jù)進行以下粗差剔除和擬合

sy=copy.copy(DATA_Y[-n:])

# 開始粗差剔除

whilenotdone:

x=numpy.array(sx)

y=numpy.array(sy)

if5 >=len(sx): # 剔除后數(shù)據(jù)不足則退出

break

k,b,r=leastsq(x,y)

xi=x

yi=xi*k+b

e=abs(y-yi)

s=math.sqrt(sum(e**2)/(len(x)-2))

ifmax(e) >= 3*s: # 不滿足3倍標準差準則

i=e.argmax() #argmax用于查找極值對應(yīng)下標

sx.pop(i) # 從數(shù)組中彈出離差最大的樣本數(shù)據(jù)

sy.pop(i)

else:

done=True

ifdoneandr>ANS[′r′]: # 擬合結(jié)果優(yōu)于已有結(jié)果

ANS.update({′k′:k, ′b′:b,′r′:r}

在粗差剔除過程中，程序中設(shè)定變量done標志擬合是否完成，不斷剔除異常值，當(dāng)異常值剔除后不具備擬合條件時，done保持未完成狀態(tài)，求解過程退出；當(dāng)離差不大于3倍標準差時，將done置于完成狀態(tài)，求解完成。如果求解得到的相關(guān)系數(shù)r優(yōu)于已有求解，將當(dāng)前求解結(jié)果作為最佳結(jié)果保留。

3 結(jié)語

本文簡要介紹了Python科學(xué)計算包的特點及其主要組成，闡述了一個帶粗差剔除功能的實驗結(jié)果處理例子，展示了Python科學(xué)計算包提供的強大功能和易用性。利用免費開源的Python科學(xué)計算開發(fā)環(huán)境，設(shè)計了數(shù)據(jù)處理算法并用程序代碼實現(xiàn)了對實際工作中的測試數(shù)據(jù)的自動處理，避免了使用昂貴的商業(yè)軟件，解放了人力，具有較高的經(jīng)濟性和實用價值。本文的介紹，能夠為計量、測試和檢驗等相關(guān)技術(shù)人員提供應(yīng)用經(jīng)驗，有利于Python科學(xué)計算包在計量測試學(xué)科和行業(yè)中的普及推廣。

[1] 張若愚.Python科學(xué)計算[M].清華大學(xué)出版社, 2012.10-12

[2]JonesE,OliphantE,PetersonP,etal.SciPy:OpenSourceScientificToolsforPython, 2001-,http://www.scipy.org/ [Online;accessed2015-05-18]

[3] 孫培強.正確選擇統(tǒng)計判別法剔除異常值[J].計量技術(shù), 2013(11)

[4] 張敏,袁輝.拉依達(PauTa)準則與異常值剔除[J].鄭州工業(yè)大學(xué)學(xué)報, 1997(1)

[5] 蔣珍美,吳先球,陳俊芳.一個具有粗差剔除功能的線性回歸程序[J].華南師范大學(xué)學(xué)報(自然科學(xué)版), 2002(2)

《計量技術(shù)》雜志歡迎大家踴躍投稿《計量技術(shù)》雜志以實用性、權(quán)威性、及時性為主要特色；堅持面向生產(chǎn)，面向基層，理論與實踐相結(jié)合的編輯方針；著重報道計量、測試、檢驗、質(zhì)量保障等方面的新技術(shù)、新產(chǎn)品、新動態(tài)、綜合評述、經(jīng)驗介紹等內(nèi)容。歡迎大家投稿，投稿要求如下：1 來稿以說明問題為主要目的，語句要精練、簡潔，全文盡量不超過4000字。2000字以上請附摘要和關(guān)鍵詞。2 來稿涉及計量單位時，請一律使用法定計量單位的名稱和符號。3 來稿時應(yīng)寫清作者姓名、單位、郵編、通訊地址及聯(lián)系電話。可通過我刊電子信箱投稿。4 編輯部收到來稿后，立即給作者“稿件回執(zhí)”，在不遲于四個月內(nèi)通知作者是否刊用。由于本刊人力有限，不刊用稿均不寄還，請作者自留底稿。

10.3969/j.issn.1000-0771.2015.07.24