何林幫,趙建虎,張紅梅,王 曉,嚴(yán) 俊
(1.武漢大學(xué)測繪學(xué)院,湖北武漢430079;2.武漢大學(xué)動力與機(jī)械學(xué)院自動化系,湖北武漢430072)
顧及姿態(tài)角的多波束聲線精確跟蹤方法
何林幫1,趙建虎1,張紅梅2,王 曉1,嚴(yán) 俊1
(1.武漢大學(xué)測繪學(xué)院,湖北武漢430079;2.武漢大學(xué)動力與機(jī)械學(xué)院自動化系,湖北武漢430072)
傳統(tǒng)聲線跟蹤方法中波束初始入射角的計(jì)算采用理想狀態(tài)下多波束提供的波束分配角,未顧及或未完全顧及船體姿態(tài)角的影響進(jìn)而導(dǎo)致波束海底投射點(diǎn)歸位計(jì)算帶入顯著性誤差問題。針對此問提出了一種顧及姿態(tài)角影響的聲線精確跟蹤方法,該方法通過研究船體姿態(tài)影響下的多波束換能器狀態(tài),給出了波束聲線實(shí)際傳播面,并推導(dǎo)出了傳播面內(nèi)波束聲線的實(shí)際初始入射角,結(jié)合聲速剖面,最終根據(jù)常梯度聲線跟蹤得到了準(zhǔn)確的測深點(diǎn)三維坐標(biāo)。通過實(shí)驗(yàn)比較傳統(tǒng)聲線跟蹤方法和精確聲線跟蹤方法計(jì)算波束在海底投射點(diǎn)的坐標(biāo)精度,驗(yàn)證了精確聲線跟蹤方法克服了傳統(tǒng)聲線跟蹤的不足,顯著地提高了多波束測深的精度。
多波束測深系統(tǒng);波束腳??;初始入射角;姿態(tài)角;常梯度聲線跟蹤
隨著海洋資源開發(fā)的不斷推進(jìn),海底勘探作為海洋資源開發(fā)的基礎(chǔ)性工作,且測深工作在海底油氣和礦產(chǎn)資源開發(fā)、海底電纜鋪設(shè)、河道清淤、工程建設(shè)方案選取等工程應(yīng)用方面發(fā)揮著重要的作用,測深精度也直接影響著工程應(yīng)用的效果。在多波束測深中,每個(gè)波束在海底投射點(diǎn)的坐標(biāo)均基于一定的聲線跟蹤模型獲得的,聲線跟蹤精度直接影響著多波束的測深精度[1]。實(shí)際數(shù)據(jù)處理發(fā)現(xiàn),聲線跟蹤常會帶來一些測深假象,并導(dǎo)致海底地形異常[2]。根據(jù)聲線跟蹤原理,影響聲線跟蹤精度的主要因素有聲速剖面誤差、波束初始入射角誤差和聲線跟蹤模型誤差。聲速剖面誤差是影響測深精度的第一誤差源,包括聲速測量誤差和聲速代表性誤差,可以借助高精度聲速剖面儀和施測一定密度的聲速剖面來消除[3?4];基于層內(nèi)常梯度的聲線跟蹤方法,由于假設(shè)與實(shí)際一致,具有較高的聲線跟蹤精度[4]。實(shí)際應(yīng)用中發(fā)現(xiàn),消除了上述2個(gè)因素影響后,測深結(jié)果仍不理想,尤其是在深水區(qū)多波束測量中[3?5]。分析認(rèn)為,主要由于目前聲線跟蹤方法中的波束初始入射角計(jì)算未顧及或未完全顧及船姿因素所致。目前,聲線跟蹤中波束初始入射角的獲取有2種方式:1)直接來自換能器提供的各波束分配角[4];2)僅顧及了橫搖角,即波束初始入射角為波束分配角與橫搖角之和[3]。以上2種初始入射角獲取方法均基于一個(gè)假設(shè),即多波束Ping測深斷面與測量船航跡方向正交,而實(shí)際多波束換能器在橫搖、縱搖姿態(tài)角作用下,Ping測深斷面與航跡方向已不再正交,而前者未顧及這一影響,后者雖顧及了橫搖影響,但忽視了縱搖;此外,在對姿態(tài)影響的處理方面,2種方法均先通過聲線跟蹤獲得波束點(diǎn)坐標(biāo),爾后再借助姿態(tài)角構(gòu)建旋轉(zhuǎn)矩陣進(jìn)行強(qiáng)制旋轉(zhuǎn)變換,獲得波束海底測深點(diǎn)在理想船體坐標(biāo)系下的坐標(biāo)。不同的是,前者借助縱搖、橫搖角構(gòu)建的旋轉(zhuǎn)矩陣來實(shí)施變換,而后者僅借助縱搖角構(gòu)建的旋轉(zhuǎn)矩陣進(jìn)行變化。綜上可知,現(xiàn)有波束測深點(diǎn)位計(jì)算均與實(shí)際存在著一定的誤差,也因此導(dǎo)致了前述測深假象,且隨著水深的增加,其影響愈加顯著。因此,本文提出了一種精確的聲線跟蹤方法,以期通過研究姿態(tài)對Ping測深斷面的影響,獲得聲線在空間的實(shí)際傳播路徑和顧及姿態(tài)角的實(shí)際波束初始入射角,從而實(shí)現(xiàn)波束測深點(diǎn)位置的正確計(jì)算,提高多波束測深精度。
1.1 常梯度聲線跟蹤方法及波束初始入射角的影響分析
海水中聲速與海水溫度、鹽度和靜壓力相關(guān),隨著深度的變化而變化[6]。由于很難獲得聲速隨深度變化的函數(shù),通常只能借助聲速剖面儀獲得一定深度間隔的聲速剖面,并根據(jù)聲速變化,借助Snell法則,沿聲線跟蹤,獲得波束海底投射點(diǎn)在船體坐標(biāo)系下的坐標(biāo)[7?8]。在每個(gè)水層中,由于只知道該層上界和下界聲速,因此通常假設(shè)聲速在該層以常梯度g傳播,采用類似處理方法處理其他各水層,并沿著聲線跟蹤到海底,即實(shí)現(xiàn)常梯度聲線跟蹤和波束海底坐標(biāo)計(jì)算[9]。
假設(shè)聲線從換能器發(fā)射,經(jīng)過N個(gè)水層,每個(gè)水層上界Zi聲速為Ci,下界Zi+1聲速為Ci+1,層內(nèi)聲速以常梯度gi變化,則層i內(nèi)聲速函數(shù)為
如圖1聲線跟蹤示意圖所示,在常梯度gi聲速變化下,聲線在第i層內(nèi)的傳播軌跡為一連續(xù)的、帶有一定曲率、半徑為Ri的弧段[4],若Snell常數(shù)為p,則Ri為
圖1 常梯度聲線跟蹤示意圖Fig.1 Schematic diagram of constant?gradient sound ray tracking
在分層聲線跟蹤時(shí),除了計(jì)算整層的垂直位移、水平位移和傳播時(shí)間外,還需要依據(jù)傳播剩余時(shí)間計(jì)算剩余層的垂直位移和水平位移。假設(shè)聲線在第i層內(nèi)傳播時(shí),聲線在該層內(nèi)r點(diǎn)處(如圖1)結(jié)束,此時(shí)剩余時(shí)間tr等于波束單程旅行時(shí)間tall減去第i層以前累計(jì)的傳播時(shí)間,則聲線在剩余層內(nèi)的垂直位移Δzr和水平位移Δyr為
則聲線傳播總的垂直位移z和水平位移y為
姿態(tài)角(橫搖角r和縱搖角p)改變了換能器的理想發(fā)射狀態(tài),為顧及該影響,傳統(tǒng)數(shù)據(jù)處理方法借助換能器的波束分配角完成上述聲線跟蹤后得到波束坐標(biāo)(x,y,z),且根據(jù)橫搖角r和縱搖角p波束坐標(biāo)的影響進(jìn)行了坐標(biāo)旋轉(zhuǎn)變換,從而獲得了波束在船體坐標(biāo)系下的坐標(biāo)(x',y',z')[10]。
以上傳統(tǒng)數(shù)據(jù)處理方法存在如下問題:
1)認(rèn)為聲線傳播面,也即Ping斷面,與航向方向正交,而實(shí)際聲線傳播面并非如此。
2)受姿態(tài)因素影響,換能器的理想狀態(tài)已發(fā)生變化,實(shí)際波束初始入射角不再是理想狀態(tài)下的波束分配角,而是受姿態(tài)角影響的入射角。根據(jù)Snell法則sin θ0/C0=sin θ1/C1=p,若初始入射角為θ0,其誤差為dθ0,則有[11]
式(8)表明,初始入射角存在誤差dθ0時(shí),會給下層入射角θ1產(chǎn)生dθ1的誤差。由于采取分層跟蹤處理,也會進(jìn)一步影響后續(xù)層的折射角θi,并給最終測深點(diǎn)坐標(biāo)計(jì)算帶來系統(tǒng)性誤差。
3)上述對波束聲線先跟蹤后旋轉(zhuǎn)的處理方法,實(shí)際是將理想狀態(tài)下的波束測深點(diǎn)旋轉(zhuǎn)到實(shí)際的過程,由于1)、2)兩方面的影響,這種旋轉(zhuǎn)變換只能給測深點(diǎn)位坐標(biāo)計(jì)算帶來新的誤差。
以上問題表明,在聲速誤差和聲線跟蹤模型正確的情況下,尋找顧及姿態(tài)影響的實(shí)際聲線傳播面及正確的初始入射角是實(shí)現(xiàn)波束點(diǎn)位坐標(biāo)精確計(jì)算的關(guān)鍵。為此,下面研究給出一種顧及姿態(tài)影響下的波束初始入射角的精確計(jì)算方法,從而實(shí)現(xiàn)測深點(diǎn)坐標(biāo)精確跟蹤計(jì)算。
1.2 顧及姿態(tài)的波束初始入射角計(jì)算
由于船體坐標(biāo)系的中心通常是以換能器為中心,因此以換能器理想水平狀態(tài)作為基準(zhǔn)面來分析船體姿態(tài)對波束初始入射角的影響,如圖2所示的換能器基陣坐標(biāo)系中,水平狀態(tài)的換能器基準(zhǔn)面位于OABC平面內(nèi),O為換能器中心,OA為基準(zhǔn)面縱軸正方向,OC為基準(zhǔn)面橫軸正方向。設(shè)OA長度為a,OC長度為c,A、B兩點(diǎn)的坐標(biāo)分別為(a,0,0)和(0,c,0)。在某一姿態(tài)(橫搖、縱搖角分別為r和p)影響下基準(zhǔn)面變化為OA1B1C1,即基陣面由水平先繞OX軸旋轉(zhuǎn)角度α(α≠r),再繞OY軸旋轉(zhuǎn)角度β形成。A點(diǎn)、C點(diǎn)經(jīng)過兩次旋轉(zhuǎn)后分別轉(zhuǎn)到A1和C1位置,A1、C1兩點(diǎn)在水平面OXY上的投影分別為A2和C2。在此狀態(tài)下,OA1與水平面夾角∠A1OA2即為縱搖角p,OC1與水平面夾角∠C1OC2即為橫搖角r。根據(jù)縱搖角、橫搖角和旋轉(zhuǎn)角定義,r和α符號一致,p和β符號一致。
圖2 換能器基陣旋轉(zhuǎn)角度模型Fig.2 Transducer array rotation model
由以上過程可知,基準(zhǔn)面OABC經(jīng)過α和β兩次旋轉(zhuǎn)得到OA1B1C1,則
則旋轉(zhuǎn)后A1點(diǎn)坐標(biāo)為
由式(10)得到旋轉(zhuǎn)后A1點(diǎn)坐標(biāo)再根據(jù)三角形正弦定理可計(jì)算出基陣面的縱搖角p(即∠A1OA2):
式中:ZA1為A1點(diǎn)在Z軸上的坐標(biāo),根據(jù)β與p符號一致可得
類似地,由式(11)得到旋轉(zhuǎn)后C1點(diǎn)坐標(biāo)再根據(jù)三角形正弦定理可計(jì)算出基陣面的橫搖角r(即∠C1OC2):
式中:ZC1為C1點(diǎn)在Z軸上的坐標(biāo),r和α符號一致,并將β=p代入式(14)得
由式(13)和(16)可知,在旋轉(zhuǎn)變換中,繞OY軸旋轉(zhuǎn)角β等于縱搖角p,而繞OX軸旋轉(zhuǎn)角α并不等于橫搖角r。因此,在目前較為精細(xì)的聲線跟蹤計(jì)算中,即使顧及姿態(tài)影響,定義初始入射角為θ0+r顯然是不正確的。
為了得到姿態(tài)(r、p)影響下真實(shí)的波束入射角,下面推導(dǎo)實(shí)際波束初始入射角θ0'的計(jì)算模型。
由以上推導(dǎo)可知,實(shí)際聲線可由理想狀態(tài)下的聲線經(jīng)α、β旋轉(zhuǎn)變換r后得到。設(shè)理想狀態(tài)下,第i個(gè)波束分配初始入射角為θi,在不失精度情況下,假設(shè)經(jīng)歷第一個(gè)水層以常聲速傳播,傳播距離為Ri,則波束在第一水層下界的落點(diǎn)Pi坐標(biāo)為(0,Risin θi,Ricos θi),而在姿態(tài)影響下的實(shí)際坐標(biāo)(xi,yi,zi)為
式(17)可由圖3解釋。假設(shè)換能器基陣水平時(shí),第i號波束的波束角度為θi,斜距為R,則點(diǎn)A坐標(biāo)為(0,Risin θi,Ricos θi),換能器基陣在橫搖r和縱搖p的影響下,A點(diǎn)旋轉(zhuǎn)到了B點(diǎn),第i號波束的實(shí)際入射角度為(即∠BOD),定義經(jīng)旋轉(zhuǎn)后第i號波束的水平角度φi=∠BDE,即為波束橫距BD與OY軸的夾角,其表達(dá)式為
圖3 波束點(diǎn)空間旋轉(zhuǎn)示意圖Fig.3 Schematic diagram of footprint rotation
由式(18)可獲得姿態(tài)影響下的波束實(shí)際初始入射角,并用于后續(xù)聲線跟蹤。
波束測深點(diǎn)在地理坐標(biāo)系下的坐標(biāo)可通過如下過程來獲得:
1)波束實(shí)際初始入射角計(jì)算。由式(13)和式(16)得到換能器基陣面的實(shí)際旋轉(zhuǎn)角α和β,并由式(18)得到波束的實(shí)際初始入射角
2)常梯度聲線跟蹤。借助實(shí)際聲速剖面、1)中獲得的波束實(shí)際初始入射角以及波束傳播時(shí)間t,并根據(jù)式(1)~(6)進(jìn)行聲線跟蹤,獲得波束在理想換能器坐標(biāo)系下坐標(biāo)(x,y,z)。
3)波束傳播垂面與理想狀態(tài)換能器坐標(biāo)系XOZ面的二面角Ab計(jì)算:
4)波束測深點(diǎn)地理坐標(biāo)計(jì)算[12]
式中:下角b?LLS為地理坐標(biāo)系下的波束測深點(diǎn)坐標(biāo),T?LLS為由GPS提供的換能器中心地理坐標(biāo),TS為船體坐標(biāo)系下的波束測深點(diǎn)坐標(biāo);h為波束垂面方位角,由羅經(jīng)提供的航向角為A,則h為
為檢驗(yàn)上述方法的正確性,本實(shí)驗(yàn)運(yùn)用常梯度聲線跟蹤方法結(jié)合聲線精確歸位模型對膠州灣測深數(shù)據(jù)進(jìn)行了處理,本次實(shí)驗(yàn)的數(shù)據(jù)是應(yīng)用EM3000多波束測深儀器采集的,作業(yè)頻率為300 kHz,波束角為1.5°× 1.5°,127個(gè)波束,波束間隔為0.9°。膠州灣測深實(shí)驗(yàn)區(qū)的水深在30~50 m范圍。通過實(shí)驗(yàn)分析本文的聲線精確跟蹤方法與傳統(tǒng)方法的差異,并分析這兩種數(shù)據(jù)處理方法下的測深精度。
選擇重疊率達(dá)50%左右的兩條相鄰條帶以及一條中央波束近似垂直前兩條的條帶,由于多波束中央波束相對于邊緣波束受姿態(tài)影響小得多,在聲速剖面正確情況下可視為實(shí)際深度,因此以第3個(gè)條帶的中央波束作為檢測線(如圖4所示),作為其他2個(gè)條帶所對應(yīng)平面位置測深點(diǎn)深度的精度評估。對這3個(gè)條帶數(shù)據(jù)分別采用本文的精確歸位模型和傳統(tǒng)方法進(jìn)行處理,獲得各波束測深點(diǎn)在地理坐標(biāo)系下的三維坐標(biāo)。在兩條相鄰條帶中,分別經(jīng)過X軸和Y軸坐標(biāo)內(nèi)查得出檢測線所在平面位置的波束測深點(diǎn)水深,以檢測線的中央波束實(shí)測深度為參考,可評估受聲線跟蹤影響較大、經(jīng)本文方法和傳統(tǒng)方法分別處理后的水深精度,進(jìn)而實(shí)現(xiàn)對給出的聲線精確跟蹤方法的正確性評估。
若公共覆蓋區(qū)某點(diǎn)P的聲線跟蹤深度為Z,對應(yīng)平面位置檢測線中央波束深度為Z0,則聲線跟蹤深度的相對精度X為
圖4 條帶與檢測線交叉示意圖Fig.4 Schematic diagram of the strip is crossing with detection line
在圖4中,右側(cè)條帶為第1個(gè)條帶,左側(cè)條帶為第2個(gè)條帶,橫向中間線為第3個(gè)條帶的中央波束(檢測線)。第1個(gè)條帶與檢測線有18個(gè)交叉點(diǎn),第2個(gè)條帶有40個(gè)交叉點(diǎn),圖5和圖6分別為本文方法、傳統(tǒng)聲線跟蹤方法計(jì)算所得第1個(gè)條帶和第2個(gè)條帶與檢測線中央波束交叉點(diǎn)的深度精度??梢钥闯觯涸诘?個(gè)條帶中,本文給出聲線跟蹤方法計(jì)算的波束深度相對誤差在1%之內(nèi),而傳統(tǒng)方法有50%左右在1%之內(nèi),最大的相對誤差在1.6%左右;在第2個(gè)條帶中,本文給出的聲線跟蹤方法計(jì)算的波束深度相對誤差也控制在1%之內(nèi),而傳統(tǒng)方法則處于0.5%~2.7%。
基于本文給出的聲線跟蹤方法計(jì)算所得的波束測深點(diǎn)深度相對誤差較小,而應(yīng)用傳統(tǒng)方法處理則存在較大的波動,分析認(rèn)為前者已全面考慮了姿態(tài)角對波束初始入射角的影響,而后者未完全顧及船體姿態(tài)的影響,因此會隨著船體姿態(tài)的變化在不同條帶中產(chǎn)生較大的測深誤差。
圖5 兩種聲線跟蹤方法在條帶1所得深度相對誤差分布Fig.5 Relative depth errors of two kinds of sound ray tracing at swath 1
圖6 兩種聲線跟蹤方法在條帶2所得深度相對誤差分布Fig.6 Relative depth errors of two kinds of sound ray tracing at swath 2
由圖5和圖6可以看出,盡管顧及了姿態(tài)角的影響,本文給出的聲線精確跟蹤方法得到的波束測深點(diǎn)深度相對誤差仍出現(xiàn)所謂的“笑臉”現(xiàn)象,分析認(rèn)為盡管采用了測量位置的聲速剖面和消除了姿態(tài)對波束初始入射角的影響,但波束測深點(diǎn)深度跟蹤精度仍會受聲速剖面測量誤差影響以及其他影響因素的殘余誤差影響。盡管如此,仍達(dá)到了相對測深誤差控制在1%范圍內(nèi)的精度,滿足了IHO S?44標(biāo)準(zhǔn),遠(yuǎn)高于傳統(tǒng)聲線跟蹤方法,表明了本文提出的聲線精確跟蹤方法的正確性。
本文提出的顧及船體姿態(tài)角的聲線跟蹤模型,首先精確地計(jì)算出多波束波束的初始入射角,然后根據(jù)層內(nèi)常聲速聲線跟蹤模型計(jì)算出波束在海底的投射點(diǎn)坐標(biāo)。通過實(shí)驗(yàn),利用檢測線的平面位置內(nèi)插得到其他兩個(gè)相鄰條帶中與檢測線波束腳印相同位置的深度,再計(jì)算內(nèi)插波束點(diǎn)的深度相對誤差,實(shí)驗(yàn)表明提出的精確聲線跟蹤方法相對于傳統(tǒng)聲線跟蹤方法,顯著地提高了波束測深點(diǎn)深度的計(jì)算精度,滿足了IHO-44精度要求,對多波束高精度測深具有更好的參考性和實(shí)用意義。
[1]CHRISTIAN D M,MARTIN C K.Bathymetric artifacts in sea beam data:how to recongnize them and what causes them[J].Journal of Geophysical Research,1986,91(B3):3407?3424.
[2]EDOUARD K.A new method for the removal of refraction artifacts in multibeam echosounder systems[D].New Bruns?wick:University of New Brunswick,2000:35?60.
[3]DIMITRI A,CHRISTIAN D M.Adaptive noise canceling ap?plied to sea beam sideloe interference rejection[J].IEEE Journal of Oceanic Engineering,1988,13(2):70?76.
[4]趙建虎,劉經(jīng)南.多波束測深及圖像數(shù)據(jù)處理[M].武漢:武漢大學(xué)出版社,2008:125?127.
[5]CHAPPELL D J,GIANI S,TANNER G.Dynamical energy analysis for built?up acoustic systems at high frequencies[J].Journal of Acoustical Society of America,2011,130(3):1420?1429.
[6]劉伯勝,雷家.水聲學(xué)原理[M].哈爾濱:哈爾濱工程大學(xué)出版社,1993:23.
[7]THOMAS E B.Geometrically derived ray?theory results and direct verification of the Pekeris solution for unbounded con?stant?gradient media[J].IEEE Journal of Oceanic Engineer?ing,2012,37(2):244?254.
[8]HANAKO O,KAZUYOSHI M,TOSHIAKI N.Reciprocal sound propagation experiment in very shallow water area of Hashirimizu port[J].Japanese Journal of Applied Phys?ics,2010(49):151?156.
[9]HAMID R,HADI J R.Target localization and tracking for an isogradient sound speed profile[J].IEEE Transactions on Signal Processing,2013,61(6):1434?1446.
[10]MCCAFFREY E K.A review of the bathymetric swath sur?vey system[J].International Hydrographic Review,1981,VIII(1):19?27.
[11]陳非凡.多波束條帶測深儀的動態(tài)測量誤差評估[J].海洋技術(shù),1999,18(1):41?45.
CHEN Feifan.Estimation of dynamic measuring error for multibeam swath bathymeter[J].Journal of Ocean Tech?nology,1998,18(1):41?45.
[12]趙建虎,劉經(jīng)南.多波束測深系統(tǒng)的歸位問題研究[J].海洋測繪,2003,23(1):6?8.
ZHAO Jianhu,LIU Jingnan.Problems on conformity to the re?al sounding points from the multi?beam sounding system[J].Hydrographic Surveying and Charting,2003,23(1):6?8.
A precise multibeam sound ray tracking method taking into account the attitude angle
HE Linbang1,ZHAO Jianhu1,ZHANG Hongmei2,WANG Xiao1,YAN Jun1
(1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.School of Power and Mechanical Engineering,Wuhan University,Wuhan 430072,China)
The assigned beam initial incident angle provided by MBS(multibeam bathymetry system)is adopted by the traditional sound ray tracking method regarding the vessel in the horizontal status,which ignores or does not consider completely the impacts of vessel attitude angles,and lead obvious error to the determination of beam foot?print position.Aiming at the problem.In order to reduce this kind of significant error,an accurate method for sound ray tracking,which considers vessel attitude angles are put forward in this paper.In the proposed method,the actu?al propagation profile is given by means of studying the actual status of MBS transducer impacted by vessel attitude and the actual initial incidence angle of MBS beam is derived by the calculation model.By integrating the sound ve?locity profile,the accurate positioning of beam footprint on seabed is achieved by constant?gradient sound ray track?ing method.By comparing the beam footprint coordinate precision calculated by traditional sound ray tracking meth?od and the precise method of sound ray tracking in experiment,a conclusion can be drawn that the precise method of sound ray tracking can overcome the deficiency of the traditional sound ray tracking.This improves significantly the precision for multi?beam bathymetry
multibeam bathymetric system(MBS);beam footprint;initial incident angle of beam;attitude angle;constant?gradient sound ray tracking
10.3969/j.issn.1006?7043.201310060
P229
A
1006?7043(2015)01?0046?05
http://www.cnki.net/kcms/doi/10.3969/j.issn.1006?7043.201310060.html
2013?10?21.網(wǎng)絡(luò)出版時(shí)間:2014?11?07.
國家自然科學(xué)基金資助項(xiàng)目(41176068,0976061,40776048,41376109).
何林幫(1981?),男,博士研究生;
趙建虎(1970?),男,教授,博士生導(dǎo)師.
何林幫,E?mail:whuhlb@163.com.