魏一搏,王 輝,李 寧,禹衛(wèi)東

(1.中國(guó)科學(xué)院電子學(xué)研究所,北京 100190; 2.中國(guó)科學(xué)院大學(xué),北京 100039)

適用于W波段ISAR成像的二維補(bǔ)償方法

魏一搏1,2,王 輝1,李 寧1,2,禹衛(wèi)東1

(1.中國(guó)科學(xué)院電子學(xué)研究所,北京 100190; 2.中國(guó)科學(xué)院大學(xué),北京 100039)

針對(duì)W波段逆合成孔徑雷達(dá)成像過(guò)程中越分辨單元徙動(dòng)嚴(yán)重的問(wèn)題,提出了一種二維補(bǔ)償成像方法.該方法首先分析了逆合成孔徑雷達(dá)成像中目標(biāo)徙動(dòng)產(chǎn)生的原因,然后使用Keystone變換對(duì)距離向徙動(dòng)進(jìn)行補(bǔ)償;對(duì)于均勻轉(zhuǎn)速目標(biāo),使用最小熵法從方位向徙動(dòng)中估計(jì)目標(biāo)的平均轉(zhuǎn)動(dòng)角速度,用相位補(bǔ)償?shù)姆椒▽?duì)均勻角速度目標(biāo)的方位向徙動(dòng)進(jìn)行補(bǔ)償;最后,對(duì)于變化角速度目標(biāo),利用子孔徑的方法降低目標(biāo)成像中的由角速度變化帶來(lái)的方位向徙動(dòng).該方法在仿真中獲得了高分辨率的成像結(jié)果,證明了該算法的有效性.

W波段;逆合成孔徑雷達(dá);Keystone變換;相位補(bǔ)償;子孔徑成像

逆合成孔徑雷達(dá)(Inverse Synthetic Aperture Radar,ISAR)是一種高分辨率微波成像系統(tǒng),它能夠提供空中、地面以及太空目標(biāo)的二維成像結(jié)果[1].近年來(lái),隨著針對(duì)W波段研究的深入以及W波段元器件的發(fā)展,W波段引起了世界各國(guó)尤其是發(fā)達(dá)國(guó)家的重視[2-3].W波段眾多的優(yōu)越特性,使得W波段ISAR技術(shù)有巨大的應(yīng)用前景.在成像中,相比傳統(tǒng)波段,W波段可以獲得厘米級(jí)的分辨率,使系統(tǒng)獲得更多的目標(biāo)細(xì)節(jié),提高對(duì)目標(biāo)的判別準(zhǔn)確率.

基于W波段的ISAR系統(tǒng),筆者針對(duì)成像過(guò)程中的徙動(dòng)影響,建立了轉(zhuǎn)動(dòng)目標(biāo)的ISAR回波相位模型.采用Keystone算法[4]進(jìn)行距離向的徙動(dòng)補(bǔ)償.對(duì)方位向的徙動(dòng),利用最小熵的方法估計(jì)得到轉(zhuǎn)動(dòng)的角速度.在此基礎(chǔ)上實(shí)現(xiàn)了方位向徙動(dòng)的補(bǔ)償.對(duì)于轉(zhuǎn)動(dòng)速度變化的目標(biāo),使用了子孔徑方法進(jìn)行成像.最后進(jìn)行了仿真驗(yàn)證.

1　ISAR成像模型

隨著W波段分辨率的提高,對(duì)于ISAR處理的補(bǔ)償精度要求進(jìn)一步提高.ISAR的補(bǔ)償可以分為運(yùn)動(dòng)補(bǔ)償和徙動(dòng)補(bǔ)償兩部分.運(yùn)動(dòng)補(bǔ)償是在對(duì)機(jī)動(dòng)目標(biāo)成像過(guò)程中產(chǎn)生的,對(duì)于機(jī)動(dòng)目標(biāo),需要對(duì)其在成像干涉過(guò)程中目標(biāo)中心點(diǎn)的距離向移動(dòng)進(jìn)行補(bǔ)償,使其等效于在同一個(gè)距離位置上進(jìn)行自旋轉(zhuǎn)動(dòng),以完成之后的成像處理.關(guān)于ISAR成像的運(yùn)動(dòng)補(bǔ)償在傳統(tǒng)波段已經(jīng)有了較為深入的研究,有許多成熟的算法[5-6],這些方法可以適用到W波段的成像處理中,因此筆者認(rèn)為輸入信號(hào)是經(jīng)過(guò)運(yùn)動(dòng)補(bǔ)償后的信號(hào).

這里根據(jù)文獻(xiàn)[2-3]中的系統(tǒng)參數(shù),假設(shè)W波段的ISAR成像系統(tǒng)分辨率為3 cm,則可以得到成像總轉(zhuǎn)角為0.05 rad.假設(shè)目標(biāo)尺寸25 m,則可以得到一階至三階的徙動(dòng)分別約為1.25 m,3.12 cm以及5.2 mm.由此可以看出,在W波段ISAR成像中,徙動(dòng)量比分辨率大許多,需要進(jìn)行補(bǔ)償處理才能得到較好的成像結(jié)果.同時(shí)進(jìn)行回波計(jì)算時(shí),采用二階的近似處理較為合理.

考慮只繞中心點(diǎn)進(jìn)行轉(zhuǎn)動(dòng)的目標(biāo),處理由于目標(biāo)轉(zhuǎn)動(dòng)帶來(lái)的越分辨單元徙動(dòng)所產(chǎn)生的圖像散焦問(wèn)題.目標(biāo)與雷達(dá)的關(guān)系如圖1所示.

圖1　ISAR成像幾何關(guān)系圖

雷達(dá)所能獲得的信息就是目標(biāo)上的點(diǎn)與雷達(dá)的距離R所帶來(lái)的回波時(shí)延信息.假設(shè)雷達(dá)與目標(biāo)中心點(diǎn)O的距離為R0,目標(biāo)上有一坐標(biāo)為(A,R)的點(diǎn),當(dāng)目標(biāo)轉(zhuǎn)動(dòng)角度為θ時(shí),采用fc+f的線性調(diào)頻信號(hào),假設(shè)目標(biāo)以恒定的角速度ω轉(zhuǎn)動(dòng)時(shí),回波的相位信息為

通過(guò)匹配濾波,由相位的第1項(xiàng)可以得到目標(biāo)點(diǎn)的距離向坐標(biāo).通過(guò)方位向的快速傅里葉變換(Fast Fourier Transform,FFT),由式(1)的第2項(xiàng),可以得到方位向的相對(duì)位置.式(1)的第1項(xiàng)和第4項(xiàng)則為目標(biāo)轉(zhuǎn)動(dòng)所帶來(lái)的距離向與方位向的越分辨單元徙動(dòng).如果不補(bǔ)償后兩項(xiàng),則會(huì)造成成像結(jié)果的模糊.

2　W波段成像特點(diǎn)及現(xiàn)有成像處理方法

W波段由于分辨率高、分辨單元小,所以越分辨單元的徙動(dòng)問(wèn)題比傳統(tǒng)波段要嚴(yán)重得多;同時(shí)隨著成像分辨率的提升,W波段圖像的數(shù)據(jù)量也比傳統(tǒng)波段大了許多.如果成像所使用算法的計(jì)算復(fù)雜度較高,會(huì)帶來(lái)計(jì)算時(shí)間的較大幅度增加.所以需要對(duì)成像算法的計(jì)算復(fù)雜度進(jìn)行考慮.

在現(xiàn)有的成像處理方法中,比較成熟的主要有距離多普勒(Range Doppler,RD)算法[1]、Keystone變換[4]、極坐標(biāo)格式算法(Polar Formation Algorithm,PFA)[7]、后向投影(Back Projection,BP)算法[8]、時(shí)頻方法[9]、運(yùn)動(dòng)參數(shù)估計(jì)法[10-11]等.但RD算法是不進(jìn)行補(bǔ)償,直接成像的方法,成像結(jié)果較差;Keystone算法能夠較好得補(bǔ)償目標(biāo)距離向的徙動(dòng),但并沒(méi)有處理方位向的徙動(dòng);PFA算法和BP算法需要得到目標(biāo)較為準(zhǔn)確的角速度才能成像,且算法的復(fù)雜度較高,耗時(shí)很長(zhǎng);時(shí)頻方法將目標(biāo)在方位向進(jìn)行時(shí)頻分析,能夠減輕方位向的徙動(dòng),但會(huì)帶來(lái)方位向分辨率的嚴(yán)重下降;低階的運(yùn)動(dòng)參數(shù)估計(jì)方法[10]只能夠?qū)D(zhuǎn)動(dòng)速度穩(wěn)定的目標(biāo)進(jìn)行補(bǔ)償,而高階的運(yùn)動(dòng)參數(shù)估計(jì)方法[11]具有較高的運(yùn)算復(fù)雜度.

3　二維補(bǔ)償成像處理

針對(duì)W波段成像的特點(diǎn),文中提出了一種在距離向與方位向分別進(jìn)行二維越分辨單元徙動(dòng)(Migration Through Resolution Cell,MTRC)補(bǔ)償?shù)某上穹椒?在距離向使用Keystone變換補(bǔ)償徙動(dòng),在方位向使用最小熵角速度估計(jì)與相位補(bǔ)償、子孔徑法相結(jié)合的方法進(jìn)行補(bǔ)償.該方法能夠在W波段得到較好的成像結(jié)果,同時(shí)該算法的運(yùn)算復(fù)雜度較低,計(jì)算所需時(shí)間較短.

3.1距離向MTRC補(bǔ)償

距離向MTRC可以利用Keystone方法進(jìn)行補(bǔ)償.Keystone變換是通過(guò)對(duì)回波信號(hào)在頻率域時(shí)間軸上的變標(biāo)操作來(lái)實(shí)現(xiàn)對(duì)距離向徙動(dòng)的補(bǔ)償?shù)?針對(duì)式(1),將tk=t(f+fc)/f帶入,得到

即補(bǔ)償了第3項(xiàng)的距離向徙動(dòng).

3.2均勻角速度目標(biāo)的方位向MTRC補(bǔ)償

為了進(jìn)行方位向的徙動(dòng)補(bǔ)償,可以進(jìn)一步展開(kāi)式(2)中的最后一項(xiàng),將tk=ak/fPRF,R=ρRr,ρR=c/2B帶入.其中的ak,r,ρr分別是圖像在Keystone變換后的方位向點(diǎn)數(shù)、距離向點(diǎn)數(shù)以及距離向的分辨率.fPRF為脈沖重復(fù)頻率(Pulse Repetition Frequency,PRF).可以得到

這里假設(shè)發(fā)射信號(hào)的相對(duì)帶寬較小,可以忽略.通過(guò)式(3)可以看出,方位向的徙動(dòng)帶來(lái)的相位差在圖像的某一點(diǎn)上,是只與目標(biāo)的轉(zhuǎn)動(dòng)速度有關(guān)的函數(shù),因此這里可以采用假定不同的轉(zhuǎn)動(dòng)角速度,得到成像結(jié)果,再利用最小熵搜索的方法來(lái)估計(jì)目標(biāo)的轉(zhuǎn)動(dòng)角速度.在搜索的過(guò)程中,可以采用多次搜索,并逐次遞增搜索精度的方法來(lái)降低計(jì)算量.角速度搜索的最大精度正比于成像結(jié)果的分辨率.在估計(jì)得到目標(biāo)的轉(zhuǎn)動(dòng)角速度后,給圖像上的每個(gè)點(diǎn)疊加由式(3)計(jì)算得到的相位,即可補(bǔ)償方位向的MTRC.

3.3均勻角速度目標(biāo)方位向MTRC補(bǔ)償誤差分析

假設(shè)角速度的估計(jì)有Δω的誤差,忽略信號(hào)相對(duì)帶寬以及Keystone變化對(duì)時(shí)間維造成的影響,由此,可得對(duì)方位向造成誤差的相位,即

其中,τ為成像時(shí)間,這里假設(shè)Δω相對(duì)ω為小量.

相位誤差在成像后造成的目標(biāo)點(diǎn)偏移為

其中,N為成像中方位向的總點(diǎn)數(shù).

由此可以看出,這一偏移是隨著目標(biāo)距離向離轉(zhuǎn)動(dòng)中心距離的增大而增大的,同時(shí)是隨著時(shí)間增大的.假設(shè)在經(jīng)過(guò)最小熵的搜索過(guò)程后,ΔN對(duì)于最大的偏移能夠達(dá)到最高的1/2格數(shù)補(bǔ)償精度.則角速度估計(jì)的誤差為

其中,RT為目標(biāo)距離向與轉(zhuǎn)動(dòng)中心的最大距離.

從式(6)可以看出,在W波段的成像中,較大的方位向MTRC為角速度的估計(jì)帶來(lái)了較高的精度.

3.4變化角速度目標(biāo)的方位向MTRC補(bǔ)償

以上所討論的都是基于目標(biāo)角速度不變的情況,然而在實(shí)際的應(yīng)用中,大部分目標(biāo)的角速度都會(huì)發(fā)生不同程度的變化,因此,需要考慮目標(biāo)角速度變化對(duì)結(jié)果帶來(lái)的影響.

假設(shè)角度隨時(shí)間變化的函數(shù)為ω(t),先對(duì)變化角速度的目標(biāo)進(jìn)行上文所述的均勻角速度目標(biāo)的補(bǔ)償處理.根據(jù)上文方法估計(jì)得到的平均角速度為可以得到在(A,R)點(diǎn),由于角速度變化所造成的相位誤差為

式(7)中的前一項(xiàng)代表著目標(biāo)角速度變化所造成的方位向的偏移,需要進(jìn)一步的補(bǔ)償.而如果能夠較為準(zhǔn)確得估計(jì)得到目標(biāo)角速度的平均值,則式(7)中后一項(xiàng)可以被之前所進(jìn)行的相位補(bǔ)償處理所補(bǔ)償.

如果目標(biāo)的角速度為ω(t),則在每個(gè)慢時(shí)間點(diǎn),目標(biāo)的方位向坐標(biāo)點(diǎn)為

從式(8)中可以看出,目標(biāo)的成像點(diǎn)會(huì)由于角速度的變化,隨著時(shí)間移動(dòng).同時(shí)這一偏移是隨著方位向的增大而增大的,在圖像的中心點(diǎn),偏移為零.

由于目標(biāo)的角速度隨著時(shí)間變化,在整個(gè)合成孔徑時(shí)間內(nèi),目標(biāo)在方位向會(huì)不斷進(jìn)行偏移.這里假設(shè)目標(biāo)的角速度是連續(xù)的,則在一小段時(shí)間內(nèi),目標(biāo)的角速度變化會(huì)較小,即如果只對(duì)一小段孔徑進(jìn)行成像,則目標(biāo)在方位向由于角速度變化所帶來(lái)的徙動(dòng)也會(huì)較小.由此可以通過(guò)將回波數(shù)據(jù)分成若干個(gè)子孔徑分別進(jìn)行成像,用再合成的方法來(lái)降低角速度變化帶來(lái)的方位向徙動(dòng).假設(shè)第K個(gè)子孔徑的平均角速度為,可以得到目標(biāo)上的(A,R)點(diǎn)方位向坐標(biāo)為

從式(10)可以看出,只需要對(duì)各子孔徑圖像進(jìn)行方位向的伸縮處理就可以得到各目標(biāo)點(diǎn)坐標(biāo)完全一致的子圖像.得到子圖像后只需將各子圖像轉(zhuǎn)換回時(shí)域,疊加后再重新成像即可得到合成的圖像.通過(guò)子孔徑的處理可以有效地降低角速度變化所帶來(lái)的影響.同時(shí)由于子孔徑方法充分利用了各子孔徑的信息進(jìn)行相干疊加,圖像的分辨率與原圖像的相同.

在實(shí)際的處理過(guò)程中,由于各子孔徑平均角速度的估計(jì)誤差較大,可以采用選取某個(gè)子孔徑為基準(zhǔn),對(duì)其他子孔徑進(jìn)行對(duì)齊的方法進(jìn)行.為了獲得某個(gè)子孔徑相對(duì)基準(zhǔn)孔徑的伸縮倍數(shù),可以設(shè)置不同的伸縮倍數(shù),對(duì)兩個(gè)孔徑進(jìn)行合成成像,根據(jù)所得圖像的熵進(jìn)行搜索的方法進(jìn)行.

通過(guò)使用子孔徑的方法,可以降低目標(biāo)角速度變化所帶來(lái)的方位向徙動(dòng).這種方法可以適用于一定幅度內(nèi)任意變化的角速度.并且這種方法是由成像結(jié)果進(jìn)行參數(shù)的選取,具有較強(qiáng)的抗干擾性能.這種方法在每次的搜索過(guò)程中,只需要對(duì)子孔徑進(jìn)行方位向的傅里葉變換以及伸縮處理.伸縮處理可以使用基于Chirp-Z變換的方法實(shí)現(xiàn),因此每次搜索的運(yùn)算復(fù)雜度較低、耗時(shí)較少.

4　仿真驗(yàn)證與分析

4.1均勻角速度目標(biāo)成像

由于W波段設(shè)備的成像數(shù)據(jù)現(xiàn)階段較難獲得,因此文中采用了仿真的方式來(lái)驗(yàn)證算法的正確性.所采用的系統(tǒng)參數(shù)如表1所示.

表1　系統(tǒng)仿真參數(shù)列表

使用RD成像算法得到的結(jié)果如圖2(a)所示.距離壓縮后,使用Keystone方法進(jìn)行了距離向補(bǔ)償?shù)玫降慕Y(jié)果如圖2(b)所示,利用角速度搜索結(jié)果進(jìn)行相位補(bǔ)償后的成像結(jié)果如圖2(c)所示.可以看出對(duì)于均勻轉(zhuǎn)速的目標(biāo),已經(jīng)完全補(bǔ)償了目標(biāo)的徙動(dòng),得到了清晰的圖像.

圖2　均勻角速度目標(biāo)二維補(bǔ)償效果圖

4.2變化角速度目標(biāo)成像

當(dāng)目標(biāo)的角速度發(fā)生變化時(shí),會(huì)造成額外的多普勒頻率變化,繼而引起方位向的成像位置移動(dòng),引起成像的方位向模糊.仿真中除目標(biāo)角速度外,均采用表1中所示的參數(shù).設(shè)置目標(biāo)的角速度具有二階以及三階的角加速度.

使用均勻角速度成像算法得到的結(jié)果如圖3(a)所示.將得到的結(jié)果分成4個(gè)子孔徑進(jìn)行處理,經(jīng)過(guò)搜索以及進(jìn)一步補(bǔ)償成像后,得到的結(jié)果如圖3(b)所示.從圖中可以看出,成像質(zhì)量雖然比沒(méi)有角加速度時(shí)的圖2(c)有一定的下降,但與圖3(a)相比,降低了角加速度所帶來(lái)的方位向模糊,得到了較好的成像結(jié)果.

圖3　子孔徑補(bǔ)償效果圖

4.3算法耗時(shí)

文中仿真計(jì)算處理的場(chǎng)景總點(diǎn)數(shù)為1 000×1 000,使用的是Core i7-2670 QM 2.2 GHz處理器.這里列出了對(duì)于變化角速度目標(biāo)計(jì)算過(guò)程的耗時(shí).算法各部分的計(jì)算時(shí)間如表2所示.

表2　各模塊計(jì)算時(shí)間s

從表2中可以看出,對(duì)于這一尺寸的場(chǎng)景來(lái)說(shuō),算法的總耗時(shí)是可以接受的.

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

筆者研究了W波段ISAR越分辨單元徙動(dòng)的產(chǎn)生原因,以及針對(duì)均勻角速度和變化角速度目標(biāo)的補(bǔ)償方法.在此基礎(chǔ)上進(jìn)行了成像,并進(jìn)行了仿真驗(yàn)證.研究表明:對(duì)于W波段的ISAR成像,距離向的徙動(dòng)可以通過(guò)Keystone算法進(jìn)行補(bǔ)償.對(duì)于均勻轉(zhuǎn)速的目標(biāo),方位向的徙動(dòng)可以通過(guò)最小熵相位補(bǔ)償方法進(jìn)行補(bǔ)償,在W波段下,補(bǔ)償?shù)恼`差較小.對(duì)于轉(zhuǎn)速變化的目標(biāo),使用子孔徑與相位補(bǔ)償相結(jié)合的方法可以得到較好的成像結(jié)果.與此同時(shí),文中所提出的算法具有較低的運(yùn)算復(fù)雜度,耗時(shí)較少.

[1]AUSHERMAN D A,KOZMA A,WALKER J L,et al.Developments in Radar Imaging[J].IEEE Transactions on Aerospace and Electronic Systems,1984,20(4):363-400.

[2]LINCOLN LABORATORY.MIT Lincoln Laboratory 2009 Annual Report 2009[EB/OL].[2015-1-12].http://www. ll.mit.edu/publications/Annual?Report?2009.pdf.

[3]DANLY B G,CHEUN J,GREGERS-HANSEN V.WARLOC:a High-power Millimeter-wave Radar[C]//Proceedings of Twenty Seventh International Conference on Infrared and Millimeter Waves.Piscataway:IEEE,2002:233-234.

[4]陳文馳,保錚,邢孟道,等.基于Keystone變換的低信噪比ISAR成像[J].西安電子科技大學(xué)學(xué)報(bào),2003,30(2): 155-159. CHEN Wenchi,BAO Zheng,XING Mengdao,et al.Keystone Transformation Based ISAR Imaging at the Low SNR Level[J].Journal of Xidian University,2003,30(2):155-159.

[5]黃大榮,張磊,邢孟道,等.機(jī)動(dòng)目標(biāo)ISAR自聚焦新方法[J].西安電子科技大學(xué)學(xué)報(bào),2014,41(3):71-78. HUANG Darong,ZHANG Lei,XING Mengdao,et al.Autofocus Method for Maneuvering Targets[J].Journal of Xidian University,2014,41(3):71-78.

[6]YU X,ZHU D Y,ZHANG J D,et al.Motion Compensation Algorithm Based on the Designing Structured Gram Matrices Method[J].IET Radar,Sonar&Navigation,2014,8(3):209-219.

[7]WALKER J L.Range-Doppler Imaging of Rotating Objects[J].IEEE Transactions on Aerospace and Electronic Systems,1980,16(1):23-52.

[8]MENSA D L,HALEVY S,WADE G.Coherent Doppler Tomography for Microwave Imaging[J].Proceedings of the IEEE,1983,71(2):254-261.

[9]李亞超,蘇軍海,邢孟道,等.利用時(shí)間-調(diào)頻率分布特性的復(fù)雜運(yùn)動(dòng)目標(biāo)ISAR成像研究[J].西安電子科技大學(xué)學(xué)報(bào),2008,35(1):1-7. LI Yachao,SU Junhai,XING Mengdao,et al.Research on ISAR Imaging of the Complexly Moving Target Based on the Time-chirp Distribution[J].Journal of Xidian University,2008,35(1):1-7.

[10]XU Z W,ZHANG L,XING M D,et al.Uniform Rotational Motion Compensation for Inverse Synthetic Aperture Radar Targets from Image Domain[C]//IEEE China Summit and International Conference on Signal and Information Processing Conference.Piscataway:IEEE,2014:166-170.

[11]LI Y C,XING M D,SU J H,et al.A New Algorithm of ISAR Imaging for Maneuvering Targets with Low SNR[J]. IEEE Transactions on Aerospace and Electronic Systems,2013,49(1):543-557.

(編輯:王 瑞)

Two dimension migration compensation method for W-band ISAR imaging

WEI Yibo1,2,WANG Hui1,LI Ning1,2,YU Weidong1
(1.Institute of Electronics,Chinese Academy of Sciences,Beijing 100190,China; 2.University of Chinese Academy of Sciences,Beijing 100039,China)

The problem of migration through the resolution cell(MTRC)in the W-band Inverse Synthetic Aperture radar(ISAR)imaging system is serious.Therefore,a novel two dimensional compensation method is proposed in this paper.Firstly,the cause of migration in ISAR imaging is analyzed.Then,the Keystone transform is used to compensate the MTRC in the range direction.After that,the minimum entropy algorithm is used to estimate the average rotate speed of the target and the rotate speed is used to compensate the cross range direction MTRC of a uniformly rotating target by using the phase correcting method.Finally,the sub-aperture imaging algorithm is used to reduce the cross range MTRC in imaging the maneuvering target.Simulation results show that a high resolution result is achieved and the effectiveness of the two dimensional compensation method is verified.

W-band;ISAR;keystone transform;phase correction;sub-aperture imaging

TN955+.1

A

1001-2400(2016)03-0155-06

10.3969/j.issn.1001-2400.2016.03.027

2015-01-19

時(shí)間:2015-07-27

魏一搏(1990-),男,中國(guó)科學(xué)院電子學(xué)研究所碩士研究生,E-mail:ghost200802@163.com.

王 輝(1972-),女,副研究員,E-mail:wanghui@mail.ie.ac.cn.

http://www.cnki.net/kcms/detail/61.1076.TN.20150727.1952.027.html