王振振

(葛洲壩集團(tuán)試驗(yàn)檢測(cè)有限公司，武漢 430015)

?

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

王振振

(葛洲壩集團(tuán)試驗(yàn)檢測(cè)有限公司，武漢 430015)

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

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

0 引言

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

1 Python科學(xué)計(jì)算包簡(jiǎn)介

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

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

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

>>> 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() # 向量的最大值、平均值、標(biāo)準(zhǔn)差

(4, 2.5, 1.1180339887498949)

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

此外，科學(xué)計(jì)算通常需要多樣的可視化工具來(lái)展示計(jì)算成果。Python科學(xué)計(jì)算常用的繪圖庫(kù)Matplotlib提供了與Matlab接口類似的交互式圖表。Matplotlib與NumPy數(shù)據(jù)結(jié)構(gòu)兼容，支持繪制曲線圖、直方圖、散點(diǎn)圖以及3D圖表，圖表中坐標(biāo)、線形、標(biāo)簽等元素都可以配置，生成的圖表還支持縮放交互，是一款非常強(qiáng)大的繪圖庫(kù)。下面的程序片段演示了Matplotlib的基本用法，運(yùn)行結(jié)果見(jiàn)圖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 帶有粗差剔除功能程序的實(shí)現(xiàn)

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

2.1 問(wèn)題描述

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

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

2.2 解決思路

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

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

粗差剔除處理流程為：

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

2)計(jì)算擬合所得離差，采用三倍標(biāo)準(zhǔn)差準(zhǔn)則判斷；

3)若標(biāo)準(zhǔn)差超過(guò)評(píng)判準(zhǔn)則，說(shuō)明存在粗差。剔除異常值后，對(duì)剩余的樣本數(shù)據(jù)再次擬合，進(jìn)入步驟1)；

4)若標(biāo)準(zhǔn)差不超過(guò)評(píng)判準(zhǔn)則，說(shuō)明粗差剔除完畢，所得擬合曲線滿足要求。

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

2.3 程序編制

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

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

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

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

defleastsq(x,y):

″″″

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

″″″

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

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#返回?cái)M合的兩個(gè)系數(shù),相關(guān)系數(shù)

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

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

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

最終實(shí)現(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ù)進(jìn)行以下粗差剔除和擬合

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倍標(biāo)準(zhǔn)差準(zhǔn)則

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

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}

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

3 結(jié)語(yǔ)

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

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

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

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

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

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

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

10.3969/j.issn.1000-0771.2015.07.24