MATLAB对静电场等势线的最小二乘描绘

刘雅彬，鲁晓东

（浙江海洋学院，浙江 舟山 316000）

在静电场描绘实验中，通常采用的方法是用一对电极产生的恒定电流模拟一对等量异种电荷产生的静电场。直接测电流场有一定的困难，而借助电压表可直接测出电场中的一系列等势点进而画出等势线，再根据等势线和电场线正交关系便可直接画出电场线［1－2］。由于在电势的探测过程中，介质分布的不稳定性以及操作时探针的抖动都会造成记录点的不准，使手绘的等势线具有较大的不确定度性，增大了电场线描绘畸变的可能性［3］。提高仪器的精度是解决问题的一种方法，但效果往往不明显，原因是流体介质的稳定性不好控制，操作时呈现出较大的随机性，因此可以运用最小二乘的方法去拟合这些点的曲线，使测量误差得到最大的弱化。而使用MATLAB则可以直接利用它提供的优化工具箱，简化编程过程，并使处理的结果能以图形方式直观地表达［4］。

1 最小二乘曲线拟合原理

实验中经常根据许多组观测数据来确定它们的函数曲线，这就是实验数据处理中的曲线拟合问题。曲线拟合的目的是就是要把数据点的固有规律表现出来，但误差的存在会使拟合产生一定程度的失真，所以必须考虑通过怎样的途径使误差影响程度最小［5］。

最小二乘准则就是使实验值与拟合值的残差平方和达到最小的条件下进行函数参数的估计。不妨设数据（xi，yi），估计参数为（α，β），拟合的函数为y＝Φ（x，α，β），最小二乘法就是使误差yi－Φ（xi，α，β）的平方和最小，即

实际上是使拟合的曲线与各数据点的距离最均匀，以此达到弱化误差对结果的影响，达到优化曲线的目的［6］。

2 Matlab等势线的描绘

2.1 等势点的输入

当实验的等势点数据记录到坐标纸后，必须把数据点的坐标值（x，y）存储到矩阵中以便于MATLAB的读取及处理。为区分不同的等势线或不同的等势线采集的点数量可能不一样，每条等势线的采集点坐标必须存储到一个独立的二维矩阵。

2.2 等势线的拟合

MATLAB软件提供了很多基于最小二乘准则的曲线拟合函数，若曲线形状是确定的，则可以用Lsqcur vefit（）。它是一个非线性最小二乘拟合函数，本质上是求解目标函数的优化问题。其使用典型格式为 ［z，resnor m］＝ lsqcur vefit（f un，beta0，xdata，ydata）其中f un是要拟合的非线性函数，简单的形式可以用inline函数进行定义，如f（x，y）＝a＊sin（b＊x＋y），f＝inline（＇a＊sin（b＊x＋y）＇，＇a＇，＇b，＇x＇，＇y＇）。如复杂可以写成m文件，具体参考相关书籍［7］；beta0是估计参数初始值，xdata，ydata是拟合点的数据即实验测得数据，z是返回的系数矩阵，即拟合的曲线特性参数，resnor m为残差平方和［8］。

对于曲线未定得的可用多项式最小二乘拟合，其命令是A＝polyfit（x，y，m）其中x表示函数中的自变量矩阵，y表示因变量矩阵，A是输出的系数矩阵，即多项式的系数。多项式在自变量x处的拟合值y可用以下命令计算：y＝polyval（A，x），进而算得拟合曲线的残差平方和。

3 应用例子及效果分析

图1 手工描绘等势点及电场线

以同轴圆柱体水槽模型形成的模拟场为例，对其进行处理。首先看由传统方法描绘的静电场，各环分布不均匀图1并且比较毛躁。由Matlab处理时，则先根据模型得到拟合方程形式并写为：这样把坐标输入点（x，y）看作函数的自变量，而f（x，y）＝0就是拟合值f（x，y，a，b，r）的目标值，优化的目标函数就简化为 ∑［f（x，y，a，b，r）－0］2。当把这些坐标纸上的数据点转存到matlab后，可直接调用函数Lsqcurvefit（），就得到拟合方程的估计参数（具体程序见附录），方程的图形化结果如图2。为比较他们的准确性，可以计算它们的百分差分布，其中理论值［9］的计算按（1）

VA，a，b为实验所测量的值，r为点到圆心的距离，百分差计算结果如表1，可以看出最小二乘法比手工方法估计的参数更接近实验值。

图2 matlab等势点描绘及电场线

表1 实验结果比较

4 结 论

从示例中看出利用Matlab拟合静电场的曲线非常方便，也适合其他实验曲线的拟合。使用时要注意的是Matlab提供的最小二乘工具都是在一定范围内寻找拟合函数的局部最优解，所以一定要结合实验条件对搜索区间、初始值等加以约束，提高最佳参数的命中率和程序运行速度。

程序附录：（拟合一个圆方程，其他等势点与此类同）

L1＝［－0.95 0；－0.6 0.6；0 0.85；0.7 0.6；0.9 0；0.6－0.65；0－0.9；－0.6－0.65］；%数据由坐标纸上读得

plot（L1（：，1），L1（：，2），＇o＇）；grid on；%打印记录点

axis（［－5，5，－5，5］）；hol d on；%坐标纸大小

beta0＝［0，0，1］；%在圆心为（0，0），半径为1附近搜索参数

z＝zeros（size（L1，1），1）；

%写出拟合方程形式

myf uncircle＝inline（＇（L1（：，1）－beta（1））.＾2＋（L1（：，2）－beta（2））.＾2－beta（3）＾2＇，＇beta＇，＇L1＇）；［beta，resnor m］＝lsqcurvefit（myf uncircle，beta0，L1，z）；%最小二乘拟合

a＝beta（1）；b＝beta（2）；r＝beta（3）；%返回拟合参数

h＝ezplot（＠（x，y）（x－a）＾2＋（y－b）＾2－r＾2，［－5 5－5 5］）；%打印拟合圆

hold on；set（h，＇Color＇，＇r＇，＇linestyle＇，＇－－＇，＇linewidt h＇，2）；

［1］张雅男，徐飞，叶影.Matlab模拟静电场与模拟静电实验的比较［J］.物理与工程，2008，18（2）：35－37.

［2］罗明海，韩亚萍，等.最小二乘法在伏安法测电阻实验中的应用［J］.大学物理实验，2012，25（2）：88－90.

［3］冯君，吕军.最小二乘法与作图法在物理实验中的应用研究［J］.黑龙江大学：自然科学版，1996，13（4）：83－87.

［4］葛一兵，王纪俊，卜敏.基于 MATLAB语言的静电场实验数据处理方法［J］.大学物理实验，2002，15（3）：70－71.

［5］陈敏.应用Matlab拟合传感器特征曲线［J］.南京师范大学学报：工程技术版，2003，3（1）：45－49.

［6］齐宝权.采用中心化最小二乘法进行测井曲线拟合［J］.测井技术，2007，31（4）：331－334.

［7］赵静，但琦.数学建模与数学实验［M］.3版.北京：高等教育出版社，2008.

［8］王广斌，刘义伦，金晓宏等.基于最小二乘原理的趋势处理以及 Matlab的实现［J］.有色设备，2005（4）：4－7.

［9］竺江峰，鲁晓东，夏雪琴.大学物理实验［M］.北京：中国水电水利出版社，2011.