显示标签为“算法”的博文。显示所有博文
显示标签为“算法”的博文。显示所有博文

2010年5月12日星期三

resenham算法是计算机图形学典型的直线光栅化算法。

resenham算法是计算机图形学典型的直线光栅化算法。

  • 从另一个角度看直线光栅化显示算法的原理
    • 由直线的斜率确定选择在x方向或y方向上每次递增(减)1个单位,另一变量的递增(减)量为0或1,它取决于实际直线与最近光栅网格点的距离,这个距离的最大误差为0.5。

     

  • 1)Bresenham的基本原理

     

    • 假定直线斜率k在0~1之间。此时,只需考虑x方向每次递增1个单位,决定y方向每次递增0或1。

      直线当前点为(xi,y)
      直线当前光栅点为(xi,yi)

      下一个直线的点应为(xi+1,y+k)
      下一个直线的光栅点
      或为右光栅点(xi+1,yi)(y方向递增量0)
      或为右上光栅点(xi+1,yi+1)(y方向递增量1)

      记直线与它垂直方向最近的下光栅点的误差为d,有:d=(y+k)–yi,且

      0≤d≤1
      当d<0.5:下一个象素应取右光栅点(xi+1,yi)>

      如果直线的(起)端点在整数点上,误差项d的初值:d0=0,
      x坐标每增加1,d的值相应递增直线的斜率值k,即:d=d + k。
      一旦d≥1,就把它减去1,保证d的相对性,且在0-1之间。

      令e=d-0.5,关于d的判别式和初值可简化成:

      e的初值e0= -0.5,增量亦为k;
      e<0时,取当前象素(xi,yi)的右方象素(xi+1,yi);>0时,取当前象素(xi,yi)的右上方象素(xi+1,yi+1);
      e=0时,可任取上、下光栅点显示。

      Bresenham算法的构思巧妙:它引入动态误差e,当x方向每次递增1个单位,可根据e的符号决定y方向每次递增 0 或 1。

      e<0,y方向不递增>0,y方向递增1
      x方向每次递增1个单位,e = e + k

      因为e是相对量,所以当e>0时,表明e的计值将进入下一个参考点(上升一个光栅点),此时须:e = e - 1

       

  • 2)Bresenham算法的实施——Rogers 版

     

    • 通过(0,0)的所求直线的斜率大于0.5,它与x=1直线的交点离y=1直线较近,离y=0直线较远,因此取光栅点(1,1)比(1,0)更逼近直线;
      如果斜率小于0.5,则反之;
      当斜率等于0.5,没有确定的选择标准,但本算法选择(1,1)

      程序

       

      • //Bresenham's line resterization algorithm for the first octal.
        //The line end points are (xs,ys) and (xe,ye) assumed not equal.
        // Round is the integer function.
        // x,y, ∆x, ∆y are the integer, Error is the real.
        //initialize variables
        x=xs
        y=ys
        ∆x = xe -xs
        ∆y = ye -ys
        //initialize e to compensate for a nonzero intercept
        Error =∆y/∆x-0.5
        //begin the main loop
        for i=1 to ∆x
        WritePixel (x, y, value)
        if (Error ≥0) then
        y=y+1
        Error = Error -1
        end if
        x=x+1
        Error = Error +∆y/∆x
        next i
        finish

       

  • 3)整数Bresenham算法

     

    • 上述Bresenham算法在计算直线斜率和误差项时要用到浮点运算和除法,采用整数算术运算和避免除法可以加快算法的速度。

      由于上述Bresenham算法中只用到误差项(初值Error =∆y/∆x-0.5)的符号

      因此只需作如下的简单变换:

      NError = 2*Error*∆x

      即可得到整数算法,这使本算法便于硬件(固件)实现。

      程序

       

      • //Bresenham's integer line resterization algorithm for the first octal.
        //The line end points are (xs,ys) and (xe,ye) assumed not equal. All variables are assumed integer.
        //initialize variables
        x=xs
        y=ys
        ∆x = xe -xs
        ∆y = ye -ys
        //initialize e to compensate for a nonzero intercept
        NError =2*∆y-∆x //Error =∆y/∆x-0.5
        //begin the main loop
        for i=1 to ∆x
        WritePixel (x, y)
        if (NError >=0) then
        y=y+1
        NError = NError –2*∆x //Error = Error -1
        end if
        x=x+1
        NError = NError +2*∆y //Error = Error +∆y/∆x
        next i
        finish

       

  • 4)一般Bresenham算法

     

    • 要使第一个八卦的Bresenham算法适用于一般直线,只需对以下2点作出改造:
      当直线的斜率|k|>1时,改成y的增量总是1,再用Bresenham误差判别式确定x变量是否需要增加1;
      x或y的增量可能是“+1”或“-1”,视直线所在的象限决定。

      程序

       

      • //Bresenham's integer line resterization algorithm for all quadrnts
        //The line end points are (xs,ys) and (xe,ye) assumed not equal. All variables are assumed integer.
        //initialize variables
        x=xs
        y=ys
        ∆x = abs(xe -xs) //∆x = xe -xs
        ∆y = abs(ye -ys) //∆y = ye -ys
        sx = isign(xe -xs)
        sy = isign(ye -ys)
        //Swap ∆x and ∆y depending on the slope of the line.
        if ∆y>∆x then
        Swap(∆x,∆y)
        Flag=1
        else
        Flag=0
        end if
        //initialize the error term to compensate for a nonezero intercept
        NError =2*∆y-∆x
        //begin the main loop
        for i=1 to ∆x
        WritePixel(x, y , value)
        if (Nerror>=0) then
        if (Flag) then //∆y>∆x,Y=Y+1
        x=x+sx
        else
        y=y+sy
        end if // End of Flag
        NError = NError –2*∆x
        end if // End of Nerror
        if (Flag) then //∆y>∆x,X=X+1
        y=y+sy
        else
        x=x+sx
        end if
        NError = NError +2*∆y
        next i
        finish

       

  • 例子

中点画圆算法

为了能以任意点为圆心画圆,我们可以把圆心先设为视点(相当于于将其平移到坐标原点),然后通过中点法扫描转换后,再恢复原来的视点(相当于将圆心平移回原来的位置)。

圆心位于原点的圆有四条对称轴x=0,y=0,x=yx=-y,从而圆上一点(x,y),可得到其关于四条对称轴的七个对称点,这称为八对称性,下面的函数就用来显示(x,y)及其七个对称点.
200772802.jpg


void CirclePoints(int x,int y,long color,CDC *pDC)

{

//第1象限

pDC->SetPixel(x,y,color);

pDC->SetPixel(y,x,color);

//第2象限

pDC->SetPixel(-x,y,color);

pDC->SetPixel(-y,x,color);

//第3象限

pDC->SetPixel(-y,-x,color);

pDC->SetPixel(-x,-y,color);

//第4象限

pDC->SetPixel(x,-y,color);

pDC->SetPixel(y,-x,color);

}

中点画圆算法就是每部单位间隔取样并且计算离圆最近的位置。在继续之前,我这里补充一个关于圆对称性的知识点,通过在圆中计算考虑使用对称性计算开销可以减小到原来的1/8。对称性质原理:(1)圆是满足x轴对称的,这样只需要计算原来的1/2点的位置;(2)圆是满足y轴对称的,这样只需要计算原来的1/2点的位置;(3)圆是满足y = x or y = -x轴对称的,这样只需要计算原来的1/2点的位置;通过上面三个性质分析得知,对于元的计算只需要分析其中1/8的点即可。例如:分析出来目标点(x,y)必然存在(x,-y),(-x,y),(-x,-y),(y,x),(y,-x),(-y,x),(-y,-x)的另外7个点。关于中心画圆算法,通过计算x = 0到 x = y的1/8圆的范围,然后通过对称原理得到其他7/8个点的信息。这里和Bresenham算法有很多相似之处,同样有一个决定下一个位置的关键值P来做权衡处理。在中点画圆算法中,通过平移的方法将假设圆心在坐标原点,然后计算,最后再平移到真实原心位置。 如果我们构造函数 F(x,y)=x2+y2-R2,则对于圆上的点有F(x,y)=0,对于圆外的点有F(x,y)>0,对于圆内的点F(x,y)<0 d="F(M)=" d="F(xp+2,yp-0.5)=" r2="">

d≥0,则应取P2为下一象素,而且下一象素的判别式为

d=F(xp+2,yp-1.5)=(xp+2)2+(yp-1.5)2-R2=d+2(xp-yp)+5

我们这里讨论的第一个象素是(0,R),判别式d的初始值为:

d0=F(1,R-0.5)=1.25-R

200772803.jpg

中点画圆算法内容:

1,输入圆心位置和圆的半径,得到圆周上的第一个点Point1;

(假设起始点为坐标原点,后面将通过坐标平移来处理非圆心在圆点)

2,计算决策关键参数的初始值,P = 5/4 - r;

3,在每个Xn的位置,从n = 0开始,更具决策值P来判断:

如果P<0,下一个点的位置为(Xn+1,Yn);

并且执行P = P + 2*x+3;

如果P>=0,下一个点的位置为(Xn+1,Yn-1);

并且执行P = P + 2.0*(x-y)+5;

4,通过对称原理计算其他7个对称相关点;

5,移动坐标到圆心点(x1,y1)

X = X + x1;

Y = Y + y1;

6,如果X重复执行35的步骤,否则结束该算法

程序如下:

void Circle::Draw(CDC *pDC)
{//中点算法画圆
int x,y;
double p;
pDC->SetViewportOrg(pMid);
x=0;
y=radis;
p=1.25-radis;
while(x<=y+1) { CirclePoints(x,y,m_lPenColor,pDC); x++; if(p>=0)
{
y--;
p+=2.0*(x-y)+5;
}
else
p+=2*x+3;
}
pDC->SetViewportOrg(0,0);
}

bresenham画圆算法

bresenham画圆算法

中点画圆算法在一个方向上取单位间隔,在另一个方向的取值由两种可能取值的中点离圆的远近而定。实际处理中,用决策变量的符号来确定象素点的选择,因此算法效率较高。

  Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页一、中点画圆算法描述

  设要显示圆的圆心在原点(0,0),半径为R,起点在(0,R)处,终点在(Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页)处,顺时针生成八分之一圆,利用对称性扫描转换全部圆。

  为了应用中点画圆法,我们定义一个圆函数

F(x,y)=x2+y2-R2(2-19)

  任何点(x,y)的相对位置可由圆函数的符号来检测:

F(x,y)Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页<0 style="line-height: 22px; margin-top: 0px; margin-right: 0px; margin-bottom: 10px; margin-left: 0px; padding-top: 0px; padding-right: 0px; padding-bottom: 0px; padding-left: 0px; ">

=0 点(x,y)位于数学圆上

>0 点(x,y)位于数学圆外

(2-20)

  如下图所示,图中有两条圆弧A和B,假定当前取点为Pi(xi,yi),如果顺时针生成圆,那么下一点只能取正右方的点E(xi+1,yi)或右下方的点SE(xi+1,yi-1)两者之一。

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页

中点画线算法

  假设M是E和SE的中点,即 Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页,则:  1、当F(M)<0时,m在圆内(圆弧a),这说明点e距离圆更近,应取点e作为下一象素点;

  2、当F(M)>0时,M在圆外(圆弧B),表明SE点离圆更近,应取SE点;

  3、当F(M)=0时,在E点与SE点之中随便取一个即可,我们约定取SE点。

  Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页二、中点画圆算法思想

  因此,我们用中点M的圆函数作为决策变量di,同时用增量法来迭代计算下一个中点M的决策变量di+1。

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页(2-21)

  下面分两种情况来讨论在迭代计算中决策变量di+1的推导。

  1、见图(a),若di<0,则选择e点,接着下一个中点就是Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页,这时新的决策变量为:

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页(2-22)

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页

(a)(di<0)>

  式(2-22)减去(2-21)得:

di+1=di+2xi+3(2-23)

  2、见图(b),若di≥0,则选择SE点,接着下一个中点就是Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页 ,这时新的决策变量为:

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页(2-24)

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页

(b)(di≥0) 中点画线算法

  式(2-24)减去(2-21)得:

di+1=di+2(xi-yi)+5(2-25)

  我们利用递推迭代计算这八分之一圆弧上的每个点,每次迭代需要两步处理:

   (1)用前一次迭代算出的决策变量的符号来决定本次选择的点。

   (2)对本次选择的点,重新递推计算得出新的决策变量的值。

  剩下的问题是计算初始决策变量d0,如下图所示。对于初始点(0,R),顺时针生成八分之一圆,下一个中点M的坐标是Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页 ,所以:

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页(2-26)

Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页

生成圆的初始条件和圆的生成方向

  Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页三、中点画圆算法实现

  1、输入:圆半径r、圆心(x0,y0);

  2、确定初值:x=0,y=r、d=5/4-r;

  3、While(x<=y)

   {

    ·利用八分对称性,用规定的颜色color画八个象素点(x,y);

    · 若d≥0

      {

       y=y-1; //wind:个人觉得这句应该置于下句

       d=d+2(x-y)+5);

      }

     否则

       d=d+2x+3;

    ·x=x+1;

   }

  Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页五、中点画圆算法完善

  在上述算法中,使用了浮点数来表示决策变量d。为了简化算法,摆脱浮点数,在算法中全部使用整数,我们使用e=d-1/4代替d。显然,初值d=5/4-r对应于e=1-r。决策变量d<0对应于e<-1/4。算法中其它与d有关的式子可把d直接换成e。又由于e的初值为整数,且在运算过程中的迭代值也是整数,故e始终是整数,所以e<-1/4等价于e<0。因此,可以写出完全用整数实现的中点画圆算法。

  要求:写出用整数实现的中点画圆算法程序,并上机调试,观看运行结果。

  Bresenham画圆和椭圆程序 - 随风倒上 - 赵志刚--Greddy的个人主页六、中点画圆算法程序

void MidpointCircle(int x0,int y0,int r,int color)

{

 int x,y;

 float d;

 x=0;

 y=r;

 d=5.0/4-r;

 while(x<=y)

 {

  putdot(x0,y0,x,y,color);

  if(d<0)

  d+=x*2.0+3;

  else

  {

   d+=2.0*(x-y)+5;

   y--;

  }

  x++;

 }

}

putdot(x0,y0,x,y,color)

{

 putpixel(x0+x,y0+y,color);

 putpixel(x0+x,y0-y,color);

 putpixel(x0-x,y0+y,color);

 putpixel(x0-x,y0-y,color);

 putpixel(x0+y,y0+x,color);

 putpixel(x0+y,y0-x,color);

 putpixel(x0-y,y0+x,color);

 putpixel(x0-y,y0-x,color);

}

2008年7月1日星期二

模糊数学(Fuzzy mathematics)及其应用

模糊数学(Fuzzy mathematics)及其应用
一 :引言
有一个古老的希腊悖论,是这样说的:“一粒种子肯定不叫一堆,两粒也不是,三粒也不是……另一方面,所有的人都同意,一亿粒种子肯定叫一堆。那么,适当的界限在哪里?我们能不能说,123585粒种子不叫一堆而123586粒就构成一堆?”
确 实,“一粒”和“一堆”是有区别的两个概念。但是,它们的区别是逐渐的,而不是突变的,两者之间并不存在明确的界限。换句话说,“一堆”这个概念带有某种 程度的模糊性。类似的概念,如“年老”、“高个子”、“年轻人”、“很大”、“聪明”、“漂亮的人”、“价廉物美”等等,不胜枚举。
经 典集合论中,在确定一个元素是否属于某集合时,只能有两种回答:“是”或者“不是”。我们可以用两个值0或1加以描述,属于集合的元素用1表示,不属于集 合的元素用0表示。然而上面提到的“年老”、“高个子”、“年轻人”、“很大”、“聪明”、“漂亮的人”、“价廉物美” 等情况要复杂得多。假如规定身高1.8米算属于高个子范围,那么,1.79米的算不算?照经典集合论的观点看:不算。但这似乎很有些悖于情理。如果用一个 圆,以圆内和圆周上的点表示集A,而且圆外的点表示不属于A。A的边界显然是圆周。这是经典集合的图示。现在,设想将高个子的集合用图表示,则它的边界将 是模糊的,即可变的。因为一个元素(例如身高1.75米的人)虽然不是100%的高个子,却还算比较高,在某种程度上属于高个子集合。这时一个元素是否属 于集合,不能光用0和1两个数字表示,而可以取0和1之间的任何实数。例如对1.75米的身高,可以说具有70%属于高个子集合的程度。这样做似乎罗嗦, 但却比较合乎实际。
精 确和模糊,是一对矛盾。根据不同情况有时要求精确,有时要求模糊。比如打仗,指挥员下达命令:“拂晓发起总攻。”这就乱套了。这时,一定要求精确:“×月 ×日清晨六时正发起总攻。”我们在一些旧电影中还能看到各个阵地的指挥员在接受命令前对对表的镜头,生怕出个半分十秒的误差。但是,物极必反。如果事事要 求精确,人们就简直无法顺利的交流思想——两人见面,问:“你好吗?”可是,什么叫“好”,又有谁能给“好”下个精确的定义?
有些现象本质上就是模糊的,如果硬要使之精确,自然难以符合实际。例如,考核学生成绩,规定满60分为合格。但是,59分和60分之间究竟有多大差异,仅据1分之差来区别及格和不及格,其根据是不充分的。
不 仅普遍存在着边界模糊的集合,就是人类的思维,也带有模糊的特色。有些现象是精确的,但是,适当的模糊化可能使问题得到简化,灵活性大为提高。例如,在地 里摘玉米,若要找一个最大的,那很麻烦,而且近乎迂腐。我们必须把玉米地里所有的玉米都测量一下,再加以比较才能确定。它的工作量跟玉米地面积成正比。土 地面积越大,工作越困难。然而,只要稍为改变一下问题的提法:不要求找最大的玉米,而是找比较大的,即按通常的说法,到地里摘个大玉米。这时,问题从精确 变成了模糊,但同时也从不必要的复杂变成意外的简单,挑不多的几个就可以满足要求。工作量甚至跟土地无关。因此,过分的精确实际成了迂腐,适当的模糊反而 灵活。
显 然,玉米的大小,取决于它的长度、体积和重量 。大小虽是模糊概念,但长度、体积、重量等在理论上都可以是精确的。然而,人们在实际判断玉米大小时,通常并不需要测定这些精确值。同样,模糊的“堆”的 概念是建立在精确的“粒”的基础上,而人们在判断眼前的东西叫不叫一堆时,从来不用去数“粒”。有时,人们把模糊性看成一种物理现象。近的东西看得清,远 的东西看不清,一般的说,越远越模糊。但是,也有例外的情况:站在海边,海岸线是模糊的;从高空向下眺望,海岸线却显得十分清晰。太高了,又模糊。精确与 模糊,有本质区别,但又有内在联系,两者相互矛盾、相互依存也可相互转化。所以,精确性的另一半是模糊。
对 模糊性的讨论,可以追溯得很早。20世纪的大哲学家罗素(B.Russel)在1923年一篇题为《含糊性》(Vagueness)的论文里专门论述过我 们今天称之为“模糊性”的问题(严格地说,两者稍有区别),并且明确指出:“认为模糊知识必定是靠不住的,这种看法是大错特错的。”尽管罗素声名显赫,但 这篇发表在南半球哲学杂志的文章并未引起当时学术界对模糊性或含糊性的很大兴趣。这并非是问题不重要,也不是因为文章写得不深刻,而是“时候未到”。罗素 精辟的观点是超前的。长期以来,人们一直把模糊看成贬义词,只对精密与严格充满敬意。20世纪初期社会的发展,特别是科学技术的发展,还未对模糊性的研究 有所要求。事实上,模糊性理论是电子计算机时代的产物。正是这种十分精密的机器的发明与广泛应用,使人们更深刻地理解了精密性的局限,促进了人们对其对立 面或者说它的“另一半”——模糊性的研究。
集 合是现代数学的基础,模糊集合一提出,“模糊”观念也渗透到许多数学分支。模糊数学的发展速度也是相当快的。从发表的论文看,几乎是指数般的增长。模糊数 学的研究可分三个方面:一是研究模糊数学的理论,以及它和精确数学、统计数学的关系;二是研究模糊语言和模糊逻辑;三是研究模糊数学的应用。在模糊数学的 研究中,目前已有模糊拓扑学、模糊群论、模糊凸论、模糊概率、模糊环论等分支。虽然模糊数学是一门新兴学科,但它已初步应用于自动控制、模式识别、系统理 论、信系检索、社会科学、心理学、医学和生物学等方面。将来还可能出现模糊逻辑电路、模糊硬件、模糊软件和模糊固件,出现能和人用自然语言对话、更接近于 人的智能的新的一类计算机。所以,模糊数学将越来越显示出它的巨大生命力。
二: 模糊数学的产生
二 十世纪六十年代,产生了模糊数学这门新兴学科。 现代数学是建立在集合论的基础上。集合论的重要意义就一个侧面看,在与它把数学的抽象能力延伸到人类认识 过程的深处。一组对象确定一组属性,人们可以通过说明属性来说明概念(内涵),也可以通过指明对象来说明它。符合概念的那些对象的全体叫做这个概念的外 延,外延其实就是集合。从这个意义上讲,集合可以表现概念,而集合论中的关系和运算又可以表现判断和推理,一切现实的理论系统都一可能纳入集合描述的数学 框架。
但是,数学的发展也是阶段性的。经典集合论只能把自己的表现力限制在那些有明确外延的概念和事物上,它明确地限定:每个集合都必须 由明确的元素构成,元素对集合的隶属关系必须是明确的,决不能模棱两可。对于那些外延不分明的概念和事物,经典集合论是暂时不去反映的,属于待发展的范 畴。
在较长时间里,精确数学及随机数学在描述自然界多种事物的运动规律中,获得显著效果。但是,在客观世界中还普遍存在着大量的模糊现象。以前人们回避它,但是,由于现代科技所面对的系统日益复杂,模糊性总是伴随着复杂性出现。
各门学科,尤其是人文、社会学科及其它“软科学”的数学化、定量化趋向把模糊性的数学处理问题推向中心地位。更重要的是,随着电子计算机、控制论、系统科学的迅速发展,要使计算机能像人脑那样对复杂事物具有识别能力,就必须研究和处理模糊性。
我们研究人类系统的行为,或者处理可与人类系统行为相比拟的复杂系统,如航天系统、人脑系统、社会系统等,参数和变量甚多,各种因素相互交错,系统很复杂,它的模糊性也很明显。从认识方面说,模糊性是指概念外延的不确定性,从而造成判断的不确定性。
在 日常生活中,经常遇到许多模糊事物,没有分明的数量界限,要使用一些模糊的词句来形容、描述。比如,比较年轻、高个、大胖子、好、漂亮、善、热、远……。 在人们的工作经验中,往往也有许多模糊的东西。例如,要确定一炉钢水是否已经炼好,除了要知道钢水的温度、成分比例和冶炼时间等精确信息外,还需要参考钢 水颜色、沸腾情况等模糊信息。因此,除了很早就有涉及误差的计算数学之外,还需要模糊数学。
人与计算机相比,一般来说,人脑具有处理模糊 信息的能力,善于判断和处理模糊现象。但计算机对模糊现象识别能力较差,为了提高计算机识别模糊现象的能力,就需要把人们常用的模糊语言设计成机器能接受 的指令和程序,以便机器能像人脑那样简洁灵活的做出相应的判断,从而提高自动识别和控制模糊现象的效率。这样,就需要寻找一种描述和加工模糊信息的数学工 具,这就推动数学家深入研究模糊数学。所以,模糊数学的产生是有其科学技术与数学发展的必然性。
模糊数学的研究内容
1965年,美国控制论专家、数学家查德发表了论文《模糊集合》,标志着模糊数学这门学科的诞生。
模糊数学的研究内容主要有以下三个方面:
第 一,研究模糊数学的理论,以及它和精确数学、随机数学的关系。察德以精确数学集合论为基础,并考虑到对数学的集合概念进行修改和推广。他提出用“模糊集合 ”作为表现模糊事物的数学模型。并在“模糊集合”上逐步建立运算、变换规律,开展有关的理论研究,就有可能构造出研究现实世界中的大量模糊的数学基础,能 够对看来相当复杂的模糊系统进行定量的描述和处理的数学方法。
在模糊集合中,给定范围内元素对它的隶属关系不一定只有“是”或“否”两种 情况,而是用介于0和1之间的实数来表示隶属程度,还存在中间过渡状态。比如“老人”是个模糊概念,70岁的肯定属于老人,它的从属程度是 1,40岁的人肯定不算老人,它的从属程度为 0,按照查德给出的公式,55岁属于“老”的程度为0.5,即“半老”,60岁属于“老”的程度0.8。查德认为,指明各个元素的隶属集合,就等于指定了 一个集合。当隶属于0和1之间值时,就是模糊集合。
第二,研究模糊语言学和模糊逻辑。人类自然语言具有模糊性,人们经常接受模糊语言与模糊信息,并能做出正确的识别和判断。
为了实现用自然语言跟计算机进行直接对话,就必须把人类的语言和思维过程提炼成数学模型,才能给计算机输入指令,建立和是的模糊数学模型,这是运用数学方法的关键。查德采用模糊集合理论来建立模糊语言的数学模型,使人类语言数量化、形式化。
如 果我们把合乎语法的标准句子的从属函数值定为1,那么,其他文法稍有错误,但尚能表达相仿的思想的句子,就可以用以0到1之间的连续数来表征它从属于“正 确句子”的隶属程度。这样,就把模糊语言进行定量描述,并定出一套运算、变换规则。目前,模糊语言还很不成熟,语言学家正在深入研究。
人们的思维活动常常要求概念的确定性和精确性,采用形式逻辑的排中律,既非真既假,然后进行判断和推理,得出结论。现有的计算机都是建立在二值逻辑基础上的,它在处理客观事物的确定性方面,发挥了巨大的作用,但是却不具备处理事物和概念的不确定性或模糊性的能力。
为了使计算机能够模拟人脑高级智能的特点,就必须把计算机转到多值逻辑基础上,研究模糊逻辑。目前,模糊罗基还很不成熟,尚需继续研究。
第 三,研究模糊数学的应用。模糊数学是以不确定性的事物为其研究对象的。模糊集合的出现是数学适应描述复杂事物的需要,查德的功绩在于用模糊集合的理论找到 解决模糊性对象加以确切化,从而使研究确定性对象的数学与不确定性对象的数学沟通起来,过去精确数学、随机数学描述感到不足之处,就能得到弥补。在模糊数 学中,目前已有模糊拓扑学、模糊群论、模糊图论、模糊概率、模糊语言学、模糊逻辑学等分支。
模糊数学的应用
模糊数学是一门新兴学 科,它已初步应用于模糊控制、模糊识别、模糊聚类分析、模糊决策、模糊评判、系统理论、信息检索、医学、生物学等各个方面。在气象、结构力学、控制、心理 学等方面已有具体的研究成果。然而模糊数学最重要的应用领域是计算机职能,不少人认为它与新一代计算机的研制有密切的联系。
目前,世界上 发达国家正积极研究、试制具有智能化的模糊计算机,1986年日本山川烈博士首次试制成功模糊推理机,它的推理速度是1000万次/秒。1988年,我国 汪培庄教授指导的几位博士也研制成功一台模糊推理机——分立元件样机,它的推理速度为1500万次/秒。这表明我国在突破模糊信息处理难关方面迈出了重要 的一步。
模糊数学还远没有成熟,对它也还存在着不同的意见和看法,有待实践去检验。
四:模糊数学的主要应用
1 模糊数学自身的理论研究进展迅速。我国模糊数学自身的理论研究仍占模糊数学及其应用学科的主导地位,所取得的研究成果在《模糊数学》、《模糊系统与数学》 等数十种学术期刊和全国高校学报中经常可见,模糊聚类分析理论、模糊神经网络理论和各种新的模糊定理及算法不断取得进展。
2.模糊数学目前在自动控制技术领域仍然得到最广泛的应用,所涉及的技术复杂繁多,从微观到宏观、从地下到太空无所不有,在机器人实时控制、电磁元件自适应控制、各种物理及力学参数反馈控制、逻辑控制等高新技术中均成功地应用了模糊数学理论和方法。
3.模糊数学在计算机仿真技术、多媒体辨识等领域的应用取得突破性进展,如图像和文字的自动辨识、自动学习机、人工智能、音频信号辨识与处理等领域均借助了模糊数学的基本原理和方法。
4. 模糊聚类分析理论和模糊综合评判原理等更多地被应用于经济管理、环境科学、安全与劳动保护等领域,如房地价格、期货交易、股市情报、资产评估、工程质量分 析、产品质量管理、可行性研究、人机工程设计、环境质量评价、资源综合评价、各种危险性预测与评价、灾害探测等均成功地应用了模糊数学的原理和方法。
5.地 矿、冶金、建筑等传统行业在处理复杂不确定性问题中也成功地应用了模糊数学的原理和方法,从而使过去凭经验和类比法等处理工程问题的传统做法转向数学化、 科学化,如矿床预测、矿体边界确定、油水气层的识别、采矿方法设计参数选择、冶炼工艺自动控制与优化、建筑物结构设计等都有应用模糊数学的成功实践。
6. 我国医药、生物、农业、文化教育、体育等过去看似与数学无缘的学科也开始应用模糊数学的原理和方法,如计算机模糊综合诊断、传染病控制与评估、人体心理及 生理特点分析、家禽孵养、农作物品种选择与种植、教学质量评估、语言词义查找、翻译辨识等均有一些应用模糊数学的实践,并取得很好效果。
模糊数学目前在自动控制技术领域仍然得到最广泛的应用,所涉及的技术复杂繁多,从微观到宏观、从地下到太空无所不有,在机器人实时控制、电磁元件自适应控制、各种物理及力学参数反馈控制、逻辑控制等高新技术中均成功地应用了模糊数学理论和方法。
模糊数学在计算机仿真技术、多媒体辨识等领域的应用取得突破性进展,如图像和文字的自动辨识、自动学习机、人工智能、音频信号辨识与处理等领域均借助了模糊数学的基本原理和方法。
模 糊聚类分析理论和模糊综合评判原理等更多地被应用于经济管理、环境科学、安全与劳动保护等领域,如房地价格、期货交易、股市情报、资产评估、工程质量分 析、产品质量管理、可行性研究、人机工程设计、环境质量评价、资源综合评价、各种危险性预测与评价、灾害探测等均成功地应用了模糊数学的原理和方法。
地 矿、冶金、建筑等传统行业在处理复杂不确定性问题中也成功地应用了模糊数学的原理和方法,从而使过去凭经验和类比法等处理工程问题的传统做法转向数学化、 科学化,如矿床预测、矿体边界确定、油水气层的识别、采矿方法设计参数选择、冶炼工艺自动控制与优化、建筑物结构设计等都有应用模糊数学的成功实践。
我 国医药、生物、农业、文化教育、体育等过去看似与数学无缘的学科也开始应用模糊数学的原理和方法,如计算机模糊综合诊断、传染病控制与评估、人体心理及生 理特点分析、家禽孵养、农作物品种选择与种植、教学质量评估、语言词义查找、翻译辨识等均有一些应用模糊数学的实践,并取得很好效果。
五:
但是一些概率论 学者认为模糊数学不过是概率论的一个应用而已。一些搞理论数学的人说这不是数学。搞应用的人则说道理说的很好,但真正的实际效果没有。然而,国际著名的应 用数学家考夫曼(A.Kauffman)教授在访华时说:“他们的攻击是毫无道理的,不必管人家说什么,我们努力去做就是。”
模糊数学是一门崭新的数学学科,它的产生不仅拓广了经典数学的基础,而且是使计算机科 学向人们的自然机理方面发展的重大突破。它在科学技术、经济发展和社会学等问题的广泛应用领域中显示了巨大的力量。它虽然只有二十多年的历史,但已被国内 外数学界以及信息、系统、计算机和自动控制科学、人员的普遍关注,它是正在迅速发展中的有着广阔应用前景的一门崭新学科。

2008年6月4日星期三

Kalman filter

1 什么是卡尔曼滤波器
What is the Kalman Filter?

在学习卡尔曼滤波器之前,首先看看为什么叫卡尔曼。跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人!

卡尔曼全名Rudolf Emil Kalman,匈牙利数学家,1930年出生于匈牙利首都布达佩斯。19531954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。如果对这编论文有兴趣,可以到这里的地址下载: http://www.cs.unc.edu/~welch/media/pdf/Kalman1960.pdf

简单来说,卡尔曼滤波器是一个“optimal recursive data processing algorithm(最优化自回归数据处理算法)。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。

2
.卡尔曼滤波器的介绍
Introduction to the Kalman Filter

为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。

在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。

假设我们要研究的对象是一个房间的温度。根据你的经验判断,这个房间的温度是恒定的,也就是下一分钟的温度等于现在这一分钟的温度(假设我们用一分钟来做时间单位)。假设你对你的经验不是100%的相信,可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声(White Gaussian Noise),也就是这些偏差跟前后时间是没有关系的而且符合高斯分配(Gaussian Distribution)。另外,我们在房间里放一个温度计,但是这个温度计也不准确的,测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。

好了,现在对于某一分钟我们有两个有关于该房间的温度值:你根据经验的预测值(系统的预测值)和温度计的值(测量值)。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。

假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值,来预测k时刻的温度。因为你相信温度是恒定的,所以你会得到k时刻的温度预测值是跟k-1时刻一样的,假设是23度,同时该值的高斯噪声的偏差是5度(5是这样得到的:如果k-1时刻估算出的最优温度值的偏差是3,你对自己预测的不确定度是4度,他们平方相加再开方,就是5)。然后,你从温度计那里得到了k时刻的温度值,假设是25度,同时该值的偏差是4度。

由于我们用于估算k时刻的实际温度有两个温度值,分别是23度和25度。究竟实际温度是多少呢?相信自己还是相信温度计呢?究竟相信谁多一点,我们可以用他们的covariance来判断。因为Kg^2=5^2/(5^2+4^2),所以Kg=0.78,我们可以估算出k时刻的实际温度值是:23+0.78*(25-23)=24.56度。可以看出,因为温度计的covariance比较小(比较相信温度计),所以估算出的最优温度值偏向温度计的值。

现在我们已经得到k时刻的最优温度值了,下一步就是要进入k+1时刻,进行新的最优估算。到现在为止,好像还没看到什么自回归的东西出现。对了,在进入k+1时刻之前,我们还要算出k时刻那个最优值(24.56度)的偏差。算法如下:((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度温度值的偏差,得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差(对应于上面的3)。

就是这样,卡尔曼滤波器就不断的把covariance递归,从而估算出最优的温度值。他运行的很快,而且它只保留了上一时刻的covariance。上面的Kg,就是卡尔曼增益(Kalman Gain)。他可以随不同的时刻而改变他自己的值,是不是很神奇!

下面就要言归正传,讨论真正工程系统上的卡尔曼。

3
卡尔曼滤波器算法
The Kalman Filter Algorithm

在这一部分,我们就来描述源于Dr Kalman 的卡尔曼滤波器。下面的描述,会涉及一些基本的概念知识,包括概率(Probability),随即变量(Random Variable),高斯或正态分配(Gaussian Distribution)还有State-space Model等等。但对于卡尔曼滤波器的详细证明,这里不能一一描述。

首先,我们先要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程(Linear Stochastic Difference equation)来描述:
X(k)=A X(k-1)+B U(k)+W(k)
再加上系统的测量值:
Z(k)=H X(k)+V(k)
上两式子中,X(k)k时刻的系统状态,U(k)k时刻对系统的控制量。AB是系统参数,对于多模型系统,他们为矩阵。Z(k)k时刻的测量值,H是测量系统的参数,对于多测量系统H为矩阵。W(k)V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(White Gaussian Noise),他们的covariance 分别是QR(这里我们假设他们不随系统状态变化而变化)。

对于满足上面的条件(线性随机微分系统,过程和测量都是高斯白噪声),卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的covariances 来估算系统的最优化输出(类似上一节那个温度的例子)。

首先我们要利用系统的过程模型,来预测下一状态的系统。假设现在的系统状态是k,根据系统的模型,可以基于系统的上一状态而预测出现在状态:
X(k|k-1)=A X(k-1|k-1)+B U(k) ……….. (1)
(1)中,X(k|k-1)是利用上一状态预测的结果,X(k-1|k-1)是上一状态最优的结果,U(k)为现在状态的控制量,如果没有控制量,它可以为0

到现在为止,我们的系统结果已经更新了,可是,对应于X(k|k-1)covariance还没更新。我们用P表示covariance
P(k|k-1)=A P(k-1|k-1) A’+Q ……… (2)
(2)中,P(k|k-1)X(k|k-1)对应的covarianceP(k-1|k-1)X(k-1|k-1)对应的covarianceA’表示A的转置矩阵,Q系统过程的covariance。式子12就是卡尔曼滤波器5个公式当中的前两个,也就是对系统的预测。

现在我们有了现在状态的预测结果,然后我们再收集现在状态的测量值。结合预测值和测量值,我们可以得到现在状态(k)的最优化估算值X(k|k)
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1)) ……… (3)
其中Kg为卡尔曼增益(Kalman Gain)
Kg(k)= P(k|k-1) H’ / (H P(k|k-1) H’ + R) ……… (4)

到现在为止,我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束,我们还要更新k状态下X(k|k)covariance
P(k|k)=
I-Kg(k) HP(k|k-1) ……… (5)
其中I 1的矩阵,对于单模型单测量,I=1。当系统进入k+1状态时,P(k|k)就是式子(2)P(k-1|k-1)。这样,算法就可以自回归的运算下去。

卡尔曼滤波器的原理基本描述了,式子12345就是他的5 个基本公式。根据这5个公式,可以很容易的实现计算机的程序。

下面,我会用程序举一个实际运行的例子。。。
4
简单例子
A Simple Example

这里我们结合第二第三节,举一个非常简单的例子来说明卡尔曼滤波器的工作过程。所举的例子是进一步描述第二节的例子,而且还会配以程序模拟结果。

根据第二节的描述,把房间看成一个系统,然后对这个系统建模。当然,我们见的模型不需要非常地精确。我们所知道的这个房间的温度是跟前一时刻的温度相同的,所以A=1。没有控制量,所以U(k)=0。因此得出:
X(k|k-1)=X(k-1|k-1) ……….. (6)
式子(2)可以改成:
P(k|k-1)=P(k-1|k-1) +Q ……… (7)

因为测量的值是温度计的,跟温度直接对应,所以H=1。式子345可以改成以下:
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1)) ……… (8)
Kg(k)= P(k|k-1) / (P(k|k-1) + R) ……… (9)
P(k|k)=
1-Kg(k)P(k|k-1) ……… (10)

现在我们模拟一组测量值作为输入。假设房间的真实温度为25度,我模拟了200个测量值,这些测量值的平均值为25度,但是加入了标准偏差为几度的高斯白噪声(在图中为蓝线)。

为了令卡尔曼滤波器开始工作,我们需要告诉卡尔曼两个零时刻的初始值,是X(0|0)P(0|0)。他们的值不用太在意,随便给一个就可以了,因为随着卡尔曼的工作,X会逐渐的收敛。但是对于P,一般不要取0,因为这样可能会令卡尔曼完全相信你给定的X(0|0)系统最优的,从而使算法不能收敛。我选了X(0|0)=1度,P(0|0)=10

系统的真实温度为25度,图中用黑线表示。图中红线是卡尔曼滤波器输出的最优化结果(该结果在算法中设置了Q=1e-6R=1e-1)。