基于高次曲線擬合的人臉特征表示

劉 宏，張樹功

（吉林大學(xué)數(shù)學(xué)研究所，長(zhǎng)春 130012）

基于高次曲線擬合的人臉特征表示

劉 宏，張樹功

（吉林大學(xué)數(shù)學(xué)研究所，長(zhǎng)春 130012）

基于數(shù)字化曲線擬合技術(shù)，提出一種新的人臉圖像不變特征表示方法.該方法先將人臉邊緣圖像分割為高次多項(xiàng)式及拋物線段構(gòu)成的特征曲線段集合，再運(yùn)用極慣性矩度量不同人臉圖像在特征曲線段兩側(cè)的局部紋理及形態(tài)差異.實(shí)驗(yàn)結(jié)果表明，該方法與直線邊緣圖法和拋物線邊緣圖法相比，能更好地描述人臉圖像的分布特征.

邊緣提取；曲線擬合；高次多項(xiàng)式；極慣性矩

1 人臉邊緣圖像的曲線表示及相應(yīng)不變特征的提取

1.1 邊緣提取及線段擬合

對(duì)于待識(shí)別的人臉圖像，首先要進(jìn)行圖像邊緣的提取.本文采用Canny邊緣提取算子［9］，該算子較其他算子對(duì)圖像隨機(jī)噪聲更不敏感，提取出的邊緣圖像較平滑和穩(wěn)定.在得到邊緣圖像后，即可利用文獻(xiàn)［10－11］中的線段擬合方法將其分割成若干直線段、拋物線段或高次多項(xiàng)式線段組成的線段圖.該線段擬合方法稱為數(shù)字曲線結(jié)構(gòu)性擬合，其思想是利用數(shù)字曲線各點(diǎn)間的拓?fù)潢P(guān)系限制，對(duì)待擬合曲線中選取的若干點(diǎn)集進(jìn)行擬合誤差估計(jì)，從而在避免對(duì)待擬合曲線全體點(diǎn)集進(jìn)行擬合誤差估計(jì)的前提下，求出擬合誤差小于預(yù)設(shè)值的多項(xiàng)式擬合曲線.該方法可簡(jiǎn)要表述如下：

對(duì)由n個(gè)點(diǎn)組成的數(shù)字曲線

令G為l次多項(xiàng)式擬合函數(shù)構(gòu)成的向量空間，l?n.其中任意函數(shù)g（x）可表示為

記g（x）對(duì)曲線C的擬合誤差為

則對(duì)曲線C的最優(yōu)擬合函數(shù)誤差可表示為

此時(shí)，擬合多項(xiàng)式為l次多項(xiàng)式，則待擬合曲線C中所有點(diǎn)數(shù)為l＋2的子集均可用來(lái)估計(jì)最優(yōu)擬合多項(xiàng)式及其擬合誤差［10］.這些子集稱為曲線C的基礎(chǔ)子集，記這些基礎(chǔ)子集的集合為M.

對(duì)某基礎(chǔ)子集D∈M，

即使用基礎(chǔ)子集可對(duì)整個(gè)曲線的擬合誤差進(jìn)行估計(jì)［10］.特別地，當(dāng)基礎(chǔ)子集M具有一定拓?fù)潢P(guān)系限制，如坐標(biāo)間距為固定值且點(diǎn)與點(diǎn)之間關(guān)系為四連通或八連通時(shí)，用基礎(chǔ)子集M的剛性子集［11］的擬合誤差~r（c）即可對(duì)r（C）進(jìn)行估計(jì).顯然，數(shù)字圖像中構(gòu)成曲線的點(diǎn)集即滿足該條件.

剛性子集的選取有多種選擇.以直線擬合為例，對(duì)有i＋1個(gè)點(diǎn)的點(diǎn)集C＝（p1，p2，…，pi，pi＋1），其中的任意三點(diǎn)均構(gòu)成該點(diǎn)集的基礎(chǔ)子集，而如（p1，pi／2＋1，pi＋1），（p1，pi，pi＋1）等則可構(gòu)成其剛性子集［11］.進(jìn)行誤差估計(jì)時(shí)，既可使用其中某一剛性子集，也可同時(shí)使用多個(gè)剛性子集.當(dāng)使用多個(gè)子集時(shí)，其中之一的擬合誤差超過(guò)準(zhǔn)許上限后即結(jié)束擬合.圖1為若干點(diǎn)構(gòu)成的一條數(shù)字曲線，如對(duì)其進(jìn)行直線段擬合，則可選取剛性子集（p1，pi，pi＋1）進(jìn)行擬合誤差估計(jì).

利用該方法，可將邊緣圖像分割為若干直線段或多項(xiàng)式曲線段.圖2（A）為人臉邊緣圖像；圖2（B）為對(duì)邊緣圖像進(jìn)行直線段擬合后的結(jié)果；圖2（C）為拋物線段擬合的結(jié)果.

比較邊緣圖像和擬合結(jié)果可見，由于面部輪廓線多呈現(xiàn)出曲線特性，僅使用直線段進(jìn)行擬合，易導(dǎo)致原本連續(xù)的曲線特征被人為截?cái)?，破壞了特征的整體性，弱化了特征的判別力.在該問(wèn)題上，拋物線擬合比直線段擬合的效果更好.但僅使用拋物線進(jìn)行擬合，仍然有明顯的局限性.圖3（A）框中為人臉鼻翼側(cè)自然褶皺形成的連貫邊緣線段，圖3（B）和圖3（C）分別為拋物線段擬合和高次多項(xiàng)式線段擬合結(jié)果.由圖3可見，高次多項(xiàng)式線段更真實(shí)地還原了原始圖像的分布特征.

對(duì)人臉圖像進(jìn)行邊緣提取后，部分區(qū)域的邊緣線曲線特征不明顯，長(zhǎng)度也較短，使用直線段或拋物線段擬合，可使特征提取簡(jiǎn)化，同時(shí)又不會(huì)產(chǎn)生特征的破碎化.而對(duì)如圖3中所示的若干區(qū)域，使用高次多項(xiàng)式曲線則能最大限度地還原原始圖像的分布特性，提高擬合精度.在實(shí)驗(yàn)中，利用面部器官分布的先驗(yàn)知識(shí)，集中對(duì)此類區(qū)域中（如鼻翼、臉頰、唇邊）的邊緣進(jìn)行高次多項(xiàng)式線段擬合，比僅采用單一的低次曲線類型更易取得較好結(jié)果.

圖1 剛性子集Fig.1 Rigid subcollection

圖2 面部邊緣圖像的不同曲線擬合結(jié)果Fig.2 Different curve fitting result for face edge images

圖3 高次曲線擬合Fig.3 High－order curve fitting

1.2 不變特征的選取

如圖4所示，對(duì)人臉圖像（A）進(jìn)行邊緣提取及線段擬合后，可得到線段圖像（B）.對(duì)圖像（B）中箭頭處所指線段，取原始圖像中該線段所在位置兩側(cè)一定寬度內(nèi)的像素點(diǎn)，構(gòu)成該線段的相關(guān)區(qū)域，如圖4（C）.然后將該相關(guān)區(qū)域均分成若干子塊，以每子塊中心位置像素值作為參考，與該子塊內(nèi)其他像素進(jìn)行比較，將與其灰度值相近的像素提取出來(lái)，即可有效刻畫出各線段附近像素的分布情況.圖5（A）為相關(guān)區(qū)域中一子塊；圖5（B）中C，S，D點(diǎn)分別表示子塊的中心點(diǎn)、與中心點(diǎn)相似的點(diǎn)及與中心點(diǎn)相異的點(diǎn)；圖5（C）為經(jīng)過(guò)分類后，子塊所體現(xiàn)出的分布特征.

圖4 原始圖像中抽取的特征Fig.4 Features extracted from original face image

稱相似點(diǎn)構(gòu)成的區(qū)域?yàn)橄嗨茀^(qū)域S，相異點(diǎn)組成的區(qū)域?yàn)橄喈悈^(qū)域D.圖6給出了子塊內(nèi)若干種可能出現(xiàn)的區(qū)域分布情況，其中灰色區(qū)域表示子塊內(nèi)的相似區(qū)域，白色區(qū)域表示相異區(qū)域.此時(shí)可將各子塊視為二維平面上非均勻分布的具有不規(guī)則外形的幾何圖形，對(duì)其特性的量化描述可使用極慣性矩的概念.

圖5 相關(guān)區(qū)域分成的兩類子塊Fig.5 Two types of blocks divided from correlative area

圖6 子塊的不同分布Fig.6 Distribution of the blocks

如圖7所示，A為某平面圖形面積（對(duì)應(yīng)上述子塊的相似區(qū)域），C為該圖形形心.對(duì)于C，可給出極慣性矩J的一般性定義易證，J為不受圖像坐標(biāo)平移、旋轉(zhuǎn)影響的不變量［12］.利用這一特性，可有效刻畫圖像的局部分布特性.由圖6可見，通過(guò)對(duì)子塊內(nèi)不同區(qū)域的劃分，實(shí)際上體現(xiàn)出了人臉圖像局部紋理的固有分布特性.顯然，在面部姿態(tài)和表情未發(fā)生顯著變化的情況下，這些分布特性不會(huì)出現(xiàn)大的改變.

在LEM方法中，不變特征的提取和比對(duì)僅考慮了直線段圖像間的差異，而未考慮原圖像的像素分布差異，無(wú)法有效利用圖像紋理中蘊(yùn)含的判別信息.而PEM方法雖然在識(shí)別過(guò)程中，加入了計(jì)算線段兩側(cè)像素點(diǎn)差異的內(nèi)容，但僅沿線段提取了若干孤立的像素點(diǎn)，未能充分體現(xiàn)線段兩側(cè)相關(guān)區(qū)域的整體分布特征.當(dāng)局部光照條件發(fā)生變化或面部表情發(fā)生小幅度變化時(shí)，即有可能導(dǎo)致采樣點(diǎn)發(fā)生較大改變，從而影響識(shí)別的準(zhǔn)確性.相較于這兩種方法，本文選取的不變特征則能較好地避免上述問(wèn)題.

圖7 極慣性矩Fig.7 Polar moment of inertia

2 算法的實(shí)現(xiàn)

在提取邊緣前，須對(duì)人臉圖像進(jìn)行歸一化處理，根據(jù)CSU人臉識(shí)別評(píng)估體系中提出的人臉圖像歸一化方法［13］，對(duì)來(lái)自于FERET［14］人臉圖片庫(kù)的人臉圖像進(jìn)行以下操作：

1）幾何歸一化.利用眼球中心坐標(biāo)將人臉圖像旋轉(zhuǎn)為水平狀態(tài)，并使圖像具有相同外尺寸；

2）橢圓輪廓化.使用一橢圓面具對(duì)人臉圖像外邊緣進(jìn)行剪裁，僅使額頭、下巴和兩頰以內(nèi)部分可見；

3）直方圖均衡化.對(duì)橢圓面具內(nèi)像素進(jìn)行直方圖均衡化處理.

處理結(jié)果如圖8所示，圖8（A）為原始圖像，圖8（B）為處理后的圖像.

然后利用Canny算子，對(duì)圖像進(jìn)行邊緣提取.在提取的邊緣圖像上，利用曲線擬合方法，完成線段分割.提取線段后，即可對(duì)原圖像中該線段對(duì)應(yīng)的相關(guān)區(qū)域進(jìn)行分析.如圖9所示，對(duì)于提取出的線段AB（圖9（A），（B），（C）中，AB分別為直線段、拋物線段和高次多項(xiàng)式線段），沿該線段兩端點(diǎn)連線的垂線方向，將AB向其兩側(cè)分別平移若干像素，得到新線段A′B′和A″B″.這兩條新線段和直線段A′A″及B′B″構(gòu)成了一個(gè)封閉圖形.該封閉圖形在原圖像中圍出的區(qū)域即為線段AB在原圖像中的相關(guān)區(qū)域.等距平移線段A′B′至A″B″，將相關(guān)區(qū)域等分為m份，再等距平移A′A″至B′B″，將區(qū)域均分為n份，即可將相關(guān)區(qū)域劃分為m×n個(gè)子塊（以圖9為例，m＝n＝6，則區(qū)域被劃分為36個(gè)子塊）.此時(shí)，當(dāng)子塊長(zhǎng)寬均為奇數(shù)個(gè)像素時(shí)，可直接得到中心像素的坐標(biāo)，并提取其灰度值；若子塊一邊為偶數(shù)，則取位于中心位置兩點(diǎn)灰度值的均值；若子塊長(zhǎng)寬均為偶數(shù)，則取位于中心位置的四點(diǎn)灰度值均值.為了克服隨機(jī)噪聲對(duì)實(shí)驗(yàn)的影響，可用中心點(diǎn)及其周圍鄰域中值濾波后的灰度值代替原灰度值.利用該灰度值，考察相關(guān)子塊中像素分布情況.再利用極慣性矩，求各子塊關(guān)于相關(guān)區(qū)域中心點(diǎn)的不變特征.用極慣性矩平行移軸公式［12］，可得該像素關(guān)于某點(diǎn)（x，y）的極慣性矩

圖8 均一化人臉圖像Fig.8 Unified face images

圖9 相關(guān)區(qū)域網(wǎng)格Fig.9 Correlative area grid

圖10 極慣性矩的計(jì)算Fig.10 Calculation of polar moment of inertia

即先求得子塊內(nèi)某像素關(guān)于自身中心的極慣性矩，再求出其關(guān)于整個(gè)相關(guān)區(qū)域中心的極慣性矩.

進(jìn)一步，考慮到像素灰度值可視為人臉在一定光照條件下的反射響應(yīng)，可給出能量函數(shù)的定義

其中：Eij為i行j列像素的能量函數(shù)；αij為與人臉物理特性相關(guān)的權(quán)重系數(shù)；Gij為i行j列像素的灰度值.對(duì)某相關(guān)區(qū)域內(nèi)第k個(gè)子塊內(nèi)所有相似點(diǎn)的能量函數(shù)求和，記為Ek.因此，以各子塊能量函數(shù)組成該相關(guān)區(qū)域的特征序列，即每個(gè)相關(guān)區(qū)域的特征序列均由m×n個(gè)特征構(gòu)成（以圖9為例，該相關(guān)區(qū)域的特征為36個(gè)），而每個(gè)相關(guān)區(qū)域的特征序列又可視為整個(gè)圖像的一個(gè)特征.

分別獲得待測(cè)圖像和樣本圖像的特征序列后，對(duì)待測(cè)圖像中某線段，可在樣本圖像中尋找與其坐標(biāo)最接近的線段作為對(duì)應(yīng)線段，比較二者間相關(guān)區(qū)域不變特征序列間的差異.考慮到這兩個(gè)相關(guān)區(qū)域未必完全重疊，且顯然兩區(qū)域重疊度越高，二者間的差異越具有判別性，故在計(jì)算最終差異時(shí)，需加入重疊度的計(jì)算.圖11（A）和圖11（B）分別表示兩個(gè)對(duì)應(yīng)線段的相關(guān)區(qū)域，圖11（C）中陰影區(qū)域即二者的重疊部分.

圖11 重疊區(qū)域Fig.11 Overlap area

考慮兩個(gè)對(duì)應(yīng)線段（即坐標(biāo)最相近的兩線段）的相關(guān)區(qū)域，兩者中較大區(qū)域的面積記為Amax，重疊區(qū)域面積記為Aol，重疊率R＝Aol／Amax.設(shè)參加比對(duì)的兩個(gè)樣本分別有M1，M2個(gè)相關(guān)區(qū)域，M1＜M2，每個(gè)相關(guān)區(qū)域有K個(gè)子塊，記某樣本第j個(gè)相關(guān)區(qū)域內(nèi)第i個(gè)子塊的相似區(qū)域Sij對(duì)應(yīng)的能量函數(shù)為eij，則第j個(gè)相關(guān)區(qū)域的能量函數(shù)為

若分別記兩樣本對(duì)應(yīng)的某相關(guān)區(qū)域的能量函數(shù)為Eaj和Ebj，則樣本間差異為

3 實(shí) 驗(yàn)

選取圖片庫(kù)中150人作為實(shí)驗(yàn)對(duì)象，對(duì)其多幅圖像提取邊緣并進(jìn)行高次多項(xiàng)式擬合及不變特征提取.通過(guò)對(duì)邊緣提取及擬合階段各閾值的調(diào)整，可控制每幅圖像相關(guān)區(qū)域的數(shù)量.此時(shí)，可視每個(gè)相關(guān)區(qū)域?yàn)閳D像的一個(gè)特征，該特征由其內(nèi)部特征序列描述.隨后，取每個(gè)實(shí)驗(yàn)對(duì)象若干圖像作為該實(shí)驗(yàn)對(duì)象的標(biāo)記樣本組，其余圖像則為待測(cè)樣本.當(dāng)待測(cè)樣本與某標(biāo)記樣本組中的一幅圖像特征差異最小，且兩者確屬同一人的圖像時(shí)，認(rèn)為識(shí)別成功.如圖12（A）所示，A，B，C三組數(shù)據(jù)分別表示使用3，2，1個(gè)圖像作為標(biāo)記樣本時(shí)，在不同特征數(shù)量（即相關(guān)區(qū)域的數(shù)量）下，算法識(shí)別率的變化情況.如圖12（B）所示，A，B，C三組數(shù)據(jù)分別表示使用直線段、高次多項(xiàng)式與拋物線混合線段、拋物線段提取特征后，在不同特征數(shù)量下算法的效果.

圖12 不同實(shí)驗(yàn)方法的結(jié)果Fig.12 Results by different experiment methods

由圖12可見，當(dāng)其他條件相同時(shí)，使用同一人的多張圖像構(gòu)成樣本能更充分體現(xiàn)其在不同光照、表情等情況下的特性，從而取得更好的識(shí)別結(jié)果.另一方面，當(dāng)樣本特征數(shù)量由5開始逐步增加時(shí)，識(shí)別率有明顯提高，但數(shù)量增加到一定程度后，識(shí)別率反而有所下降.經(jīng)分析認(rèn)為，當(dāng)線段及其相關(guān)區(qū)域達(dá)到一定數(shù)量后，即能有效刻畫人臉特征的分布情況，如繼續(xù)增加特征數(shù)量，相關(guān)區(qū)域間的重疊情況會(huì)明顯增加，過(guò)多的冗余信息反而弱化了最具判別力的若干單個(gè)特征在特征比對(duì)中的作用，從而影響了算法的效果.因此，特征數(shù)量的控制及特征的選擇能在很大程度上控制識(shí)別結(jié)果的優(yōu)劣.本文目前主要從圖像邊緣的強(qiáng)弱和線段的長(zhǎng)度等方面對(duì)特征進(jìn)行篩選，當(dāng)圖像小范圍內(nèi)出現(xiàn)多條長(zhǎng)度和強(qiáng)度都相近的邊緣線段時(shí)，這些線段都將參與最終的特征比對(duì)，而最有效的方法則是在其中挑選出最具代表性的線段，從而避免信息冗余對(duì)判別帶來(lái)的干擾.

綜上所述，本文采用高次多項(xiàng)式曲線擬合算法分割圖像邊緣，在獲得圖像較精確描述的同時(shí)，與直線段相比，提高了算法復(fù)雜度.實(shí)驗(yàn)結(jié)果表明，在相同條件下，分別使用兩種擬合方法對(duì)同一圖片進(jìn)行擬合運(yùn)算時(shí)，高次多項(xiàng)式曲線擬合與直線段擬合的時(shí)間開銷相比雖有一定程度增加，但仍屬10ms量級(jí).因此，其計(jì)算成本是可接受的.

［1］ Patil A M，Kolhe S R.2DFace Recognition Techniques：A Survey［J］.International Journal of Machine Intelligence，2010，2：74－83.

［2］ Biederman I，Ju G.Surface versus Edge－Based Determinants of Visual Recognition［J］.Cognitive Psychology，1988，20（1）：38－64.

［3］ Bruce V，Hanna E，Dench N，et al.The Importance of“Mass”in Line Drawings of Faces［J］.Applied Cognitive Psychology，1992，6：619－628.

［4］ Takács B.Comparing Face Images Using the Modified Hausdorff Distance［J］.Pattern Recognition，1998，31：1873－1881.

［5］ ZHAO Wenyi.Face Processing：Advanced Modeling and Methods［M］.Burlington：Elsevier，2006：254－286.

［6］ GAO Yongsheng，Leung M K H.Face Recognition Using Line Edge Map［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2002，24（6）：764－779.

［7］ DU Cheng，SU Guangda，LIN Xinggang.Face Recognition Using a Modified Line Segment Hausdorff Distance［C］／／Proceedings of the Second International Conference on Machine Learning and Cybernetics.Xi’an：IEEE，2003：3016－3020.

［8］ Deboeverie F，Veelaert P.Face Recognition Using Parabola Edge Map［J］.Advanced Concepts for Intelligent Vision Systems Lecture Notes in Computer Science，2008，5259：994－1005.

［9］ Canny J.A Computational Approach to Edge Detection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（6）：679－685.

［10］ Veelaert P，Teelen K.Fast Polynomial Segmentation of Digitized Curves［J］.Discrete Geometry for Computer Imagery，2006，4245：482－493.

［11］ Veelaert P.Constructive Fitting and Extraction of Geometric Primitives［J］.Graphical Models and Image Processing，1997，59（4）：233－245.

［12］ 范欽珊.工程力學(xué)教程［M］.北京：高等教育出版社，1998.（FAN Qinshan.Engineering Mechanics［M］.Beijing：Higher Education Press，1998.）

［13］ Bolme D S，Beveridge J R.The CSU Face Identification Evaluation System：Its Purpose，F(xiàn)eatures and Structure［J］.Computer Vision Systems，2003，2626：304－313.

［14］ Phillips P J，Moon H J.The FERET Evaluation Methodology Face－Recognition Algorithms［EB／OL］.2011－02－15.http：／／www.nist.gov／itl／iad／ig／colorferet.cfm／.

（責(zé)任編輯：韓 嘯）

Face Eigen Expression Based on High－Order Curve Fitting

LIU Hong，ZHANG Shugong
（Institute of Mathematics，Jilin University，Changchun130012，China）

We proposed a new method to extract invariant features from human face image.With the help of digital curve fitting，we represented face binary edge map with high－order polynomial and parabola segments.Then，the polar moment of the inertia of pixel block beside the curve segments was used to measure the difference between local texture and shape.From the analysis of the experimental result，it is concluded that compared with line edge map and parabola edge map，this method can give a better way to describe the feature distribution of face images.

edge extraction；curve fitting；high－order polynomial；polar moment of inertia

TP391.4

A

1671－5489（2014）04－0746－07

在人臉識(shí)別領(lǐng)域，相對(duì)于三維圖像和紅外圖像等數(shù)據(jù)來(lái)源，二維圖像具有更便于獲取和處理的優(yōu)勢(shì)，因此，基于二維圖像的人臉識(shí)別技術(shù)得到了廣泛應(yīng)用［1］.其中特征臉?lè)?、彈性模板匹配法和Bayes法等分別利用了二維圖像中蘊(yùn)含的矩陣特征、小波分布特征和概率分布特征.研究表明［2－5］，面部邊緣圖像中已蘊(yùn)含了豐富的、可供辨識(shí)的信息.因此，除上述幾種方法外，還有一種算法，試圖通過(guò)提取邊緣圖像中相應(yīng)的幾何分布特征完成對(duì)人臉的分類或識(shí)別.如Gao等［6］提出的LEM（line edge map，即直線邊緣圖）方法，利用人臉圖像的邊緣圖像，生成相應(yīng)的直線段圖集，并由此提取其中具有辨識(shí)力的特征［7］.Deboeverie等［8］提出了PEM（parabola edge map，即拋物線邊緣圖）方法，運(yùn)用拋物曲線完成人臉邊緣圖像的分割和特征提取，使這種算法有了進(jìn)一步的發(fā)展.以上兩種方法均使用單一類型的線段提取圖像特征，雖然便于數(shù)據(jù)處理和存儲(chǔ)，但在處理某些相對(duì)復(fù)雜的面部曲線時(shí)，會(huì)不可避免地將原來(lái)具有整體性的特征碎片化，從而使算法的誤差增大、判別能力降低.本文針對(duì)該問(wèn)題，充分利用人臉器官輪廓線分布規(guī)律的先驗(yàn)知識(shí)，引入高次曲線對(duì)邊緣圖像進(jìn)行分割，使提取出的線段特征更具整體性和穩(wěn)定性，并在此基礎(chǔ)上，引入新的不變特征提取方法，取得了較好的實(shí)驗(yàn)效果.

10.13413／j.cnki.jdxblxb.2014.04.22

2013－10－18.

劉 宏（1982—），男，漢族，博士研究生，從事數(shù)值計(jì)算和模式識(shí)別的研究，E－mail：hongliu07＠m(xù)ails.jlu.edu.cn.通信作者：張樹功（1958—），男，漢族，博士，教授，博士生導(dǎo)師，從事數(shù)值計(jì)算和計(jì)算機(jī)代數(shù)的研究，E－mail：sgzh＠m(xù)ail.jlu.edu.cn.

國(guó)家自然科學(xué)基金（批準(zhǔn)號(hào)：11171003）和吉林省自然科學(xué)基金（批準(zhǔn)號(hào)：20101597）.