Python科学计算包在实验数据处理中的应用

王振振

(葛洲坝集团试验检测有限公司，武汉 430015)



Python科学计算包在实验数据处理中的应用

王振振

(葛洲坝集团试验检测有限公司，武汉 430015)

Python科学计算包以其开源免费、面向多用途等优点，已成功应用于众多科学计算领域。本文介绍了通过工程应用中一个带粗差剔除功能的线性回归程序例子，简要介绍了Python科学计算包在计量测试工作中实验结果处理的应用，展示了其强大的功能和易用性，为其在该领域的普及推广提供了应用经验。

Python；粗差剔除；数据处理；拉依达准则

0 引言

Python 是一种开源、跨平台、通用型高级动态语言。作为一种广泛使用的通用型程序设计语言，Python语言应用领域涵盖GUI编程、网络通讯、科学计算、硬件通讯和多媒体编程等各个层面，同时还可以通过C语言进行扩展。Python弱类型、解释型语言的特性，使得开发者可以交互式运行命令，方便及时检验、调试数据，非常适合科学计算编程。近些年来，随着NumPy、SciPy、matplotlib等众多函数库的完善，Python已经具备了足够的功能，能够充分满足常见的科学计算需求[1]。Python语言形成的科学计算包，与流行的商业软件Matlab相比，具备开源免费、面向多用途等优点。本文将简要介绍Python科学计算包，并通过一个带粗差剔除功能的线性回归程序例子，阐述Python科学计算包在实验数据处理中的应用。

1 Python科学计算包简介

Python 科学计算包泛指Python语言下一系列面向科学计算的程序库集合。它没有确切的定义标准，使用者根据自己的需要选择工具组合。当前也存在高度整合、便于安装的Python科学计算发行版如Python(x,y)、Anaconda等。这些发行版提供了基于图形界面的工程管理、数据存取、数据可视化、编码调试等一体化功能，安装、使用非常方便，很大程度上推动了Python在数学、科学和工程计算领域的应用普及。

计算机数值运算离不开向量计算，正是由于多数运算都是基于向量和矩阵而非单值，Matlab才具备了实现编制快速计算程序的能力。在Python语言中，第三方库NumPy提供了向量计算所必需的数据结构与基础函数，从而成为了科学计算包的基础。

数据结构ndarray(n-dimensional array，多维数组，简写为array)及其相关运算函数即是NumPy的核心功能。ndarray是一种多维数组类型，通过下标访问内部元素，并支持切片(slice)，用以实现向量和矩阵。有别于Python 内置类型list，它要求内部元素具备相同的数据类型。NumPy中内置了很多向量和矩阵的运算辅助函数，这些函数大多数都由C语言优化实现快速计算，使用者无需担心效率问题。以下代码片段演示了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 # 向量元素与常数的加法

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

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

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

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

(4, 2.5, 1.1180339887498949)

NumPy提供了基本向量元素的三角函数、反三角函数、统计函数、随机数等一系列常用功能。得益于Python弱类型动态语言的特性，正如上例所示，NumPy进行向量元素与常数运算、向量与向量运算没有特别的语法机制进行区别，使用起来非常简单方便。作为NumPy的扩充，SciPy函数库提供了更为广泛的科学计算函数，包括图像处理、信号分析、线性方程组、插值与拟合、数理统计等一系列功能，限于篇幅本文不再作深入探讨[2]。

此外，科学计算通常需要多样的可视化工具来展示计算成果。Python科学计算常用的绘图库Matplotlib提供了与Matlab接口类似的交互式图表。Matplotlib与NumPy数据结构兼容，支持绘制曲线图、直方图、散点图以及3D图表，图表中坐标、线形、标签等元素都可以配置，生成的图表还支持缩放交互，是一款非常强大的绘图库。下面的程序片段演示了Matplotlib的基本用法，运行结果见图1。

>>> import matplotlib.pyplot as plt

>>> import numpy as np

>>> data = np.loadtxt(′D:data.csv′) # loadtxt载入csv数据文件

>>> plt.plot(data)

>>> plt.show()

图1 Matplotlib绘制曲线图形

2 带有粗差剔除功能程序的实现

下面将结合某款正在研制中的混凝土性能自动测试系统的程序编制过程，介绍Python科学计算包在实验数据处理中的应用。

2.1 问题描述

该系统通过混凝土试件自动进行激振测试，将获取的振动信号进行分析，计算特征值，研究特征值与混凝土凝结程度的关系。由于受材料不确定性的影响，振动测试结果一般都具有模糊性，然而经实验发现，单次特征值计算结果的模糊性，对混凝土性能指标随时间变化的整体趋势没有影响。从图2中可以发现，虽然振动测试特征值单次计算结果数值分布较为分散，但是整体上具备一定区间内的线性关系。经多次实验验证，通过人工作图法在散点图中剔除粗差，绘制穿过最多散点的直线段，就能得到测量值与混凝土凝结程度的对应关系。现在需要通过编制程序实现这一过程，即自动对带有粗差的实验结果进行处理，寻找线性区间并对线性曲线直线段进行拟合。

图2 实验所得特征值与混凝土凝结时长的关系曲线

2.2 解决思路

为了从图2中离散样本的实验结果求解出线性关系，需要解决曲线拟合、粗差剔除、寻找最佳拟合区间三个问题。

针对这种简单的直线段拟合，本例采用最小二乘法实现。考虑到用于统计计算的实验结果数据较多(一般样本数量n>50)，可以采取拉依达准则(即三倍标准差准则)对样本数据进行粗差判别[3]。“用拉依达准则判断粗大误差的基本思想是以给定的置信概率99.7%为标准，以三倍测量列的标准偏差限为依据，凡超过此界限的误差，就认为它不属于随机误差的范畴，而是粗大误差。含有粗大误差的测量值称为异常值,异常值是不可取的，应该从测量数据中剔除”[4]。

粗差剔除处理流程为：

1)对样本数据使用最小二乘法进行初步拟合；

2)计算拟合所得离差，采用三倍标准差准则判断；

3)若标准差超过评判准则，说明存在粗差。剔除异常值后，对剩余的样本数据再次拟合，进入步骤1)；

4)若标准差不超过评判准则，说明粗差剔除完毕，所得拟合曲线满足要求。

为了寻找最佳拟合区间，即对应作图法中寻找最长连续直线段，可采用简单的遍历机制，从小样本容量逐步扩大至所有样本，逐一进行粗差剔除和拟合运算，以拟合所得相关系数为标准判断，保留最佳拟合结果。

2.3 程序编制

下面将简述如何运用Python科学计算包实现程序编制。在本例的真实应用环境中，由于采集所得样本数量足够大，为了更快速的得到线性区间的拟合曲线，程序首先对样本数据进行了简单的极值筛选。数值过小与过大的样本不列入粗差判别的范围。

2.3.1 最小二乘法拟合函数

类似于Matlab，Numpy函数库中已经提供了最小二乘法多项式拟合函数polyfit。调用polyfit(x, y, deg)函数可以对x、y两个数组进行指定阶数deg 的多项式拟合，返回相应阶数的多项式各项系数。在本例中，由于需要同时获取多项式拟合计算过程中标准差，从运行效率考虑，没有使用已有的polyfit函数，而是单独编写了可返回相关系数r的最小二乘法拟合函数leastsq。计算公式如下：

以下代码演示了实现过程。

defleastsq(x,y):

″″″

一阶多项式y=kx+b拟合，返回k、b以及相关系数r

″″″

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

meanx,meany=x.mean(),y.mean() # 求x、y数组的平均值

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#返回拟合的两个系数,相关系数

从实现过程可以看出，使用Python根据对照计算公式编制程序直观明了，非常适合科研与工程技术人员等进行数据计算与处理。

2.3.2 剔除粗差与保留最佳拟合区间

该系统的实际应用要求程序即时对已采集到的数据进行分析，即振动测试与数据处理自动执行。同时考虑三倍标准差准则对样本数量的要求，程序中设定当样本数量n大于15时才进行数据粗差剔除与拟合处理[5]。从样本中末端15笔数据开始，逐一扩大计算的样本数量，依次对末端15、16、17……直至全部数据进行粗差剔除和拟合运算，保留该过程中最佳拟合结果作为最终的拟合结果。

最终实现代码片段如下：

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

if15 >len(sx):

return

forninrange(15,len(DATA_X)):

done=False

sx=copy.copy(DATA_X[-n:]) # 取末端数据进行以下粗差剔除和拟合

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

# 开始粗差剔除

whilenotdone:

x=numpy.array(sx)

y=numpy.array(sy)

if5 >=len(sx): # 剔除后数据不足则退出

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用于查找极值对应下标

sx.pop(i) # 从数组中弹出离差最大的样本数据

sy.pop(i)

else:

done=True

ifdoneandr>ANS[′r′]: # 拟合结果优于已有结果

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

在粗差剔除过程中，程序中设定变量done标志拟合是否完成，不断剔除异常值，当异常值剔除后不具备拟合条件时，done保持未完成状态，求解过程退出；当离差不大于3倍标准差时，将done置于完成状态，求解完成。如果求解得到的相关系数r优于已有求解，将当前求解结果作为最佳结果保留。

3 结语

本文简要介绍了Python科学计算包的特点及其主要组成，阐述了一个带粗差剔除功能的实验结果处理例子，展示了Python科学计算包提供的强大功能和易用性。利用免费开源的Python科学计算开发环境，设计了数据处理算法并用程序代码实现了对实际工作中的测试数据的自动处理，避免了使用昂贵的商业软件，解放了人力，具有较高的经济性和实用价值。本文的介绍，能够为计量、测试和检验等相关技术人员提供应用经验，有利于Python科学计算包在计量测试学科和行业中的普及推广。

[1] 张若愚.Python科学计算[M].清华大学出版社, 2012.10-12

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

[3] 孙培强.正确选择统计判别法剔除异常值[J].计量技术, 2013(11)

[4] 张敏,袁辉.拉依达(PauTa)准则与异常值剔除[J].郑州工业大学学报, 1997(1)

[5] 蒋珍美,吴先球,陈俊芳.一个具有粗差剔除功能的线性回归程序[J].华南师范大学学报(自然科学版), 2002(2)

《计量技术》杂志欢迎大家踊跃投稿《计量技术》杂志以实用性、权威性、及时性为主要特色；坚持面向生产，面向基层，理论与实践相结合的编辑方针；着重报道计量、测试、检验、质量保障等方面的新技术、新产品、新动态、综合评述、经验介绍等内容。欢迎大家投稿，投稿要求如下：1 来稿以说明问题为主要目的，语句要精练、简洁，全文尽量不超过4000字。2000字以上请附摘要和关键词。2 来稿涉及计量单位时，请一律使用法定计量单位的名称和符号。3 来稿时应写清作者姓名、单位、邮编、通讯地址及联系电话。可通过我刊电子信箱投稿。4 编辑部收到来稿后，立即给作者“稿件回执”，在不迟于四个月内通知作者是否刊用。由于本刊人力有限，不刊用稿均不寄还，请作者自留底稿。

10.3969/j.issn.1000-0771.2015.07.24