CN102054180A - 一种自动检测结霜的方法 - Google Patents

Abstract

本发明公开了一种自动检测结霜的方法，利用对设置在户外的多个载玻片的图像序列I_i的处理，检测出是否结霜及结霜的时刻，其中，所述图像序列I_i由每隔固定时间采集的载玻片图像所组成，所述载玻片呈四边形，i是正整数，表示图像帧数序号。本发明的方法能够自动检测结霜现象是否发生，并且检测出结霜发生的时刻，而且检测过程中除了初始手工标记外，不需要人工干预，自动化程度高，对于不同形式的结霜都能很好的检测。

Description

一种自动检测结霜的方法 

技术领域

本发明属图像处理和气象观测交叉技术领域，具体涉及一种基于图像处理的自动检测结霜的方法。 

背景技术

霜的观测是地面气象观测的一项重要内容。在《地面气象观测规范》(标准编号QX/T 46-2007)中，对结霜这一现象的定义是：霜是水汽在地面和近地面物体上凝华而成的白色松脆的冰晶；或由露冻结而成的冰珠。 

结霜需要在一定的风速、湿度、温度条件下才能形成。风速过大，不利于水蒸气的聚集凝华；温度过高，湿度过大将会先形成大量的露，而露的凝结过程必将伴随着凝结潜热的释放，将部分补偿了因辐射而损失的热量，不利于霜的形成。霜一般容易在晴朗小风且寒冷的夜间生成。 

由于霜形成的天气条件比较复杂，所以目前霜的观测和记录主要通过有经验的观测员人工完成。因此结霜的自动检测对提高检测准确率，减少观测员工作量具有重要意义。 

所谓拟合是指已知某函数的若干离散函数值，通过调整该函数中若干待定参数，使得该函数与已知点集的差别(最小二乘意义)最小。当拟合的函数形式是非线性函数时，就叫做非线性拟合。列文伯格-马夸尔特非线性拟合算法是本领域常见的一种求解函数值最小化的数值方法，在文献TheLevenberg-Marquardt algorithm：implementation and theory(J.More， Numerical analysis，Vol.630(1978)，pp.105-116)中有详细的介绍。 

发明内容

本发明的目的在于提供一种自动检测结霜出现的方法，基于户外载玻片图像序列，采用图像处理方法对图像序列进行处理来实现自动检测结霜。该方法能准确地检测到结霜的时间，并且操作简便，工作效率高。 

户外载玻片图像包含1个或多个载玻片，载玻片放置在户外离地面较近的地方。对这固定在户外的载玻片间隔固定的时间连续拍摄就得到了户外载玻片图像序列。 

结霜现象发生的时候，载玻片表面会发生突变，本发明通过标记载玻片在图像中的位置、计算相关序列函数、非线性拟合、相对亮度变化的判断来这些步骤来实现结霜的自动检测。 

具体方案如下： 

一种自动检测结霜的方法，利用对设置在户外的多个载玻片的图像序列I_i的处理，检测出是否结霜及结霜的时刻，其中，所述图像序列I_i由每隔固定时间采集的载玻片图像所组成，所述载玻片呈四边形，i是正整数，表示图像帧数序号，该方法具体步骤如下： 

(1)对于所述图像序列I_i中的第1帧图像I₁，标记I₁中载玻片的上底边和下底边，其中上底边的线段标记为 

下底边的线段标记为 m表示上底边，b表示下底边，k＝1，2，L，M，M为载玻片个数， 

分别为上底边两端点， 分别为下底边两端点； 

(2)当采集的图像序列超过N帧的时候，开始计算任意第t帧图像对 应的序列相关函数 

的值，获得对应的M组N个数据点 N为正整数，N≥10，t≥N+1，所述序列相关函数 

定义为： 

$g_t^k\left(\mathrm{\Δt}\right)=\frac{S\left(A_k\right)\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{I}_t(i,j)I_{t-\mathrm{\Δt}}(i,j)-\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{I}_t(i,j)\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{I}_{t-\mathrm{\Δt}}(i,j)}{\sqrt{S\left(A_k\right)\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{\left[I_t(i,j)\right]}^2-{\left[\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{I}_t(i,j)\right]}^2}\sqrt{S\left(A_k\right)\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{\left[I_{t-\mathrm{\Δt}}(i,j)\right]}^2-{\left[\underset{(i,j)\∈A_k}{\mathrm{\Σ}}{I}_{t-\mathrm{\Δt}}(i,j)\right]}^2}}$

其中，Δt＝{1，2，L，N}，I_t(i，j)表示第t帧图像中第j行第i列的灰度值，S(A_k)是载玻片所在的四边形区域A_k的面积； 

(3)对M组的每一组N个数据点 分别进行拟合，得到的任意第k组数据点 

的拟合参数记为 

拟合后的数据点 

通过下式计算得到： 

${\hat{g}}_t^k\left(\mathrm{\Δt}\right)=a_t^k\arctan (b_t^k\mathrm{\Δt}+c_t^k)+d_t^k$

(4)计算拟合后数据 

和拟合前数据 

的相关系数 

计算公式如下： 

$r_t^k=\frac{N\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{g}_t^k\left(\mathrm{\Δt}\right){\hat{g}}_t^k\left(\mathrm{\Δt}\right)-\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{g}_t^k\left(\mathrm{\Δt}\right)\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{\hat{g}}_t^k\left(\mathrm{\Δt}\right)}{\sqrt{N\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{\left[g_t^k\left(\mathrm{\Δt}\right)\right]}^2-{\left[\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{g}_t^k\left(\mathrm{\Δt}\right)\right]}^2}\sqrt{N\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{\left[{\hat{g}}_t^k\left(\mathrm{\Δt}\right)\right]}^2-{\left[\underset{\mathrm{\Δt}=1}{\overset{N}{\mathrm{\Σ}}}{\hat{g}}_t^k\left(\mathrm{\Δt}\right)\right]}^2}},$

并定义函数h_k(t)为： 

$h_k\left(t\right)=\left\{\begin{array}{cc}0,& r_t^k<\mathrm{TR}\\ t-\frac{c_t^k}{b_t^k},& r_t^k\&GreaterEqual;\mathrm{TR}\end{array}\right.,$

其中TR是阈值，TR∈[0.95，0.99]，如果h_k(t)＝0，那么返回步骤(2)继续检测下一帧图像，否则，进入步骤(5)； 

(5)设h_k′(t)为该帧图像中不为0的函数h_k(t)，h_k′(t)≠0， k′＝1，2，L，M，如果h_k′(t-w)≠0，w＝0，1，L，T，其中T是参数，为整数，T≥3，那么进入步骤(6)，否则返回步骤(2)继续检测下一帧图像； 

(6)计算载玻片区域相对亮度变化L_k′(t)，计算公式如下： 

$L_{k^{\&prime;}}\left(t\right)=\frac{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)-\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)+\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}$

其中Ω表示整个图像区域，A_k′为h_k′(t)对应的载玻片区域， 

(7)检测判断：如果L_k′(t)≥TL，那么第 

帧发生了结霜；否则无结霜，其中TL是阈值，TL∈[0.03，0.07]。 

本发明具有以下的特点： 

1、自动检测结霜现象是否发生，并且检测出结霜发生的时刻； 

2、检测过程中除了初始手工标记外，不需要人工干预，自动化程度高； 

3、对于不同形式的结霜都能很好的检测。 

附图说明

图1和图2示出了一帧户外载玻片图像，图像中包含了3块载玻片。其中图1是拍摄的户外载玻片原始图像；图2是进行手工标记后的结果。 

图3和图4示出了结霜过程中一个载玻片图像，其中图3表示第t帧图像，以及第t-1帧，t-2帧，…，t-10帧第2个载玻片图像；图4是图3中对应的一组数据对以及非线性拟合的结果。 

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。 

载玻片是呈长方形的玻璃片，户外载玻片图像序列就是把载玻片放置 在户外离地面很近的地方，在载玻片的上方放置摄像头拍摄到的图像序列。摄像头每间隔一段固定的时间(例如10分钟)就拍摄一帧灰度图像，摄像头启动以后就会一直拍摄下去，从而使得我们可以实时检测结霜的情况。本实施例以包括3个载玻片的载玻片图像序列为例。 

下面的说明将假定我们已经获得了户外载玻片图像序列，图像序列采集时间间隔为10分钟。 

下面详细说明本发明的步骤： 

(1)对于户外载玻片图像序列I_i，i是正整数，中的第1帧图像I₁，手工标记图像I₁中3个载玻片的上底边和下底边，标记的方式是勾画线段的2个端点，一共需要勾画6条线段，标记上底边的3条线段从左至右依次记为 

m表示分界线的标记，k＝1，2，3，标记下底边的3条线段从左至右依次标记为 

b表示下底边的标记，k＝1，2，3，如附图2所示。 

这样，可以得到3个载玻片所在的四边形区域： 

$A_k\{(x_{m_k,1},y_{m_k,1}),(x_{m_k,2},y_{m_k,2}),(x_{b_k,2},y_{b_k,2}),(x_{b_k,1},y_{b_k,1})\},k=\mathrm{1,2,3}.$

因为载玻片在序列图像中的位置是固定不变的，所以我们只需要手工标记一次即可。 

(2)当采集的图像序列超过N帧的时候，开始计算任意第t帧图像(t≥N+1)对应的序列相关函数 

的值，N为正整数，N≥10， 

定义为： 

$\begin{array}{c}{g}_t^k\left(\mathrm{\&Delta;t}\right)=\frac{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)I_{t-\mathrm{\&Delta;t}}(i,j)-\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_{t-\mathrm{\&Delta;t}}(i,j)}{\sqrt{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{\left[I_t(i,j)\right]}^2-{\left[\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)\right]}^2}\sqrt{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{\left[I_{t-\mathrm{\&Delta;t}}(i,j)\right]}^2-{\left[\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_{t-\mathrm{\&Delta;t}}(i,j)\right]}^2}},\\ k=\mathrm{1,2,3}\end{array}$

其中，Δt＝{1，2，L，N}，I_t(i，j)表示该第t帧图像中第j行第i列的灰度值，S(A_k)是四边形区域A_k的面积，计算公式为： 

$S\left(A_k\right)=\frac{1}{2}|\det \left(\begin{array}{cc}{x}_{m_k,1}& y_{m_k,1}\\ x_{m_k,2}& y_{m_k,2}\end{array}\right)+\det \left(\begin{array}{cc}{x}_{m_k,2}& y_{m_k,2}\\ x_{b_k,2}& y_{b_k,2}\end{array}\right)+\det \left(\begin{array}{cc}{x}_{b_k,2}& y_{b_k,2}\\ x_{b_k,1}& y_{b_k,1}\end{array}\right)+\det \left(\begin{array}{cc}{x}_{b_k,1}& y_{b_k,1}\\ x_{m_k,1}& y_{m_k,1}\end{array}\right)|$

(3)在步骤(2)中，对于第t帧图像，我们通过计算函数 

得到了对应的3组N个数据点 

Δt＝1，2，L，N，k＝1，2，3，对每一组N个数据点对进行曲线拟合，拟合函数均为y＝a arctan(bΔt+c)+d，其中a，b，c，d表示待拟合参数。采用列文伯格-马夸尔特(Levenberg-Marquardt)非线性拟合算法对每一组N个数据点 

分别进行拟合，得到的第k组数据点 

的拟合参数记为 

那么拟合后的数据点 

通过下式计算得到： 

${\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)=a_t^k\arctan (b_t^k\mathrm{\&Delta;t}+c_t^k)+d_t^k,\mathrm{\&Delta;t}=\mathrm{1,2}L,N$

(4)计算拟合后数据 和拟合前数据 的相关系数 

计算公式如下： 

$r_t^k=\frac{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right){\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)-\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right)\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)}{\sqrt{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[g_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2-{\left[\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2}\sqrt{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2-{\left[\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2}},k=\mathrm{1,2,3}$

定义函数h_k(t)： 

$h_k\left(t\right)=\left\{\begin{array}{cc}0,& r_t^k<\mathrm{TR}\\ t-\frac{c_t^k}{b_t^k},& r_t^k\&GreaterEqual;\mathrm{TR}\end{array}\right.,k=\mathrm{1,2,3}$

其中TR是阈值，TR∈[0.95，0.99]。如果h_k(t)＝0，k＝1，2，3，那么返回步骤(2)继续检测下一帧图像，否则，进入步骤(5)。 

(5)设h_k′(t)为该帧图像中不为0的函数h_k(t)，h_k′(t)≠0，k′∈{1，2，3}，如果h_k′(t-w)≠0，w＝0，1，L，T，其中T是参数，为整数，T≥3，那么进入步骤(6)，否则返回步骤(2)继续检测下一帧图像。 

(6)计算载玻片区域相对亮度变化L_k′(t)，计算公式如下： 

$L_{k^{\&prime;}}\left(t\right)=\frac{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)-\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)+\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}$

其中Ω表示整个图像区域，A_k′为h_k′(t)对应的载玻片区域。 

如果L_k′(t)≥TL，那么第 

帧发生了结霜；否则无结霜。其中TL是阈值，TL∈[0.03，0.07]。 

(7)返回步骤(2)，继续检测下一帧图像。 

Claims (4)

1.一种自动检测结霜的方法，利用对设置在户外的多个载玻片的图像序列I_i的处理，检测出是否结霜及结霜的时刻，其中，所述图像序列I_i由每隔固定时间采集的载玻片图像所组成，所述载玻片呈四边形，i是正整数，表示图像帧数序号，该方法具体步骤如下：

(1)对于所述图像序列I_i中的第1帧图像I₁，标记I₁中载玻片的上底边和下底边，其中上底边的线段标记为

下底边的线段标记为

m表示上底边，b表示下底边，k＝1，2，L，M，M为载玻片个数，

分别为上底边两端点，

分别为下底边两端点；

(2)当采集的图像序列超过N帧的时候，开始计算任意第t帧图像对应的序列相关函数

的值，获得对应的M组N个数据点

N为正整数，N≥10，t≥N+1，所述序列相关函数

定义为：

$g_t^k\left(\mathrm{\&Delta;t}\right)=\frac{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)I_{t-\mathrm{\&Delta;t}}(i,j)-\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_{t-\mathrm{\&Delta;t}}(i,j)}{\sqrt{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{\left[I_t(i,j)\right]}^2-{\left[\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_t(i,j)\right]}^2}\sqrt{S\left(A_k\right)\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{\left[I_{t-\mathrm{\&Delta;t}}(i,j)\right]}^2-{\left[\underset{(i,j)\&Element;A_k}{\mathrm{\&Sigma;}}{I}_{t-\mathrm{\&Delta;t}}(i,j)\right]}^2}}$

其中，Δt＝{1，2，L，N}，I_t(i，j)表示第t帧图像中第j行第i列的灰度值，S(A_k)是载玻片所在的四边形区域A_k的面积；

(3)对M组的每一组N个数据点分别进行拟合，得到的任意第k组数据点

的拟合参数记为

拟合后的数据点

通过下式计算得到：

${\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)=a_t^k\arctan (b_t^k\mathrm{\&Delta;t}+c_t^k)+d_t^k$

(4)计算拟合后数据

和拟合前数据的相关系数

计算公式如下：

$r_t^k=\frac{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right){\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)-\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right)\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)}{\sqrt{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[g_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2-{\left[\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2}\sqrt{N\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2-{\left[\underset{\mathrm{\&Delta;t}=1}{\overset{N}{\mathrm{\&Sigma;}}}{\hat{g}}_t^k\left(\mathrm{\&Delta;t}\right)\right]}^2}},$

并定义函数h_k(t)为：

$h_k\left(t\right)=\left\{\begin{array}{cc}0,& r_t^k<\mathrm{TR}\\ t-\frac{c_t^k}{b_t^k},& r_t^k\&GreaterEqual;\mathrm{TR}\end{array}\right.,$

其中TR是阈值，TR∈[0.95，0.99]，如果h_k(t)＝0，那么返回步骤(2)继续检测下一帧图像，否则，进入步骤(5)；

(5)设h_k′(t)为该帧图像中不为0的函数h_k(t)，h_k′(t)≠0，k′＝1，2，L，M，如果h_k′(t-w)≠0，w＝0，1，L，T，其中T是参数，为整数，T≥3，那么进入步骤(6)，否则返回步骤(2)继续检测下一帧图像；

(6)计算载玻片区域相对亮度变化L_k′(t)，计算公式如下：

$L_{k^{\&prime;}}\left(t\right)=\frac{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)-\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}{\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_t(i,j)+\underset{(i,j)\&Element;A_{k^{\&prime;}}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)+\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_{t-N}(i,j)-\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{I}_t(i,j)}$

其中Ω表示整个图像区域，A_k′为h_k′(t)对应的载玻片区域，

(7)检测判断：如果L_k′(t)≥TL，那么第

帧发生了结霜；否则无结霜，其中TL是阈值，TL∈[0.03，0.07]。

2.根据权利要求1所述的方法，其特征在于，所述的载玻片的四边形区域A_k的面积S(A_k)的计算公式为：

$S\left(A_k\right)=\frac{1}{2}|\det \left(\begin{array}{cc}{x}_{m_k,1}& y_{m_k,1}\\ x_{m_k,2}& y_{m_k,2}\end{array}\right)+\det \begin{array}{c}\left(\begin{array}{cc}{x}_{m_k,2}& y_{m_k,2}\\ x_{b_k,2}& y_{b_k,2}\end{array}\right)+\det \end{array}\left(\begin{array}{cc}{x}_{b_k,2}& y_{b_k,2}\\ x_{b_k,1}& y_{b_k,1}\end{array}\right)+\det \left(\begin{array}{cc}{x}_{b_k,1}& y_{b_k,1}\\ x_{m_k,1}& y_{m_k,1}\end{array}\right)|$

3.根据权利要求1或2所述的方法，其特征在于，拟合函数为y＝a arctan(bΔt+c)+d，其中a，b，c，d表示待拟合参数。

4.根据权利要求1-3之一所述的方法，其特征在于，所述拟合采用列文伯格-马夸尔特(Levenberg-Marquardt)非线性拟合算法。

Images