李芳芳 胡東輝 丁赤飚 仇曉蘭
(中國科學院空間信息處理與應用系統技術重點實驗室 北京 100190)
(中國科學院電子學研究所 北京 100190)
機載雙天線InSAR對飛數據處理與分析
李芳芳*胡東輝丁赤飚仇曉蘭
(中國科學院空間信息處理與應用系統技術重點實驗室北京 100190)
(中國科學院電子學研究所北京100190)
干涉合成孔徑雷達(InSAR)技術具有高精度的地形測繪能力。然而在山區地形條件下,受SAR側視成像的影響,存在較多的幾何畸變區域,干涉相位表現為不連續或者缺少有效信息的情況,對單一的干涉圖像對處理難以得到精確的數字高程模型(DEM)。融合兩個或多角度的干涉數據可以用于解決這一問題。該文根據這一思路利用機載雙天線InSAR對飛數據進行實驗,針對山區地形條件下參考高程近似影響運動補償精度的問題提出了基于高程迭代的運動補償方法,對于陰影、疊掩區域容易導致相位解纏誤差的問題提出基于地形特征的相位解纏方法,從而盡可能降低單一角度數據的DEM誤差,以達到消除對飛數據重疊區域3維定位不一致的目的,拼接實驗結果驗證了處理方法的有效性。
機載InSAR;運動補償;相位解纏;對飛數據融合
干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar,InSAR)是利用同一地區多幅SAR復圖像獲取相位信息,進行地形測繪和地表形變監測的一項技術[1]。然而,由于SAR側視成像的特點,在地形起伏較大的山區,疊掩、陰影所占比例較大,獲取的干涉相位會出現不連續甚至全是噪聲的情況,無法利用常規的InSAR技術重建精確的數字高程模型(Digital Elevation Model,
DEM)[2,3]。解決這一問題的一種思路是融合同一區域的多角度觀測數據,從而在某一角度下的幾何畸變區域可以利用其它角度的數據補償[4-6]。
然而,在實際數據處理過程中,如果得到的單一角度DEM誤差較大,則在多角度數據拼接時會出現“拼接縫”的現象。而在山區地形條件下,人為地布設控制點非常困難,無法通過不同角度圖像中控制點之間的仿射變換消除目標3維定位的不一致性。因此,要達到在無控制點情況下的多角度數據“無縫拼接”,需要在單一角度數據處理的過程中盡可能地減小誤差,從而盡可能消除重疊區域DEM的不一致現象。
一方面,運動補償是機載干涉SAR處理的關鍵步驟之一。在山區地形條件下,參考高程的誤差是引起運動補償誤差的主要原因,從而影響干涉相位和反演DEM的精度。另一方面,受幾何畸變的影響,山區地形的相位解纏也是一個復雜的問題。本文利用機載雙天線干涉SAR對飛數據進行融合處理實驗,以實現山區地形DEM重建,重點針對上述兩個關鍵問題展開研究,分別提出基于DEM迭代的運動補償方法和基于地形特征的相位解纏方法。實驗結果驗證了處理方法的有效性。
機載雙天線干涉SAR對飛數據拼接實現山區地形干涉測量的流程如圖1所示。本節將重點介紹運動補償和相位解纏兩個關鍵步驟的處理方法,以提高各個角度的DEM精度,進而實現對飛數據的“無縫拼接”。
2.1運動補償
高精度的運動補償是高分辨率機載SAR成像的關鍵。機載干涉SAR系統對運動補償的精度要求更為嚴苛。對于山區地形干涉處理來說,運動補償存在誤差時,一方面引起單一角度數據的干涉相位誤差,從而使DEM測量精度下降,另一方面在多角度數據拼接時導致重疊區域的DEM不一致,出現“拼接縫”現象。

圖1 對飛數據拼接實現山區地形測繪流程圖Fig.1 Flow chart of processing two anti-parallel aspects of InSAR data
傳統的運動補償算法通常存在參考高程的近似[7,8],即利用平面高程作為該區域所有目標的高程,在山區等起伏較大的地形條件下,該假設與實際差別較大,會導致較大的運動補償殘余誤差,進而影響干涉相位和DEM的精度。要補償高程不準引起的誤差,需要使參考DEM盡可能接近真實值。因此,本節提出通過高程迭代的方法減小高程誤差引起的運動補償殘余誤差。
補償流程如圖2所示,首先按照平地假設進行運動補償及成像處理,經過雙天線圖像配準、相位濾波、相位解纏等干涉處理步驟后,根據幾何關系反演生成粗精度的DEM;然后將該DEM作為新的參考高程,重新進行運動補償。
下面通過定量計算分析高程迭代的精度。如圖3為存在高程誤差時的運動補償幾何示意圖,設天線1的運動誤差為[Δx(t),0,Δz(t )],天線2的運動誤差為[Δx(t)+δx(t),0,Δz(t)+δz(t )]。文獻[9]給出了存在高程誤差時,兩天線運動補償殘余誤差的不一致分量為:

其中,B,α分別為基線長度和基線角,φ1,φ2分別為天線1和天線2中目標A的視角,Δφ=φ2-φ1,α(t)如圖4所示,Δhinit為初次運動補償時的參考高程誤差,由此得到高程測量誤差為:

假設干涉系統的參數如表1,運動誤差為10 m,橫滾角變化為θr=3°,則由其引起的兩天線不一致運動誤差分別為:δx(t)=B [cos(α+θr)-cosα],δz(t)=B[ sin(α+θr)-sinα],由此可以計算出A≈0.05?1。由此可見,盡管該DEM并不準確,但是與平地假設相比,已大大減小了高程誤差,在高程迭代時不會收斂到相反方向。通常通過一到兩次迭代即可滿足精度要求。

圖2 高程迭代的運動補償流程圖Fig.2 Flow chart of MOCO based on iterative DEM

表1 干涉系統參數Tab.1 InSAR system parameters

圖3 存在高程誤差的運動補償幾何示意圖Fig.3 MOCO geometry with reference height error

圖4 斜距平面運動誤差示意圖Fig.4 Motion error in the slant range plane
2.2相位解纏
2.2.1陰影區域山區地形存在較大面積的陰影區域,干涉相位呈噪聲特性,殘差點密集,相位解纏時容易引起相位誤差的傳播[10,11]。本小節在分析陰影區域干涉相位特性的基礎上,利用相位補償的方法避免了相位誤差的傳播。
干涉相位沿距離向的局部條紋頻率可表示為:

在陰影起始點和結束點之間的地形坡度角為?=θ-90°,代入式(3)即有Δf≈0,表明陰影區域沿距離向的局部條紋頻率接近為零,由此可知兩個臨界端點的干涉相位之差近似為零。為了進一步確定實際系統陰影區域相位變化的范圍,下面以中科院電子所X波段機載雙天線干涉系統為例,仿真計算該相位差,在不同視角和山體高度下的相位差如圖5所示,可以看出該相位差很小,絕對值均在π以內,也就是說陰影區域兩個臨界端點之間不存在相位周期的改變。這一特征可以用來確定陰影區域兩側干涉相位的相對值,從而可以避免相位解纏路徑穿過陰影區域時引起的相位誤差的傳播。
根據上述分析,我們利用陰影區域起始點和結束點處的干涉相位線性擬合出陰影區域內部的干涉相位,不會影響相位解纏時相位模糊數的計算。具體的相位解纏步驟如下:
(1)根據相干系數圖和SAR幅度圖像檢測并提取陰影區域,并利用形態學方法使區域邊緣規整化;
(2)在每一方位向陰影區域兩側選取相位擬合的起始點和結束點,干涉相位分別為φtop,φbot,這里取距離陰影邊界最近的相干系數較高的點,以保證所選起始點和結束點相位的可靠性;
(3)線性擬合出陰影區域的干涉相位,即

其中,N為陰影區域沿距離向的寬度,d為待補償位置距離起始點的像素單元;
(4)對相位補償后的干涉相位圖,利用基于質量指導的路徑跟蹤法進行解纏,得到解纏相位。

圖5 不同高度和視角對應陰影起始和結束點干涉相位差Fig.5 Interferometric difference between the beginning and end of shadow with different height and look angle
按照上述的相位解纏步驟對實測的包含陰影區域的干涉相位進行相位補償并解纏,驗證上述方法的正確性。所用數據為中科院電子所機載X波段雙天線InSAR系統在綿陽地區的實驗數據。實驗數據及結果如圖6所示。其中圖6(a)為SAR幅度圖像,圖6(b)為濾波后的干涉相位圖,圖6(c)為相干系數圖,可以看出圖中有較大面積的陰影區域,在幅度圖像中表現為較暗的區域,干涉相位中表現為噪聲,相干系數很低。


圖6 陰影區域相位解纏結果Fig.6 Phase unwrapping result of shadow area
圖6(d)為檢測出的陰影區域的掩膜圖,掩膜后進行直接相位解纏可以避免陰影區域相位噪聲引起的誤差傳播,然而從掩膜圖可以看到掩膜后使干涉相位左右兩側形成不連通區域,因此直接解纏得到圖6(e)的結果,陰影區域兩側解纏相位有明顯跳變,相位的相對值與理想值不一致。根據本節方法對陰影區域進行相位補償,實際上并不需要對所有陰影區域均做補償,只要使補償后高質量區域連通即可。因此對陰影區域部分進行相位補償后的結果如圖6(f),解纏后的結果如圖6(g),此時補償相位將兩側區域連通,使得解纏相位連續。圖6(h)給出了利用直接解纏和相位補償后解纏得到的DEM在圖6(g)黑色直線所示的方位向上的對比,可見不進行相位補償時,DEM會有較大跳變,而相位補償后的DEM基本上是連續的,與實際情況相符。
2.2.2疊掩區域疊掩是SAR圖像中的另一種幾何畸變現象,由于SAR側視成像的特點,當地面場景的迎坡坡度過大時,會導致不同高度區域的回波投影到同一個距離-多普勒單元內,形成疊掩現象。疊掩區域的干涉相位反映的是多個散射源矢量疊加的結果,干涉相位存在跳變,解纏時無法傳遞正確的相位纏繞周期。多基線技術利用信號譜估計等方法可以實現對疊掩區域多個散射源的分辨[12-14],然而目前多基線干涉數據的缺乏限制了該研究的開展。本小節通過分析疊掩區域的干涉相位特性,提出單基線條件下基于頻率估計的疊掩區域相位解纏方法。
圖7給出了簡單山體形成疊掩的SAR成像幾何關系示意圖??梢钥闯?,山頂點P1到SAR的斜距為rP1,小于山底點P2到SAR的斜距rP2,因此在斜距rP1和rP2之間的每一個距離向像素單元內均由來自區域A,區域B和區域C中的3個散射源疊加形成,從而產生了疊掩區域。而在斜距小于rP1的區域D和斜距大于rP2的區域E中,回波信號都僅包含單個散射源,是非疊掩區域。
通過仿真計算出圖7場景沿距離向的干涉相位如圖8所示,圖中兩條豎直黑色虛線之間表示了疊掩區域的干涉相位,這里分別計算了形成疊掩的3部分區域的干涉相位??梢钥闯?,區域B的相位頻率與區域A和區域C相反。而區域A的相位與區域D的相位是連續的,同樣,區域C的相位與區域E的相位連續。此外,在疊掩區域的起始點P1處,區域B和區域C重合為一點,因此具有相同的干涉相位φ4;與之類似,在疊掩區域的結束點P2處,區域B和區域A具有相同的干涉相位φ2。通常情況下認為疊掩區域中迎坡面的散射強度占主導地位,因此疊掩區域的干涉相位通常表現為區域B的相位特征,即具有與非疊掩區域相反的相位頻率,且在起始點和結束點處存在相位的跳變。在這一假設前提下,圖8中的絕對相位纏繞后的相位如圖9所示。觀察實際的干涉數據可以發現疊掩區域確實具有與相鄰區域相反的相位頻率,圖10給出了一幅四川綿陽地區含疊掩區域的實測干涉相位圖,黑色方框內部的區域為疊掩區域。

圖7 SAR疊掩區域成像幾何關系示意圖Fig.7 SAR imaging geometry of layover area

圖8 不同區域的絕對干涉相位Fig.8 Absolute phase of different areas

圖9 纏繞干涉相位Fig.9 Wrapped phase under the assumption
根據式(3)的干涉相位沿距離向的局部頻率,可以從理論上解釋上述仿真結果。當迎坡面的坡度角大于雷達視角時,在SAR圖像中產生疊掩現象。也就是對于圖7所示的區域B而言,有?B>θ,由式(3)可知,此時區域B對應干涉相位的距離向局部頻率ΔfB<0。而對于與區域B中的散射源投影到相同位置的區域A和區域C來說,區域A的坡度角0°<?A<θ,區域C的坡度角θ-90°<?C<0°,代入式(3)分別有ΔfA>0,ΔfC>0。與疊掩區域相鄰的區域D和區域E分別具有與區域A和區域C相同的坡度角,因此同樣有ΔfD>0,ΔfE>0。當迎坡面散射強度占主導地位時,疊掩區域的干涉相位在距離向呈現出與其它區域相反的頻率特性。這就從理論上解釋了以上仿真得出的結論。

圖10 實測疊掩區域干涉相位圖Fig.10 Real Interferogram with layover areas
根據上述分析,利用圖8所示的區域A和區域C分別與區域D及區域E之間的相位連續性,可以通過估計區域D和區域E的干涉相位頻率,從而獲得區域A和區域C的相位頻率。這樣在解纏得到區域D的相位φ1或區域E的相位φ3后,就可以進一步估計出區域B邊界的解纏相位φ2或φ4,從而消除相位跳變的影響。根據這一思想,本文提出了基于局部頻率估計的疊掩區域相位解纏方法,具體實現步驟如下:
(1)首先利用相干系數圖、干涉相位圖以及SAR幅度圖像檢測并提取出疊掩區域;
(2)將疊掩區域掩膜后,對其它區域利用傳統方法進行相位解纏,此時得到如圖8中區域D和區域E的相位分別為φDr和φEr,其中,rP1和rP2處的相位分別為φ1r和φ3r,此時φ1r和φ3r之間的相對值與實際情況存在偏差(相位φ的下標r表示relative,即相對相位,下同);
(3)將非疊掩區域掩膜,僅對疊掩區域利用傳統方法進行相位解纏,得到如圖8中區域B的解纏相位,其邊界值分別為φ2r和φ4r,此時的解纏相位并不是反映地形高程的絕對相位,而與絕對相位之間相差一個常數值;
(4)通過局部頻率估計[15,16]得到疊掩附近區域,即如圖8中的區域D及區域E沿距離向的干涉相位頻率分別為fr_D和fr_E;
(5)假設區域D的解纏相位即為絕對相位,即φDa=φDr,φ1a=φ1r(相位φ的下標a表示absolute,即絕對相位,下同),則根據區域A和區域D之間的相位連續性,可以得出區域B在rP2處的絕對相位為:

其中,N為疊掩區域沿距離向的像素數。進一步,區域B在rP1處以及其內部的絕對相位可分別表示為:

再利用區域C和區域E之間的相位連續性,則可以得出區域E在rP2處的絕對相位為:

因此,區域E的絕對相位可表示為:

通過以上步驟,整個區域的相位得到展開,而且可以消除疊掩區域邊界相位不連續的影響,避免相位誤差的傳播。
按照上述相位解纏方法對圖10中所示的實測數據進行干涉解纏,從而驗證本文方法的有效性。

圖11 含疊掩區域實測數據Fig.11 Real InSAR data including layover area
圖11(a)和圖11(b)給出了該實測數據對應的SAR幅度圖像和相干系數圖,由此檢測出疊掩區域的掩膜圖如圖11(c)。分別用本文方法和簡單的質量指導解纏方法對圖10的干涉相位進行解纏,得到的解纏相位分別如圖12(a)和圖12(b)所示。圖12(c)對比了圖12(a)和圖12(b)中黑色直線所示方位門的兩種方法的解纏相位,圖12(d)給出了與之相對應的DEM解算結果對比。可以看出,由于迎坡面的相位占主導地位,因而本文方法的解纏相位在疊掩起始點和結束點均出現了相位的跳變,相應的DEM也存在跳變而且在疊掩區域呈遞減趨勢,也就是疊掩表現出的“頂底倒置”現象,這與前面的理論分析是一致的。而直接用質量指導法解纏由于未考慮疊掩區域邊界相位不連續的影響,造成疊掩區域解纏相位的誤差,在積分路徑穿過疊掩區域時進而將誤差傳播到疊掩右側區域。

圖12 實測疊掩區域相位解纏實驗結果Fig.12 Phase unwrapping result of layover area
進一步,將本文方法解算出的DEM投影到地距坐標系中,并與外源粗DEM數據ASTER DEM進行比較如圖12(e)所示,可以看出二者基本相符。由于本文方法僅反演出了疊掩區域中迎坡面的高程信息,因而在地距坐標系中存在無高程信息的區域,即相當于圖7中示意的區域A和區域C的高程沒有重建。要得到疊掩區域多個散射源的高程信息,可以利用多基線干涉技術進行分辨不同的散射源,這里由于缺乏多基線數據,因而不再進行討論。
本文采用的實驗數據是由中科院電子所研制的X波段機載雙天線干涉SAR系統于2011年7月在四川綿陽地區采集的數據,實驗通過從兩個方向對同一場景進行照射,獲得兩組干涉數據,其中一組數據為由東向西飛行,另一組數據為由西向東飛行,航跡示意圖如圖13所示。兩個角度得到的圖像差異較大,尤其在山區地形條件下,將會形成不同的幾何畸變區域,利用這一特點通過拼接互補實現幾何畸變區域的干涉測量。系統參數如表2所示。
本節利用上述方法對一幅7 km×10 km大小的圖像分別進行兩個角度的處理得到各自的正射DEM和DOM,并進行拼接,結果如圖14所示。可以看出,在拼接圖中沒有明顯的拼接縫現象,達到了較好的拼接效果,但仍然存在少量無法測量高程的區域(圖中高程置零的藍色區域),這是由于這些區域在兩個角度的圖像中均表現為相干性較差的區域,后續還可以進一步利用多個角度數據進行補充。

圖13 實驗數據航跡示意圖Fig.13 Anti-parallel flight trajectories

表2 系統參數Tab.2 Specifications of CAS InSAR system
山區地形測繪是干涉SAR處理的一個難點,本文利用機載雙天線干涉SAR對飛數據融合拼接實現山區DEM重建。針對運動補償和相位解纏兩個關鍵步驟,分別給出了基于高程迭代的運動補償方法和基于相位補償的陰影區域相位解纏方法,從而獲得精確的單一角度DEM。實驗結果驗證了對飛數據融合的有效性。


圖14 實驗結果Fig.14 Experimental results
[1]Rosen P A,Hensley S,Joughin I R,et al..Synthetic aperture radar interferometry[J].Proceedings of the IEEE,2000,88(3):331-382.
[2]Eineder M and Holzner J.Interferometric DEMs in Alpine terrain-limits and options for ERS and SRTM[C].Proceedings of International Geoscience and Remote Sensing Symposium,Honolulu,Hawaii USA,2000:3210-3212.
[3]Eineder M.Interferometric DEM reconstruction of Alpine areas-experiences with SRTM data and improved strategies for future missions[OL].www.eproceedings.org,2005.
[4]Schmitt M and Stilla U.Utilization of airborne multi-aspect InSAR data for the generation of urban ortho-images[C].International Geoscience and Remote Sensing Symposium,2010:3937-3940.
[5]Lachaise M,Balss U,Fritz T,et al..The dual-baseline interferometric processing chain for the TanDEM-X mission[C].Proceedings of International Geoscience and Remote Sensing Symposium,Munich,Germany,2012:5562-5565.
[6]Gruber A,Wessel B,Huber M,et al..The approach for combing DEM acquisitions for the TanDEM-X DEM mosaic[C].Proceedings of International Geoscience and Remote Sensing Symposium,Melbourne,Australia,2013:2970-2973.
[7]Fornaro G.Trajectory deviations in airborne SAR:analysis and compensation[J].IEEE Transactions on Aerospace and Electronic Systems,1999,35(3):997-1009.
[8]Reigber A,Alivizatios E,Potsis A,et al..Extended wavenumber-domain synthetic aperture radar focusing with integrated motion compensation[J].IEE Proceedings-Radar,Sonar and Navigation,2006,153(3):301-310.
[9]李芳芳,仇曉蘭,孟大地,等.機載雙天線InSAR運動補償誤差的影響分析[J].電子與信息學報,2013,35(3):559-567.
Li F F,Qiu X L,Meng D D,et al..Effects of motion compensation errors on performance of airborne dual-antenna InSAR[J].Journal of Electronics & Information Technology,2013,35(3):559-567.
[10]王健,向茂生,李紹恩.一種基于InSAR相干系數的SAR陰影提取方法[J].武漢大學學報(信息科學版),2005,30(12):1063-1066.
Wang J,Xiang M S,and Li S E.A method for extracting the SAR shadow form InSAR coherence[J].Geomatics and Information Science of Wuhan University,2005,30(12):1063-1066.
[11]索志勇,李真芳,吳建新,等.干涉SAR陰影提取及相位補償方法[J].數據采集與處理,2009,24(3):264-269.
Suo Z Y,Li Z F,Wu J X,et al..Method for shadow detection and compensation on interferometric synthetic aperture radar[J].Journal of Data Acquisition & Processing,2009,24(3):264-269.
[12]Gini F,Lombardini F,and Montanari M.Layover solution in multibaseline SAR interferometry[J].IEEE Transactions on Aerospace and Electronic Systems,2002,38(4):1344-1356.
[13]Baselice F,Budillon A,Ferraioli G,et al..Joint SAR imaging and DEM reconstruction from multichannel layover-affected SAR data[C].Proceedings of International Geoscience and Remote Sensing Symposium,Cape Town,South Africa,2009:25-28.
[14]Baseline F,Ferraioli G,and Pascazio V.DEM reconstruction in layover areas from SAR and auxiliary input data[J].IEEE Geoscience and Remote Sensing Letters,2009,6(2):253-257.
[15]Li Hai,Li Zhen-fang,Liao Gui-sheng,et al..An estimation method for InSAR interferometric phase combined with image auto-coregistration[J].Science in China,Series F:Information Sciences,2006,49(3):386-396.
[16]Trouvé E,Nicolas J M,andMa?tre H.Improving phase unwrapping techniques by the use of local frequency estimates[J].IEEE Transactions on Geoscience and Remote Sensing,1998,36(6):1963-1972.

李芳芳(1986-),女,山西人,助理研究員,博士,主要研究領域為干涉SAR信號處理。
E-mail:liff86@gmail.com
胡東輝(1970-),男,北京人,副研究員,碩士,主要研究領域為SAR信號處理、SAR定標等。
丁赤飚(1969-),男,陜西人,研究員,博士,主要研究領域為合成孔徑雷達、遙感信息處理及應用系統等。
仇曉蘭(1982-),女,江蘇人,副研究員,博士,主要研究領域為雙站SAR、SAR目標特性等。
Antiparallel Aspects of Airborne Dual-antenna InSAR Data Processing and Analysis
Li Fang-fangHu Dong-huiDing Chi-biaoQiu Xiao-lan
(Key Laboratory of Technology in Geo-spatial Information Processing and Application System,Institute of Electronics,
Chinese Academy of Sciences,Beijing 100190,China)
(Institute of Electronics,Chinese Academy of Sciences,Beijing 100190,China)
Interferometric Synthetic Aperture Radar(InSAR)is a powerful technique for precise topographic mapping.However,owing to the side-looking SAR imaging geometry,geometry distortions appear in mountainous scenarios.Because of phase discontinuities or the absence of a valid phase,it is difficult to recover accurate DEM in such areas with single-aspect InSAR data.Fusion of two or more different aspects of InSAR data can deal with this problem in practice.Experiments using two antiparallel aspects of airborne InSAR data are carried out based on this idea.To decrease the processing error in single-aspect data and fuse them seamlessly,a MOtion COmpensation(MOCO)method using iterative DEM is used to reduce the MOCO error.Besides,phase-unwrapping methods based on terrain characteristics are proposed to avoid phase-unwrapping error owing to phase discontinuities in areas of shadow and layover.Experimental results verify the effectiveness of the processing methods.
Airborne Interferometric SAR(InSAR); MOtion Compensation(MOCO); Phase unwrapping; Fusion of anti-parallel data
TP391
A
2095-283X(2015)01-0038-11
10.12000/JR14135
李芳芳,胡東輝,丁赤飚,等.機載雙天線InSAR對飛數據處理與分析[J].雷達學報,2015,4(1):38-48.http://dx.doi.org/10.12000/JR14135.
Reference format:Li Fang-fang,Hu Dong-hui,Ding Chi-biao,et al..Antiparallel aspects of airborne dual-antenna InSAR data processing and analysis[J].Journal of Radars,2015,4(1):38-48.http://dx.doi.org/10.12000/JR14135.
2014-11-20收到;2014-12-18改回;2015-01-30網絡優先出版
國家863計劃(2007AA120302)和國家自然科學基金(61401428)資助課題
李芳芳liff86@gmail.com