基于子孔徑參數(shù)估計(jì)的雙基地ISAR圖像融合方法研究

許 然 李亞超 邢孟道

(西安電子科技大學(xué)雷達(dá)信號(hào)處理國(guó)家重點(diǎn)實(shí)驗(yàn)室 西安 710071)

1 引言

逆合成孔徑雷達(dá)(ISAR)具有全天候、全天時(shí)和遠(yuǎn)距離觀察等特性，可以對(duì)目標(biāo)進(jìn)行高分辨率成像和識(shí)別[1,2]，獲取相關(guān)信息，在現(xiàn)代戰(zhàn)爭(zhēng)中發(fā)揮著關(guān)鍵的作用。它通過(guò)發(fā)射大帶寬信號(hào)獲得距離向高分辨，對(duì)來(lái)自不同角度的回波進(jìn)行相干處理獲得方位向高分辨。為了提高對(duì)目標(biāo)的成像和識(shí)別能力，常常需要較長(zhǎng)的相干處理時(shí)間(或更大的相干積累角)，但這樣會(huì)引入非理想因素，使信號(hào)模型復(fù)雜化，給處理帶來(lái)較大難度。

與單基地雷達(dá)相比，雙基地雷達(dá)具有更遠(yuǎn)的作用距離、更強(qiáng)的抗干擾和獲取信息能力等特點(diǎn)[3-5]。文獻(xiàn)[5]中對(duì)利用水面反射回波進(jìn)行雙基ISAR成像的可行性進(jìn)行了研究；張亞波等人[6]對(duì)雙基地ISAR成像算法理論進(jìn)行了詳細(xì)討論；文獻(xiàn)[7]分析了雙基地 ISAR數(shù)據(jù)的波數(shù)譜。通過(guò)多個(gè)雷達(dá)系統(tǒng)，將在不同視角獲取的目標(biāo)圖像進(jìn)行融合處理，可以獲取質(zhì)量更高的圖像以及更加豐富的目標(biāo)信息[8,9]。雖然利用已經(jīng)存在的很多圖像融合算法[10,11]可以對(duì)兩幅圖像進(jìn)行融合，但這些算法大多數(shù)都需要在圖像域先對(duì)兩幅圖像進(jìn)行伸縮旋轉(zhuǎn)完成配準(zhǔn)，再使用小波變換或其他方法進(jìn)行融合，并沒(méi)有運(yùn)用到任何ISAR圖像特性。文獻(xiàn)[9]則最先提出了一種基于MFT(Matrix Fourier Transform)的多視角ISAR圖像融合算法，填補(bǔ)了這項(xiàng)空白，但其采用的并不是一發(fā)雙收的雙基地ISAR雷達(dá)體制。如果對(duì)雙基地雷達(dá)的收發(fā)分置效應(yīng)加以利用，便可以得到不同視角下的目標(biāo)回波，對(duì)其進(jìn)行運(yùn)動(dòng)補(bǔ)償后分別成像，再將兩幅圖像相干融合疊加，可以減輕噪聲的影響，改善信噪比與圖像質(zhì)量，得到更精確的目標(biāo)特征信息；而且在散射各向異性和存在遮擋的情況下，將不同視角圖像融合相加后可以得到更為真實(shí)全面的目標(biāo)形狀。但同時(shí)也正是這種分置效應(yīng)，使得兩幅目標(biāo)像的形狀和尺寸各異，無(wú)法直接在圖像域疊加，造成融合困難。

本文采用一發(fā)雙收的雙基地雷達(dá)配置，首先建立回波信號(hào)模型并進(jìn)行成像及融合原理分析，然后利用孤立強(qiáng)散射點(diǎn)在不同子孔徑間的多普勒差異估計(jì)出目標(biāo)轉(zhuǎn)角速度和半雙基地角，在數(shù)據(jù)域完成對(duì)兩個(gè)雷達(dá)各自獲取的兩幅圖像的融合。最后，通過(guò)仿真驗(yàn)證了該方法的可行性和有效性。

2 雙基地ISAR雷達(dá)模型

假設(shè)雷達(dá)與目標(biāo)之間的平動(dòng)分量可以被精確補(bǔ)償，令逆時(shí)針旋轉(zhuǎn)方向?yàn)檎较?，建立如圖1所示的雙基地ISAR雷達(dá)模型。雷達(dá)1發(fā)射信號(hào)，雷達(dá)1和雷達(dá)2同時(shí)接收信號(hào)。圖中RT,RR分別是雷達(dá)1和雷達(dá)2到目標(biāo)旋轉(zhuǎn)參考點(diǎn)O的距離矢量，其模值分別為RT,RR，單位矢量分別為和且與地面夾角分別為θT和θR，收發(fā)視線的和向量REq=RT+RR，也即為角平分線向量，為半雙基角，為雙基地雷達(dá)等效視線方向角。

以雷達(dá)1視線方向和等效雷達(dá)視線方向分別建立直角坐標(biāo)系XOY和UOV。設(shè)P為目標(biāo)物體上任意散射點(diǎn)，rp為其位置矢量，且與U軸夾角大小為θ。則該散射點(diǎn)到雷達(dá)1和雷達(dá)2的距離分別為

圖1 雙基地ISAR雷達(dá)轉(zhuǎn)臺(tái)模型

式中tm是方位時(shí)間。由于一般目標(biāo)物體尺寸遠(yuǎn)小于到雷達(dá)的距離，即rp?RR,RT，故式(1)，式(2)可近似成：

其中up=rpcosθ,vp=rpsinθ，表示P點(diǎn)在直角坐標(biāo)系UOV中的坐標(biāo)。

3 雙基地ISAR雷達(dá)成像分析

假設(shè)雷達(dá)1發(fā)射周期性線性調(diào)頻信號(hào)為

式中t=+tm為全時(shí)間，為快時(shí)間，tm=mTr(m=0,1,2,…)為方位慢時(shí)間，Tr表示方位脈沖時(shí)間間隔，Tp為脈沖寬度，γ=B/Tp為調(diào)頻率，B為信號(hào)帶寬，σc為信號(hào)復(fù)振幅，則雷達(dá)k(k=1,2)接收到來(lái)自P點(diǎn)的回波基帶信號(hào)可表示成

式中Rk(tm)為雷達(dá)k與P點(diǎn)之間的斜距，將P點(diǎn)在XOY中的坐標(biāo)xp=rpcos(θ+β),yp=rpsin (θ+β)代入式(1)可以得到

同時(shí)，R2(tm)=(Rt(tm)+Rr(tm))/2表示雷達(dá)2接收到P點(diǎn)回波的等效斜距，代入式(3)，式(4)化簡(jiǎn)得

在對(duì)兩個(gè)雷達(dá)采集的回波基帶信號(hào)進(jìn)行距離向匹配濾波后得到

在式(7)，式(8)中，RT和(RT+RR)/2分別表示目標(biāo)相對(duì)于雷達(dá)視線的平動(dòng)分量，對(duì)其補(bǔ)償后僅剩下轉(zhuǎn)動(dòng)分量，此時(shí)兩個(gè)回波信號(hào)在忽略常數(shù)相位項(xiàng)后在距離頻域-方位時(shí)域可表示為

需要指出的是，一般實(shí)際情況中RT,RR均是隨方位時(shí)間在改變的，因此雙基地雷達(dá)成像的距離軸是時(shí)變的，并且與方位軸不正交[12]。但如若選取較短的成像時(shí)間，可認(rèn)為RT,RR變化較小，雙基地角的時(shí)變效應(yīng)可以忽略，即認(rèn)為整個(gè)成像時(shí)間段中β近似常數(shù)，則將平動(dòng)分量補(bǔ)償?shù)艉蟊憧蛇\(yùn)用以上轉(zhuǎn)臺(tái)模型進(jìn)行成像。而且由于相干積累角較小，其轉(zhuǎn)動(dòng)角速度ω也可近似認(rèn)為是常量，并不需要考慮其高階項(xiàng)，因此大大簡(jiǎn)化了雙基地ISAR雷達(dá)模型且縮短了數(shù)據(jù)量，方便于快速成像處理。但是由于選取了較短的成像時(shí)間，造成分辨率較低，最后圖像質(zhì)量不夠高；如果將兩個(gè)雷達(dá)單獨(dú)錄取的回波數(shù)據(jù)進(jìn)行相干融合處理，便可達(dá)到增大方位積累角的目的，最終提高成像質(zhì)量，并獲取更多視角下目標(biāo)的散射特性，大大增加對(duì)目標(biāo)的特征提取和識(shí)別能力。

4 雙基地ISAR的圖像融合

4.1 圖像融合原理

根據(jù)式(10)，結(jié)合以上分析，可以得到雷達(dá)1的目標(biāo)基頻總回波為

將式(13)代入式(12)化簡(jiǎn)后得到

從圖2中可以看出，半雙基地角會(huì)直接導(dǎo)致雷達(dá)2回波數(shù)據(jù)的波數(shù)譜支撐區(qū)域大小在方位向和距離向上均多出一個(gè)cosβ因子而變小，使得雙基ISAR圖像的距離向和方位向分辨率分別退化為這可以通過(guò)將雷達(dá)2數(shù)據(jù)的波數(shù)譜進(jìn)行插值加以解決。而直接運(yùn)用RD算法得到兩個(gè)雷達(dá)的兩幅圖像，由于分辨率以及雷達(dá)觀測(cè)視線的不一致將會(huì)導(dǎo)致圖像中目標(biāo)的大小與角度不同，這些由半雙基地角β帶來(lái)的問(wèn)題必須通過(guò)分析雙基地ISAR系統(tǒng)特性，在數(shù)據(jù)域就加以解決，以得到更好的ISAR圖像融合結(jié)果。實(shí)際處理中，考慮到噪聲相位的隨機(jī)性和非相干性，以及雷達(dá)系統(tǒng)帶寬以及方位積累角的限制，如果把兩個(gè)雷達(dá)獲取的數(shù)據(jù)利用BP(Back Projection)算法對(duì)回波相位完成補(bǔ)償后進(jìn)行相干融合，那么噪聲便可以相互抵消，使得成像質(zhì)量得到改善，處理過(guò)程可由式(15)表示為

式(15)中 Map[·]表示將雷達(dá)2數(shù)據(jù)從UOV坐標(biāo)系中映射到XOY坐標(biāo)系的變換函數(shù)，(x,y)為重構(gòu)得到的目標(biāo)散射函數(shù)。該融合處理需要β和ω兩個(gè)參數(shù)，下面將介紹一種適用于本文圖像融合的β和ω的估計(jì)方法。該方法首先將同一幅圖像分為兩個(gè)子孔徑，通過(guò)提取強(qiáng)孤立散射點(diǎn)在兩個(gè)子孔徑下的多普勒頻率差異，估計(jì)得到ω并完成對(duì)兩個(gè)雷達(dá)所獲目標(biāo)圖像的方位向重新定標(biāo)，再利用該散射點(diǎn)在兩幅圖像中的坐標(biāo)計(jì)算出β。

圖2 雷達(dá)1與雷達(dá)2接收回波的波數(shù)譜區(qū)域

4.2 基于子孔徑的參數(shù)估計(jì)

對(duì)經(jīng)過(guò)距離壓縮和運(yùn)動(dòng)補(bǔ)償后的ISAR回波數(shù)據(jù)進(jìn)行方位 FFT便可得到目標(biāo)的 RD(Range Doppler)像，縱坐標(biāo)反映了目標(biāo)上各散射點(diǎn)的距離向位置，橫坐標(biāo)則代表了各散射點(diǎn)的多普勒頻率。國(guó)內(nèi)外很多學(xué)者提出了不同的方法進(jìn)行轉(zhuǎn)動(dòng)角速度的估計(jì)，完成對(duì)ISAR圖像的方位向重新定標(biāo)[13,14]。

將目標(biāo)的回波數(shù)據(jù)在時(shí)域分為前后兩個(gè)子孔徑然后運(yùn)用RD算法，可以得到兩幅ISAR圖像。因?yàn)槌上駮r(shí)間段較短，這兩幅圖像具有極強(qiáng)的相似性，但子孔徑間的時(shí)間差會(huì)帶來(lái)散射點(diǎn)多普勒頻率的些許偏移，且該偏移是由目標(biāo)轉(zhuǎn)動(dòng)引起的，那么利用其多普勒頻率在兩個(gè)子孔徑間的差異就可以對(duì)目標(biāo)轉(zhuǎn)動(dòng)角速度進(jìn)行估計(jì)。圖3表明了散射點(diǎn)在前后子孔徑下的多普勒頻率差異。

圖3 散射點(diǎn)前后子孔徑下多普勒頻率差異示意圖

上文中式(7)給出了目標(biāo)上任意散射點(diǎn)P與雷達(dá)1之間的斜距隨方位慢時(shí)間的變化關(guān)系，其中的轉(zhuǎn)動(dòng)分量正是方位多普勒頻率的來(lái)源：

式(16)中，由于整個(gè)成像過(guò)程轉(zhuǎn)動(dòng)角較小，sinωtm≈ωtm,cosωtm≈ 1 。將整個(gè)成像時(shí)間段分為前后兩個(gè)子孔徑，令子孔徑間的時(shí)間差為Δtm,Δtm=N/ 2Prf ，并且有sinωΔtm≈ωΔtm,cosωΔtm≈1，則R1Ω在后孔徑下的表達(dá)式為

分別得到該散射點(diǎn)在兩個(gè)子孔徑中的多普勒頻率：

那么該散射點(diǎn)在兩個(gè)子孔徑下的多普勒頻率差值為

式(19)右端各參數(shù)中， Δtm,λ為已知系統(tǒng)參量，而散射點(diǎn)的縱向位置yp則已經(jīng)在RD平面圖像中獲得，那么通過(guò)提取該散射點(diǎn)兩個(gè)孔徑間的Δfd便可以利用式(20)對(duì)ω進(jìn)行估計(jì)：

由于在式(16)，式(17)中使用了近似，而實(shí)際上：

其中對(duì)xpsin (ωΔtm)c os (ωtm)一項(xiàng)使用近似后在式(19)中無(wú)法體現(xiàn)出其影響，雖然其值很小但仍會(huì)帶來(lái)測(cè)量誤差。故在使用該方法時(shí)應(yīng)盡量選取零多普勒頻率附近的散射點(diǎn)使xpsin (ωΔtm)c os (ωtm)盡量小以減小其影響。同時(shí)為了提高ω的估計(jì)精度，可以對(duì)該散射點(diǎn)所在距離單元的數(shù)據(jù)進(jìn)行插值，插值倍數(shù)越高，獲取的多普勒頻率信息越精確；而yp的測(cè)量精度依賴(lài)于系統(tǒng)帶寬，帶寬越大，yp值越準(zhǔn)確，ω的估計(jì)值也越接近真實(shí)值。

在估計(jì)得到ω后，就可以利用fd=2xpω/λ對(duì)兩個(gè)雷達(dá)獲取的圖像進(jìn)行方位向重新定標(biāo)，得到目標(biāo)上各散射點(diǎn)分別在兩個(gè)坐標(biāo)系XOY和UOV中的真實(shí)尺寸與形狀信息。需要注意的是，由于兩幅圖像的分辨率相差了一個(gè)cosβ大小的因子，故獲得的該散射點(diǎn)在UOV中的坐標(biāo)為其實(shí)際值的cosβ倍，即=upcosβ,=vpcosβ。利用兩個(gè)坐標(biāo)系的變換關(guān)系式(13)可以得到

通過(guò)以上分析可知，該方法只需要對(duì)孤立強(qiáng)散射點(diǎn)的多普勒頻率信息進(jìn)行提取，不需要進(jìn)行復(fù)雜的迭代搜索，很適合用于本文圖像融合，其效果很大程度上依賴(lài)于對(duì)半雙基地角的估計(jì)精度(實(shí)際上也將轉(zhuǎn)動(dòng)角速度的估計(jì)包含其中)，估計(jì)誤差Δβ應(yīng)滿足：

其中φxmax和φymax分別表示目標(biāo)在兩個(gè)坐標(biāo)軸上的最大尺寸半徑。式(23)的意義在于，由估計(jì)誤差帶來(lái)的兩幅圖像相對(duì)的散射點(diǎn)偏移應(yīng)小于融合坐標(biāo)系下的分辨單元大小，否則融合就無(wú)法保證相干性，對(duì)應(yīng)的散射點(diǎn)錯(cuò)位疊加，使最終圖像散焦甚至嚴(yán)重失真。為了減小誤差帶來(lái)的影響，可以通過(guò)對(duì)對(duì)比度或熵等圖像整體信息的最優(yōu)化[15]來(lái)使估計(jì)值逼近實(shí)際值。

5 仿真結(jié)果與分析

本文采用圖1所示的雙基地ISAR轉(zhuǎn)臺(tái)模型進(jìn)行仿真驗(yàn)證，令坐標(biāo)原點(diǎn)為成像參考點(diǎn)，錄取數(shù)據(jù)過(guò)程中對(duì)參考點(diǎn)進(jìn)行精確跟蹤。假設(shè)半雙基地角近似為常數(shù)且β=18°，雷達(dá)發(fā)射信號(hào)中心頻率fc=9 GHz，帶寬B=300 MHz，采樣頻率Fs=400 MHz ，脈沖重復(fù)頻率Prf=1 00 Hz ，脈沖積累數(shù)N=1 28，目標(biāo)轉(zhuǎn)動(dòng)角速度ω=0 .03 rad/s ，則整個(gè)成像過(guò)程中目標(biāo)總轉(zhuǎn)角 Δθ≈ 2 .2°，那么雷達(dá)1，雷達(dá)2的距離向分辨率分別為ρy1=0.5 m,ρy2≈0 .526 m ，方位向分辨率分別為ρx1=0.434 m,ρx2≈ 0 .456 m 。在原始回波信號(hào)中添加高斯分布的白噪聲，使得回波信噪比為-3 dB。圖4所示為59點(diǎn)仿真目標(biāo)模型。

圖4 仿真目標(biāo)模型

將雷達(dá)1數(shù)據(jù)在時(shí)域分為兩個(gè)子孔徑并分別獲得其RD圖像，然后使用本文方法對(duì)轉(zhuǎn)動(dòng)角速度和半雙基地角進(jìn)行估計(jì)，結(jié)果如表1和表2所示。

對(duì)估計(jì)值分別取平均得到=0 .0301 rad/s,=1 8.1151°，可見(jiàn)估計(jì)精度較高。再利用估計(jì)得到的參數(shù)，按照式(15)完成融合處理，得到如下結(jié)果：

表1 目標(biāo)中3個(gè)散射點(diǎn)單元的轉(zhuǎn)角速度估計(jì)

表2 目標(biāo)中3個(gè)散射點(diǎn)單元的半雙基地角估計(jì)

圖5(a)和圖5(b)分別給出了雷達(dá)1和雷達(dá)2數(shù)據(jù)的BP算法等高線成像結(jié)果圖，可見(jiàn)目標(biāo)在噪聲背景的干擾下不能得到很好的辨識(shí)。而在對(duì)其數(shù)據(jù)進(jìn)行旋轉(zhuǎn)插值的校正以后，其圖像已與雷達(dá)1的圖像相互配準(zhǔn)，直接將數(shù)據(jù)相加就可以完成圖像融合。從以上結(jié)果可以看出，融合后的圖像噪聲與副瓣被抑制，信號(hào)能量得以積累加強(qiáng)，較大程度地改善了圖像的信噪比，令目標(biāo)特征信息更準(zhǔn)確。為了進(jìn)一步體現(xiàn)圖像融合的有效性，表3給出了幾幅圖像的熵和對(duì)比度的對(duì)比結(jié)果加以驗(yàn)證。

表3 融合前后圖像質(zhì)量對(duì)比

從表3可以看出，融合后的圖像對(duì)比度更高，圖像熵比融合前的圖像小，說(shuō)明其“銳化”程度較高，圖像質(zhì)量更好，證明了本文方法的有效性。

6 結(jié)束語(yǔ)

雙基地雷達(dá)系統(tǒng)能夠提供更大的相干積累角，獲得多視角下的目標(biāo)圖像和更豐富的目標(biāo)信息。本文采用雙基地 ISAR雷達(dá)轉(zhuǎn)臺(tái)模型，提出一種基于子孔徑參數(shù)估計(jì)的雙基地ISAR雷達(dá)一發(fā)一收體制下的圖像融合方法，該方法首先將其中一個(gè)雷達(dá)的回波分為兩個(gè)子孔徑，并提取孤立強(qiáng)散射單元在兩個(gè)孔徑下的多普勒頻率偏差，以此估計(jì)出目標(biāo)的轉(zhuǎn)動(dòng)角速度和半雙基地角，然后將兩個(gè)雷達(dá)數(shù)據(jù)進(jìn)行映射變換后完成融合，得到了質(zhì)量更高的目標(biāo)圖像。但本文的模型只適用于非機(jī)動(dòng)目標(biāo)和成像時(shí)間段較短的情況，并沒(méi)有考慮雙基地角時(shí)變效應(yīng)以及轉(zhuǎn)動(dòng)角速度的高階項(xiàng)的影響，在這種復(fù)雜情況下的雙基地ISAR成像以及圖像融合算法仍需要進(jìn)一步研究。

圖5 基于本文方法的雙基ISAR圖像融合結(jié)果

[1]Chen V C and Hao Ling.Time-Frequency Transforms for Radar Imaging and Signal Analysis[M].First Edition,Norwood,MA,Artech House,2002：1-6.

[2]Martorella M,Giusti E,Demi L,et al..Target recognition based by means of polarimetric ISAR images[J].IEEE Transactions on Aerospace and Electronic Systems,2011,47(1)：225-239.

[3]Chen V C,des Rosiers A,and Lipps R.Bi-static ISAR range-doppler imaging and resolution analysis[C].2009 IEEE Radar Conference,Pasadena,California,United States,2009：1-5.

[4]Martorella M,Palmer J,Homer J,et al..On bistatic inverse synthetic aperture radar[J].IEEE Transactions on Aerospace and Electronic Systems,2007,43(3)：1125-1134.

[5]Palmer J,Homer J,Longstaff I D,et al..ISAR imaging using an emulated multistatic radar system[J].IEEE Transactions on Aerospace and Electronic Systems,2005,41(4)：1464-1472.

[6]張亞標(biāo),朱振波,湯子躍,等.雙站逆合成孔徑雷達(dá)成像理論研究[J].電子與信息學(xué)報(bào),2006,28(6)：969-972.Zhang Ya-biao,Zhu Zhen-bo,Tang Zi-yue,et al..Bistatic inverse synthetic aperture radar image formation[J].Journal of Electronics＆Information Technology,2006,28(6)：969-972.

[7]高昭昭,梁毅,邢孟道,等.雙基地逆合成孔徑雷達(dá)成像分析[J].系統(tǒng)工程與電子技術(shù),2009,31(5)：1055-1059.Gao Zhao-zhao,Liang Yi,Xing Meng-dao,et al..Analysis of ISAR imagery for bistatic radar[J].Systems Engineering and Electronics,2009,31(5)：1055-1059.

[8]Hume A L and Baker C J.Netted radar sensing[C].2001 CIE International Conference on Radar,Beijing,China,2001：110-114.

[9]Li Zhi-xi,Papson S and Narayanan R M.Data-level fusion of multilook inverse synthetic aperture radar images[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(5)：1394-1406.

[10]毛士藝,趙巍.多傳感器圖像融合技術(shù)綜述[J].北京航空航天大學(xué)學(xué)報(bào),2002,28(5)：512-518.Mao Shi-yi and Zhao Wei.Comments on multisensor image fusion techniques[J].Journal of Beijing University of Aeronautics and Astronautics,2002,28(5)：512-518.

[11]Kim Yong-hyun,Lee Chang-no,Han Dong-yeob,et al..Improved additive-wavelet image fusion[J].IEEE Geoscience and Remote Sensing Letters,2011,8(2)：263-267.

[12]董健,尚朝軒,高梅國(guó),等.雙基地 ISAR 成像平面研究及目標(biāo)回波模型修正[J].電子與信息學(xué)報(bào),2010,32(8)：1855-1862.Dong Jian,Shang Chao-xuan,Gao Mei-guo,et al..The image plane analysis and echo model amendment of bistatic ISAR[J].Journal of Electronics＆Information Technology,2010,32(8)：1855-1862.

[13]Yeh C M,Xu J,Peng Y N,et al..Rotational motion estimation for ISAR via triangle pose difference on two range-Doppler images[J].IET Radar,Sonar and Navigation,2010,4(4)：528-536.

[14]Wang Yong and Jiang Yi-cheng.A novel algorithm for estimating the rotation angle in ISAR imaging[J].IEEE Geoscience and Remote Sensing Letters,2008,5(4)：608-609.

[15]Cao Pan,Xing Meng-dao,Sun Guang-cai,et al..Minimum entropy via subspace for ISAR autofocus[J].IEEE Geoscience and Remote Sensing Letters,2010,7(1)：205-209.