郭勝龍李 洋楊士林周勇勝洪 文

①(中國科學(xué)院電子學(xué)研究所 北京 100190)

②(微波成像技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室 北京 100190)

③(中國科學(xué)院大學(xué) 北京 100190)

④(中國科學(xué)院光電研究院 北京 100094)

基于模型的π/4模式簡縮極化干涉數(shù)據(jù)目標(biāo)分解

郭勝龍①②③李 洋①②③楊士林①②③周勇勝④洪 文*①②

①(中國科學(xué)院電子學(xué)研究所 北京 100190)

②(微波成像技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室 北京 100190)

③(中國科學(xué)院大學(xué) 北京 100190)

④(中國科學(xué)院光電研究院 北京 100094)

基于模型的目標(biāo)分解是極化合成孔徑雷達(dá)(SAR)的一個重要應(yīng)用，基于模型目標(biāo)分解依賴于極化數(shù)據(jù)僅能獲得各散射機(jī)制的能量信息。該文將基于模型的分解技術(shù)應(yīng)用到π/4模式的簡縮極化干涉SAR(C-PolInSAR)數(shù)據(jù)中，對互相關(guān)矩陣進(jìn)行分解，在獲得各散射機(jī)制功率的同時獲得對應(yīng)的散射相位中心。該文首先推導(dǎo)出3種散射機(jī)制π/4簡縮極化SAR干涉觀測下散射模型，然后運(yùn)用數(shù)值計(jì)算方法進(jìn)行目標(biāo)分解，最終求解出各散射機(jī)制的功率貢獻(xiàn)及相位中心高度信息。仿真數(shù)據(jù)驗(yàn)證了該算法的有效性，分析了不同波段及不同地表參數(shù)對分解結(jié)果的影響。

合成孔徑雷達(dá)(SAR)；π/4模式簡縮極化；簡縮極化干涉SAR；模型目標(biāo)分解

1 引言

2005年文獻(xiàn)[1]首次提出π/4簡縮極化SAR模式，該模式為發(fā)射1路由線性水平極化(H)與垂直極化(V)合成的45°線性極化波，同時接收H, V兩路極化波。該模式因其只發(fā)射1種極化波且僅獲取2通道數(shù)據(jù)成功克服全極化合成孔徑雷達(dá)系統(tǒng)帶寬小，負(fù)荷重等[2]的缺點(diǎn)，同時又獲得比單極化、雙極化較多的極化信息。文獻(xiàn)[1]根據(jù)自然地物呈現(xiàn)的散射特性提出反射對稱性及同極化與交叉極化之間的關(guān)系，成功地重建出L波段全極化(偽全極化)數(shù)據(jù)[1]，且地物分類結(jié)果的精確性與全極化數(shù)據(jù)一致。文獻(xiàn)[3]采用P波段π/4簡縮極化數(shù)據(jù)及極化干涉理論對森林區(qū)域進(jìn)行分類，獲得與全極化數(shù)據(jù)相似的分類結(jié)果。此外還有一系列針對該模式的研究。因此對簡縮極化研究具有深遠(yuǎn)意義。除此模式外，簡縮極化還有圓極化發(fā)射線極化接收模式[4,5]及圓極化發(fā)射圓極化接收模式[6]。本文運(yùn)用模式進(jìn)行分析與檢驗(yàn)該模式下的極化干涉目標(biāo)分解能力，同時也推導(dǎo)出該模式下體散射功率求解方法。

該文將基于模型分解技術(shù)[7,8]應(yīng)用到π/4模式簡縮極化干涉數(shù)據(jù)中。首先推導(dǎo)出3種散射機(jī)制(單散射，二次散射及體散射)在π/4模式簡縮極化干涉觀測下的散射模型，然后對簡縮極化干涉數(shù)據(jù)的互相關(guān)觀測矩陣進(jìn)行數(shù)值方法分解。數(shù)值求解前首先根據(jù)簡縮極化Stokes矢量的地表上均勻體(Random Volume on Ground, RVoG)散射模型[9]求解體散射幅度，其次根據(jù)外部數(shù)據(jù)求解地表相位，最后數(shù)值求解出3種散射模型的加權(quán)系數(shù)。加權(quán)系數(shù)即包含每種散射機(jī)制的功率信息又包括散射相位中心高度的信息。因此通過上述分解技術(shù)在獲取3種散射機(jī)制的功率貢獻(xiàn)的同時還能獲得各散射機(jī)制的散射相位中心高度。

2 π/4模式簡縮極化干涉

π/4簡縮極化模式是指發(fā)射45°線極化波(可通過控制H, V極化波相位差為0得到)，接收H, V極化波。該模式的簡縮極化數(shù)據(jù)可由全極化復(fù)散射矩陣S得到[1]，如式(1)所示。

其中，i=1,2分別代表基線兩端主圖像與從圖像。因此，對應(yīng)的互相關(guān)矩陣J可以表示為

其中，J為2×2互相關(guān)復(fù)矩陣既包含極化信息又包含干涉信息。將J表示成3個2×2矩陣之和，如式(4)所示，3個矩陣分別對應(yīng)π/4簡縮極化干涉下單散射Js，二次散射Jd與體散射Jvol。

其中fs,fd,fvol分別3種散射的加權(quán)系數(shù)。加權(quán)系數(shù)為復(fù)數(shù)，幅度代表該散射機(jī)制的功率貢獻(xiàn)，相位則代表該散射相位中心。在已知地表相位及干涉的垂直波束kz情況下，可進(jìn)一步求出各種散射機(jī)制發(fā)生的豎直高度。

3 模型建立

接下來，本文首先重新推導(dǎo)3種散射機(jī)制在π4簡縮極化模式下的散射模型，模型著重引入一些由干涉處理引起的新的變量及形式。

3.1 單散射模型

在全極化中，單散射模型是建立在Bragg 條件或者小尺度模型(Small Perturbation Model, SPM)[9]下，小尺度模型為假如一個表面滿足k· s ＜0.3和k· l＜0.3兩個條件，我們稱該表面相對光滑。這里k指雷達(dá)波束，s指表面起伏高度均方誤差，l指相關(guān)長度。該小尺度模型假定交叉通道hv=0。因此Bragg參數(shù)僅取決于地表介電常數(shù)與波入射角。由此假定主圖像與從圖像相同極化通道幅值相同，相位差僅由干涉相位引起。即

單散射機(jī)制的全極化散射矩陣Ss可表示為

其中i=1,2。于是由式(1)～式(3)可得單散射模型互相關(guān)矩陣為式(7)。

根據(jù)式(5)的假設(shè)可得Js矩陣的4個元素為

3.2 二次散射模型

二次散射模型建立在二面角反射器的散射特性

基礎(chǔ)上[10]?；€兩端二次散射機(jī)制全極化散射矩陣Sd為

其中i=1,2,Rth, Rgh, Rtv, Rgv分別為不同極化下豎直物體及地表的Fresnel系數(shù)。γv,γh為不同極化波傳輸復(fù)系數(shù)，反映波在傳輸過程中能量的衰減及相位改變。由于Fresnel系數(shù)僅依賴于地物及入射波，則對主、從圖像來說相同極化通道Fresnel系數(shù)相同。又假定基線兩端h,v極化波傳播過程能量衰減相同及相同極化通道相位改變僅由干涉相位引起，由此可見基線兩端二次散射仍滿足式(5)假定。

由式(10)可知，其與單散射機(jī)制的散射矩陣結(jié)構(gòu)相同又滿足式(5)假定。同理，可得簡縮極化干涉觀測模式下的歸一化二次散射模型，如式(11)所示。

3.3 體散射模型

帶自樹木冠層的體散射模型是建立在一團(tuán)隨機(jī)指向圓柱類型的散射子的散射貢獻(xiàn)基礎(chǔ)上[10]。體散射的全極化散射矩陣Svol如式(12)所示。其中θ角為圍繞雷達(dá)視線的旋轉(zhuǎn)角。

考慮到體散射元素可近似為極細(xì)的水平向小圓柱體[10]即為方便推導(dǎo)可令由此式(12)可化簡為

對體散射而言，所有極化通道內(nèi)體散射相位中心都相同，主、從圖像相位差僅由干涉相位引起且相同極化通道的體散射幅度相同。

式中，m=1,2,n=1,2分別代表矩陣的行與列。最終推導(dǎo)出在簡縮極化干涉觀測模式下的歸一化均勻體散射模型為

4 基于模型的目標(biāo)分解

在第3節(jié)，本文推導(dǎo)出3種散射機(jī)制在π4簡縮極化干涉觀測下的模型。將3種散射模型代入式(4)可得

該方程組中包括5個復(fù)未知量(10個實(shí)數(shù)未知量)：fs,fd,fvol,α, β，而僅有4個復(fù)觀測量(8個實(shí)數(shù)觀測量)，即J的4個元素：J11,J12,J21,J22。由引可見該方程組是一個欠定非線性方程組，為此，在用數(shù)值方法求解之前，先求解2個量，剩余8個未知量通過8個方程數(shù)值求解。需要說明J11是基線兩端HL(線極化波發(fā)射，H通道接收)通道互相關(guān)數(shù)據(jù)，J22是基線兩端VL(線極化波發(fā)射，V通道接收)通道互相關(guān)數(shù)據(jù)，因此下文用HL通道，VL通道分別代指J11,J22兩通道。

4.1 體散射幅度的求解

體散射幅度|fvol|可根據(jù)簡縮極化Stokes矢量的RVoG模型理論求得。根據(jù)RVoG理論，觀測數(shù)據(jù)的Stokes矢量可表示成體散射模型和秩為1的純目標(biāo)(單散射目標(biāo)或者二次散射目標(biāo))模型之和。

4.2 二次散射相位中心求解

在森林區(qū)域，偶次散射主要來自回波在地表與樹干之間的相互作用，偶次散射相位中心近似位于地形表面[11]，因此我們通過估計(jì)地形相位φtopo來近似為fd相位φd。在有先驗(yàn)信息情況下或者有外部數(shù)據(jù)時，如數(shù)字高程模型(Digital Elevation Model, DEM)數(shù)據(jù)，地形相位可以直接得到。當(dāng)缺乏有關(guān)地形信息時，可利用全極化數(shù)據(jù)依靠RVoG模型[12]進(jìn)行地形相位估計(jì)。該文實(shí)驗(yàn)數(shù)據(jù)為PolSARpro[13]仿真全極化干涉數(shù)據(jù)，仿真數(shù)據(jù)地表設(shè)置為平地?zé)o起伏，因此森林周邊裸地相位作為整個森林?jǐn)?shù)據(jù)的形相位。

4.3 非線性方程組求解

其中，F(xiàn)d為二次散射的加權(quán)系數(shù)fd的幅度,φvol為體散射加權(quán)系數(shù)fvol的相位。通過求解滿足式(19)最小誤差的數(shù)值解，可得fs,fd,fvol,α, β。根據(jù)求解出的5個復(fù)參數(shù)，可求解出3種散射機(jī)制分別在HL 互相關(guān)通道和VL互相關(guān)通道的功率貢獻(xiàn)。以地表為參考點(diǎn)，在已知豎直波束參數(shù)kz條件下，3種散射機(jī)制的相位中心高度由式(20)計(jì)算出[14]。式(20)中type代表單散射，二次散射與體散射3種散射機(jī)制。

接下來本文分別用L波段森林仿真數(shù)據(jù)(不同地表參數(shù))及實(shí)測數(shù)據(jù)對上述理論進(jìn)行驗(yàn)證與分析。

5 實(shí)驗(yàn)結(jié)果與分析

5.1 L波段，地表粗糙

本文用PolSARProSim[13]軟件仿真L波段(1.3 GHz)全極化干涉數(shù)據(jù)，然后根據(jù)式(1)～式(3)合成π/4模式簡縮極化干涉數(shù)據(jù)進(jìn)行分解實(shí)驗(yàn)。仿真森林高度為10 m。地面無坡度地表參數(shù)高粗糙度。水平基線10 m，豎直基線1 m。仿真場景的Pauli圖如圖1所示。圖中紅線代表下文分析所取的一條方位線數(shù)據(jù)。

圖 1 仿真場景Pauli分解圖

3種散射機(jī)制的幅度分解結(jié)果如圖2所示。由圖2(a)可知，VL通道內(nèi)，森林區(qū)域體散射占主導(dǎo)地位，單散射較低于體散射，個別點(diǎn)處單散射與體散射功率相同，統(tǒng)計(jì)分析單散射平均功率值低于體散射 5.80 dB 如表1所示，二次散射能量最低(-6.73 dB)。由圖2(b)HL通道內(nèi)，體散射依然占主導(dǎo)地位，單散射次之，二次散射最低。與VL通道相比單散射與二次散射功率都有所提高(單散射提高1.61 dB，二次散射提高2.49 dB)如表1前2行數(shù)據(jù)所示。這一結(jié)果與Fresnel系數(shù)定義相符合，即H通道的貢獻(xiàn)要高于V通道貢獻(xiàn)[14]。需要注意的是，由于地表粗糙，VL通道二次散射能量較低，且HL通道內(nèi)也不占主導(dǎo)地位[15](這與圖3形成鮮明對比，圖3由于地表光滑則HL通道二次散射占主導(dǎo)地位)。3種散射機(jī)制的相位中心高度分解結(jié)果如圖4所示。由圖4可以清晰看到，在森林區(qū)域單散射相位中心在森林冠層，這就說明單散射主要來自于森林冠層而非來自地表。在裸地區(qū)域，散射相位中心都在地表處。值得注意的是，由于在求解方程時已假設(shè)二次散射相位中心就是地表相位[11]，所以，以地表為參考，二次散射中心值為0。體散射相位中心位于樹高一半與樹高之間(即5～10 m)，符合RVoG模型樹高估計(jì)[16]。

表1 π/4 模式簡縮極化干涉目標(biāo)分解森林區(qū)域功率值統(tǒng)計(jì)

5.2 L波段，地表光滑

圖2 3種散射機(jī)制功率貢獻(xiàn)參數(shù)：L波段，地表粗糙

圖3 3種散射機(jī)制功率貢獻(xiàn)，參數(shù)：L波段，地表光滑

現(xiàn)在改變仿真參數(shù)，將地表參數(shù)改為光滑，其它仿真參數(shù)不變。分解結(jié)果如圖3所示。與圖2相比，圖3中，對于VL通道森林區(qū)域體散射依然占主導(dǎo)地位，單散射能量下降(下降5.78 dB，表1中1, 3行所示)，這是因?yàn)榈乇砉饣R面反射較強(qiáng)，從而引起單次散射能量下降；由于地表光滑，二次散射量也下降，符合Fresnel對光滑地表的定義[15]。對于HL通道，體散射與單散射沒有大的變化，然而二次散射量(13.21 dB)明顯增高，這主要是因?yàn)楣饣乇砼c垂直樹干之間引起強(qiáng)烈的二次散射。在裸地，體散射占主導(dǎo)地位主要由噪聲引起，由于噪聲也是雜亂無章與均勻體散射一樣表現(xiàn)為去極化特性，因此噪聲分量被分解為體散射分量。3種散射機(jī)制相位中心高度如圖5所示。與圖4相比無明顯改變，森林單散射來自森林冠層，二次散射發(fā)生在地表，體散射相位中心位于5～10 m之間。裸地處3種機(jī)制均發(fā)生在地表。

5.3 L波段實(shí)測數(shù)據(jù)

最后用L波段BioSAR2008數(shù)據(jù)集驗(yàn)證以上算法。數(shù)據(jù)采集地點(diǎn)為瑞典Vindeln北部森林區(qū)，我們選擇水平基線為6 m的一對干涉數(shù)據(jù)。測試區(qū)域的Pauli分解偽彩色合成圖如圖6所示。首先根據(jù)式(1)～式(3)合成混合極化干涉數(shù)據(jù)，然后進(jìn)行分解試驗(yàn)。目標(biāo)分解結(jié)果如圖7，圖8所示。

3種散射機(jī)制的幅度貢獻(xiàn)分解結(jié)果如圖8所示，功率值統(tǒng)計(jì)如表1所示。對VL互相關(guān)通道，體散射占主導(dǎo)地位(12.41 dB)，單散射次于體散射，整個森林區(qū)域二次散射都非常低(-5.83 dB)，這由實(shí)際地表粗糙引起的，實(shí)測數(shù)據(jù)的分解結(jié)果與仿真數(shù)據(jù)結(jié)果一致。在像素160附近，對應(yīng)圖 8中為森林中間的一片空地，此處體散射下降，單散射有一個較小的上升。在像素310附近是森林中的道路，同樣體散射下降。此外，也發(fā)現(xiàn)在道路邊緣，二次散射功率上升，這是由于地表與森林樹干產(chǎn)生較強(qiáng)的二次散射。對于HL互相關(guān)通道，在森林區(qū)域，體散射仍占主導(dǎo)地位。與VL通道相比，二次散射功率明顯上升(上升11.28 dB左右)，這與Fresnel系數(shù)定義一致，二次散射的Fresnel系數(shù)H通道要高于V通道，也與仿真數(shù)據(jù)結(jié)果一致。

圖8為獲取的3種散射機(jī)制相位中心高度圖。如圖8所示，二次散射相位中心高度在地表，因此二次散射主要發(fā)生在地表，單散射主要來源于森林上層，除了個別像素點(diǎn)160, 210, 310, 510以外。在像素160附近是森林中裸地，像素310附近為道路，其它像素位置應(yīng)為森林中的空隙，因此在這些像素點(diǎn)附近單散射發(fā)生在地表。在大片裸地(像素160附近)以及林中道路(像素310附近)，由于無樹木，則體散射相位中心落在地表。

圖4 以地表為參考3種散射機(jī)制相位中心高度，參數(shù)：L波段，地表粗糙

圖5 以地表為參考3種散射機(jī)制相位中 心高度，參數(shù)：L波段，地表光滑

圖6 L波段Biomass數(shù)據(jù)Pauli分解圖

圖7 3種散射機(jī)制功率貢獻(xiàn)， 參數(shù)：L波段實(shí)測數(shù)據(jù)

圖8 以地表為參考3種散射機(jī)制相位 中心高度，參數(shù)：L段波實(shí)測數(shù)據(jù)

6 結(jié)束語

本文將基于模型的目標(biāo)分解技術(shù)應(yīng)用到π/4簡縮極化干涉數(shù)據(jù)中，通過對簡縮極化干涉互相關(guān)矩陣進(jìn)行分解，獲得3種散射機(jī)制的功率的同時還能獲得各散射機(jī)制的散射相位中心高度信息。通過分析不同地表參數(shù)不同波段數(shù)據(jù)以及實(shí)測數(shù)據(jù)的分解結(jié)果，驗(yàn)證了該分解技術(shù)的正確性與可行性。這也說明該分解技術(shù)能夠精確地適應(yīng)于不同地表參數(shù)及不同波段的數(shù)據(jù)。

[1] Souyris J C, Imbo P, Fjortoft R, et al.. Compact polarimetry based on symmetry properties of geophysical media: the pi/4 mode[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(3): 634-646.

[2] 陳琳. 簡縮極化SAR信息處理模型與方法研究[D]. [博士論文], 中國科學(xué)院電子學(xué)研究所, 2013. Chen Lin. Investigation on models and methods of compact polarimetric SAR information processing[D]. [Ph.D. dissertation], Institute of Electronics, Chinese Academy of Sciences, 2013.

[3] Dubois-Fernandez P, Anglliaume S, Souyris J C, et al.. The specificity of P-band PolInSAR data over vegetation[C]. Proceedings of ESA PoLinSAR Workshop, Frascati, 2007: 256-259.

[4] Raney R K. Hybrid-polarity SAR architecture[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(11): 3397-3404.

[5] Collins M J, Denbina M, and Atteia G. On the reconstruction of quad-Pol SAR data from compact polarimetry data for ocean target detection[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 591-600.

[6] Nord M E, Ainsworth T L, Lee J S, et al.. Comparison of compact polarimetric synthetic aperture radar modes[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 174-188.

[7] Ballester-Berman J D and Lopez-Sanchez J M. Applying the Freeman-Durden decomposition concept to polarimetric SAR Interferometry[J]. IEEE Transactions on Geoscience Remote Sensing, 2010, 48(1): 466-479.

[8] Jiao Zhi-hao, Yang Jian, Yeh C, et al.. Modified threecomponent decomposition method for polarimetric SAR data[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(1): 200-204.

[9] Chen H, Cloude S R, and Goodenough D G. Compact radar polarimetry for fire scar mapping and land cover classification[C]. Proceedings of the IET International Radar Conference, Xi’an, China, 2013: 1-4.

[10] Tsang L, Kong J A, and Shin R T. Theory of Microwave Remote Sensing[M]. New York: Wiley Interscience, 1985: 4-6.

[11] Wilsen C B, Sarabandi K, and Lin Y C. The effect of treearchitecture on the polarimetric and interferometric radar response[C]. Proceedings of the International Geoscience and Remote Sensing Symposium, Seattle, USA, 1998: 1499-1501.

[12] 白璐. 基于多層模型的PolInSAR 垂直結(jié)構(gòu)參數(shù)估計(jì)方法研究及其在森林覆蓋區(qū)的應(yīng)用[D]. [博士論文], 中國科學(xué)院大學(xué), 2009. Bai Lu. Multiple-layer model based PolInSAR methods to estimate vertical structure parameters and their applications in forest areas[D]. [Ph.D. dissertation], University of Chinese Academy of Science, 2009.

[13] Pottier E. ProSARPro_v4.2.0 [OL]. http://earth.esa.int/ polsarpro/default.html. 2012.

[14] Cloude S R. Polarisation: Applications in Remote Sensing [M]. London: UK, Oxford University Press, 2009: 211-212.

[15] Lecki D G and Ranson K J. Manual of Remote Sensing: Principles and Applications of Imaging Radar[M]. 3rd edition. New York: Wiley, 1998: 9.

[16] Lavalle M. Full and compact polarimetric radar interferometry for vegetation remote sensing[D]. [Ph.D. dissertation], University of Rennes, 2009.

郭勝龍： 男，1987年生，博士生，研究方向?yàn)闃O化SAR/簡縮極化SAR理論研究與應(yīng)用.

李 洋： 男，1983年生，博士生，研究方向?yàn)闃O化SAR、簡縮極化SAR、相干層析理論研究與應(yīng)用.

楊士林： 男，1989年生，碩士生，研究方向?yàn)榛旌蠘O化SAR系統(tǒng)仿真與設(shè)計(jì).

周勇勝： 男，1982年生，博士，研究方向?yàn)闃O化干涉SAR機(jī)理、機(jī)制與植被信息提取.

洪 文： 女，1968年生，教授，博士生導(dǎo)師，研究方向?yàn)镾AR成像與系統(tǒng)及其應(yīng)用、極化/極化干涉SAR數(shù)據(jù)處理及應(yīng)用、3維微波成像新概念新體制新方法等.

Model-based Target Decomposition with π/4Mode Compact Polarimetric SAR Interferometry Data

Guo Sheng-long①②③Li Yang①②③Yang Shi-lin①②③Zhou Yong-sheng④Hong Wen①②①(Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China)
②(Science and Technology on Microwave Imaging Laboratory, Beijing 100190, China)
③(University of Chinese Academy of Sciences, Beijing 100190, China)
④(Academy of Opto-Electronics, Chinese Academy of Sciences, Beijing 100094, China)

The model-based target decomposition is one of the most fundamental applications to the polarimetry SAR. However, it can only retrieve the power information of three scattering mechanisms with polarimetric SAR data. By applying this model-based decomposition technique to the π/4mode Compact Polarimetric SAR Interferometry (C-PolInSAR) data and decomposing the cross correlation matrix, this study can obtain the power contributions of three scattering mechanisms as well as their scattering phase center. Firstly, this study introduces the models corresponding to the three scattering mechanisms under π/4mode C-PolInSAR observation. Then, the targets can be decomposed by the numerical method, and finally both power contributions and scattering center of mechanisms are retrieved. The simulation data validate the effectiveness of the decomposition algorithm and the impacts of different wave-band and different ground parameters on the decomposition results are analyzed.

SAR; π/4mode compact polarimetry; Compact Polarimetric SAR Interferometry (C-PolInSAR); Model-based decomposition

TN958

A

1009-5896(2015)01-0022-07

10.11999/JEIT140550

2014-04-28收到，2014-09-18改回

國家自然科學(xué)基金(61431018, 41101335)資助課題

*通信作者：洪文 wendy_iecas@163.com