王愛學(xué),趙建虎,尚曉東,張紅梅
(1.武漢大學(xué) 測繪學(xué)院,湖北 武漢430079; 2.東華理工大學(xué) 江西省數(shù)字國土重點實驗室,江西 南昌 330013; 3.武漢大學(xué) 動機學(xué)院,湖北 武漢 430072)
?
單波束水深約束的側(cè)掃聲吶圖像三維微地形反演
王愛學(xué)1,2,趙建虎1,尚曉東1,張紅梅3
(1.武漢大學(xué) 測繪學(xué)院,湖北 武漢430079; 2.東華理工大學(xué) 江西省數(shù)字國土重點實驗室,江西 南昌 330013; 3.武漢大學(xué) 動機學(xué)院,湖北 武漢 430072)
高分辨率海底微地形可呈現(xiàn)更豐富的海床細節(jié)特征,對水下目標(biāo)正確判讀具有重要意義。除去底質(zhì)因素影響,側(cè)掃聲吶圖像明暗變化與海底地形起伏變化有關(guān),挖掘側(cè)掃聲吶圖像中的微地形信息,可以彌補現(xiàn)有測深手段分辨率上的不足。針對SFS方法從側(cè)掃聲吶圖像反演出的相對地形缺少起伏尺度約束和起算基準(zhǔn)等問題,本文通過小波分析方法從單波束數(shù)據(jù)中分別提取與側(cè)掃聲吶圖像反演地形相關(guān)的高頻地形及表征地形趨勢變化的低頻地形,以此構(gòu)建側(cè)掃聲吶反演地形的尺度約束及基準(zhǔn)改正模型。試驗獲得的反演地形從整體趨勢到細節(jié)特征均與真實地形具有較好的一致性,統(tǒng)計結(jié)果表明,反演地形內(nèi)、外符合均方根誤差均小于15 cm,反演地形整體精度與單波束地形相當(dāng)。
側(cè)掃聲吶圖像;微地形反演;明暗形狀恢復(fù);約束方法;地形頻譜分析
隨著我國海洋開發(fā)向遠海深海推進,各種海洋活動對海床地形的質(zhì)量和分辨率也提出了更高的要求[1-2]。精細的海床地形可表征海床表面“麻坑”或“丘頂”等微觀特征[3-4],對發(fā)現(xiàn)天然氣水合物疑似埋藏區(qū)域、分析海底礦藏成因機制、研究海床沉積物搬運機理均具有重要意義。目前,高分辨率海底地形主要通過大比例尺單波束測量和多波束測量來獲取[5]。單波束獲取高分辨率海底地形,需要密集布線,測量效率低、成本高;多波束系統(tǒng)測量效率高,但隨著測量深度增加,其測深點間距增大,邊緣地形分辨率快速下降。與上述測深系統(tǒng)不同,側(cè)掃聲吶采用拖曳模式貼近海底,可獲得到厘米級分辨率海底地貌圖像,盡管其不具有直觀的地形信息,但去除底質(zhì)因素的影響,聲吶圖像的陰影及明暗強度均與地形變化存在明顯依賴關(guān)系。若能挖掘這些微地形信息,并將其與常規(guī)測深數(shù)據(jù)相結(jié)合,即可彌補現(xiàn)有測深系統(tǒng)在分辨率上的缺陷。計算機視覺領(lǐng)域采用明暗修復(fù)形狀(shape from shading,SFS)方法由圖像明暗強度反演目標(biāo)形狀[6],該方法基于Lambert法則和一系列的約束條件求解目標(biāo)在特定尺度下的最優(yōu)相對高度。為了求解SFS反演地形的絕對尺度并實現(xiàn)其與實測地形的融合,本文采用小波分析方法來尋找側(cè)掃聲吶圖像反演地形與實際地形的相關(guān)性,并以此為基礎(chǔ)提出了一種基于單波束水深約束的側(cè)掃聲吶圖像微地形反演算法。
SFS地形反演方法以Lambert漫反射法則為基礎(chǔ),認為入射強度一定的情況下,能量在粗糙表面的漫反射強度僅與能量入射方向、物體表面法線形成夾角θ的余弦有關(guān)[7-9]。物體或海底表面高度可用空間曲面函數(shù)z=f(x,y)描述(如圖1),物體表面法向量N(nx,ny,nz)可表示為高度的梯度形式N(p,q, -1),根據(jù)Lambert法則可列出聲能在海床表面的輻照度方程:
I(x,y) =R(p(x,y),q(x,y))=I0cosθ=
(1)
式中:I0為入射能量強度,S(ps,qs, -1)為能量入射方向,式(1)表征了聲能反射強度與地形梯度之間的關(guān)系,SFS問題的基本任務(wù)就是對給定的圖像I,根據(jù)每個像素灰度值I(x,y)求定被測物體表面上對應(yīng)位置的表面相對地形參數(shù),進而重構(gòu)物體的表明形狀z(x,y)。該問題求解還需一系列約束條件,式(2)為Zheng等[10]綜合了強度約束、強度梯度約束及可積性約束等條件形成的一種經(jīng)典最小化SFS泛函方程:
J=?F(p,q,z)=?Ω[(I-R)2+(Rx-Ix)2+(Ry-Iy)2+μ((zx-p)2+(zy-q)2))]dxdy=min
(2)
上述泛函極值問題可通過變分法轉(zhuǎn)化為歐拉方程,并采用有限差分方法建立迭代過程,近而求得圖像各處的最優(yōu)地形。大量試驗表明,SFS方法在計算機合成圖像的三維恢復(fù)中可反演出滿意的目標(biāo)形狀[11],但在真實的光學(xué)圖像或聲學(xué)影像三維恢復(fù)中,反演結(jié)果往往與真實形狀有一定的差異。其原因除了漫反射模型的理想化外,聲吶發(fā)射聲波的頻率和脈沖寬度一定,其采集的聲吶圖像明暗強度難以反映海底地形的整體變化趨勢,若無真實地形做初始約束時,SFS最小化法僅能解算出特定尺度(起伏范圍)和基準(zhǔn)(起算基值)下的相對微地形。
圖1 海床表面形狀模型Fig.1 Sea floor surface model
若將真實海床地形看成一個全頻段非平穩(wěn)隨機信號,那么該信號可經(jīng)小波變換在空間域分解為多個頻段地形,在底質(zhì)相同區(qū)域,側(cè)掃聲吶圖像的反演地形應(yīng)與其中某一頻段地形具有相關(guān)性。確定小波變換的最優(yōu)分解層數(shù),是尋找與反演地形最相關(guān)的真實地形頻段,并利用外部水深數(shù)據(jù)構(gòu)建反演地形約束模型的關(guān)鍵。
(3)
在真實地形逐層分解過程中,當(dāng)提取的高頻部分中心頻率與聲吶圖像反演地形的主要頻率相當(dāng)時,二者地形起伏變化達到最大相關(guān)性。因此,可以先對聲吶圖像的反演地形進行傅里葉分析,獲得其地形頻譜分布特征,并選其半功率譜對應(yīng)頻率作為小波變換預(yù)提取高頻信息的中心頻率f′,進而利用式(4),可求出需要分解的尺度a:
a=f*·fs/f′
(4)
N=lba=lb(f*fs/f′)
(5)
根據(jù)分解層數(shù)與分解尺度的關(guān)系按式(5)可求出最佳的分解層數(shù)N。
單波束地形分辨率主要受測線布設(shè)寬度限制,但其測線剖面具有較高的采樣率,可較好描述海底地形的整體趨勢和細節(jié)特征。因此,可選取少量單波束測線,采用小波分解方法分別提取與側(cè)掃聲吶圖像反演地形相關(guān)的高頻地形和低頻地形,利用高頻地形為同區(qū)域反演地形提供尺度約束,而低頻地形作為起算基準(zhǔn),將上述尺度和基準(zhǔn)用于無單波束水深的區(qū)域,從而得到整個區(qū)域的恢復(fù)地形。該地形反演約束方法的主要流程如圖2所示。
圖2 單波束地形約束的SFS地形反演方法流程圖Fig.2 Flow chart about SFS terrain recovering method constrained by single-beam soundings
3.1 尺度約束
忽略海底底質(zhì)、噪聲等因素的影響,側(cè)掃聲吶圖像的反演地形與單波束高頻地形的關(guān)系可用線性方程描述,若令單波束采樣序列為X,側(cè)掃聲吶圖像反演地形為f(X),單波束高頻地形為hhig(X),g(X)為尺度約束后的反演地形,三者的關(guān)系可用下式表示:
g(X)≈hhig(X)=kf(X)+d
(6)
尺度約束的目的是使反演地形與單波束高頻地形的起伏程度一致,而地形的標(biāo)準(zhǔn)差是衡量地形起伏程度的重要指標(biāo)[13-14]。因此,可利用地形的標(biāo)準(zhǔn)差來求解尺度縮放系數(shù)k。兩組地形序列在統(tǒng)計規(guī)律上應(yīng)該滿足如下關(guān)系:
(7)
3.2 基準(zhǔn)改正
結(jié)合式(6)、(7),可以求出反映真實尺度的反演地形。該地形還需添加正確的起算基準(zhǔn),才能得到最終的反演地形。根據(jù)小波分解原理,真實地形H(X)經(jīng)小波多層分解后,得到目標(biāo)頻段的高頻地形hhig(X),同時也獲得反映真實地形趨勢變化的低頻地形hlow(X),該低頻地形即為反演地形的起算基準(zhǔn):
(8)
整個約束反演的過程可用式(8)表示,H′(X)為經(jīng)過尺度和基準(zhǔn)約束后的最終反演地形。上述尺度與基準(zhǔn)求解方法適用于與單波束測線同區(qū)域的側(cè)掃聲吶圖像地形反演,在沒有單波束測線的區(qū)域,則需要采用插值方法,如移動曲面擬合、反距離加權(quán)、三角剖分等,對前文求得的尺度參數(shù)和起算基準(zhǔn)進行插值或延拓。
本文利用實測的單波束水深數(shù)據(jù)和同區(qū)域側(cè)掃聲吶圖像設(shè)計了外部水深約束的SFS地形反演方法試驗。選用深圳以南某水域為試驗區(qū),區(qū)域內(nèi)海底底質(zhì)相近,地形特征變化明顯;采用單波束測量對該區(qū)域海底地形進行精細測量,單波束采樣間隔為0.5 m,測線間隔為10 m,圖3為試驗區(qū)單波束測量地形圖;采用Edgetech 4200側(cè)掃聲吶采集同區(qū)域海底表面的聲吶圖像,圖4為一塊150 m×85 m的側(cè)掃聲吶圖像,圖像分辨率為300×170,像素尺寸為0.5 m,該聲吶圖像已經(jīng)過必要的畸變改正、灰度均衡等處理,聲吶圖像在試驗區(qū)中范圍如圖3所示。
圖3 所選區(qū)域海底地形及測線分布情況Fig.3 Seabed topography of test site and layout of survey lines
圖4 實測側(cè)掃聲吶圖像Fig.4 Side-scan sonar image measured in test area
在與聲吶圖像重合的區(qū)域內(nèi)共有9條單波束測線,如圖3所示,選取其中5條,即圖中的虛線測線,作為內(nèi)符合測線,用于構(gòu)建側(cè)掃聲吶圖像反演地形的約束模型,剩下的4條線,即圖中的實線測線,作為外符合檢測線。
圖5為側(cè)掃聲吶圖像經(jīng)最小化SFS數(shù)值算法計算得到的反演地形圖,從中可以看出,SFS算法根據(jù)圖像的明暗變化得到圖像區(qū)域的相對地形,反演地形的高度在0上下浮動,尺度浮動范圍為-1×10-4~1×10-4m,顯然該反演地形并不具有正確的起伏尺度和絕對的起算基準(zhǔn)。接下來用5條內(nèi)符合測線對其進行尺度和起算基準(zhǔn)修正。
圖5 側(cè)掃聲吶圖像的SFS反演地形圖Fig.5 SFS recovering topography from side-scan sonar image
根據(jù)5條內(nèi)符合單波束測線及每條測線上各點的坐標(biāo),在反演地形中插值同位置處的反演值,從而得到5對位置對應(yīng)的實測單波束地形剖面和反演地形剖面序列。分別對每條單波束地形剖面進行小波分解,本文選用bior3.3小波函數(shù),按前所述方法,對各單波束地形分解至第7層。圖6分別列出了5條內(nèi)符合測線的實際地形剖面、7層小波分解后的低頻地形、7層小波分解后的高頻地形以及對應(yīng)位置上的反演地形剖面,圖中橫坐標(biāo)為剖面采樣序號。
從圖6中可看出,經(jīng)7層小波分解,從5條內(nèi)符合單波束測線中均能較好的分離出高、低頻地形。從各單波束高頻地形及其對應(yīng)的反演地形起伏變化來看,在同一位置處,反演地形與高頻地形具有一致的變化趨勢,這說明側(cè)掃聲吶圖像近似反映了海底地形的高頻細節(jié)特征,基于聲吶圖像反演的海底地形與實際地形的高頻部分存在較強的線性相似關(guān)系。
利用式(6)、(7)求出尺度縮放參數(shù)k及尺度平移參數(shù)d,如表1所示。將上述5條測線求解得到的尺度參數(shù),通過合理的插值方法延拓至整個反演區(qū)域。本試驗采用線狀單波束測深作為外部約束數(shù)據(jù),測線與測線之間近似平行分布,測線間隔約為20 m,鑒于外源數(shù)據(jù)的分布特征,可以根據(jù)插值點至相鄰測線之間的垂直距離,采用反距離加權(quán)法,求解相鄰測線之間區(qū)域的尺度參數(shù)。
表1 各測線的尺度參數(shù)
根據(jù)各位置上求得的尺度參數(shù),對整個反演地形進行尺度修正,尺度約束后的反演地形圖如圖7所示。從該圖可以看出,經(jīng)過尺度約束以后,反演地形的浮動范圍從-1×10-4~1×10-4m修正至-0.4~0.4 m的正常尺度范圍。
接下來,利用從5條單波束水深中提取的低頻地形數(shù)據(jù),采用帶線性插值的三角剖分法構(gòu)建整個反演區(qū)域的起算基準(zhǔn)面,如圖8所示,該起算基準(zhǔn)面反映了海底地形的整體變化趨勢。
將圖7中經(jīng)過尺度修正后的反演地形與起算基準(zhǔn)面相加,得到最終的海底反演地形,如圖9所示。圖10為經(jīng)過密集單波束精細測量獲得的試驗區(qū)域海底地形圖,其分辨率與側(cè)掃聲吶圖像相同,可以認為該地形是真實的海底的地形,與之相比,圖9中的反演地形不僅具備正確的整體變化趨勢,而且保留了側(cè)掃聲吶圖像豐富的細節(jié)特征,反演地形中的細節(jié)特征與真實地形的細節(jié)特征具有相同的起伏程度和較好的相似性。
為了定量的分析該算法的地形恢復(fù)效果,對上述9條測線序列,計算約束反演地形與同位置處單波束實測地形之間的差值,并統(tǒng)計每條測線和所有測線總體的均方差、平均誤差、最大誤差、最小誤差等相關(guān)精度參數(shù),內(nèi)、外符合精度統(tǒng)計情況分別見表2和表3。
圖6 單波束測線小波分解及與反演地形剖面的比較Fig.6 Wavelet decomposition of single-beam soundings and comparison with recovering terrain
圖7 尺度約束后的反演地形圖Fig.7 Recovering topography constrained by scale correction
圖8 低頻地形生成的起算基準(zhǔn)面Fig.8 Reference surface generated by low frequency terrain
圖9 尺度和基準(zhǔn)約束后的最終恢復(fù)地形圖Fig.9 Final recovering topography with scale and reference correction
圖10 同區(qū)域單波束精細測量地形圖Fig.10 Elaborate terrain surveyed by single-beam sonar
表3 外符合精度統(tǒng)計表
從內(nèi)符合精度統(tǒng)計表可以看出(見表2),最終的反演地形在5條內(nèi)符合測線處的均方根誤差皆小于15 cm,平均誤差均為零,這主要是因為本文采用地形的標(biāo)準(zhǔn)差來計算尺度縮放參數(shù),采用線性變換進行尺度修正,這保證了尺度約束前后反演地形的統(tǒng)計規(guī)律不會發(fā)生變化。從外符合精度統(tǒng)計表可以看出(見表3),在沒有單波束水深區(qū)域,通過尺度和基準(zhǔn)的延拓,最終的反演地形精度與內(nèi)符合區(qū)域相當(dāng),均方根誤差同樣小于15 cm,每條外符合測線的平均誤差也較小,達到了厘米級甚至毫米級。圖11和12分別為內(nèi)、外符合概率分布曲線,從中可以看出,整個反演地形的誤差近似服從標(biāo)準(zhǔn)正態(tài)分布,基本上不存在系統(tǒng)誤差。上述統(tǒng)計指標(biāo)均滿足IHO海道測量規(guī)范的要求。
圖11 內(nèi)符合概率密度曲線Fig.11 Internal probability density curve
圖12 外符合概率密度曲線Fig.12 External probability density curve
1)試驗及統(tǒng)計結(jié)果表明,本文所提出的采用少量單波束測深數(shù)據(jù)對聲吶圖像的反演地形進行約束改正的思路是可行的、正確的,經(jīng)過尺度和基準(zhǔn)約束以后的反演地形從整體趨勢到細節(jié)特征均與真實地形保持了較好的一致性和準(zhǔn)確性,達到了由側(cè)掃聲吶圖像恢復(fù)高分辨率海底微地形的目的。
2)該方法僅需要少量外部水深數(shù)據(jù),便可利用側(cè)掃聲吶圖像反演出精度可靠的、高分辨率的海底地形,大大較少了大比例尺單波束測量的成本,也可彌補多波束邊緣測深分辨率隨水深增加而下降的缺陷。
由于側(cè)掃聲吶圖像的明暗變化除受地形因素影響外,還與底質(zhì)、入射角、環(huán)境噪聲等有關(guān),海底并非是理想的漫反射體,SFS方法用于聲吶圖像海底地形反演有其局限性,尤其是側(cè)掃聲吶近場區(qū)域圖像并不滿足朗伯體法則。因此,基于側(cè)掃聲吶圖像開展微地形恢復(fù),應(yīng)避開天底大掠射角區(qū)域。采用文中所提方法開展大區(qū)域微地形反演時,側(cè)掃聲吶數(shù)據(jù)測線間隔應(yīng)保證50%以上覆蓋率,而單波束測線間隔應(yīng)根據(jù)實際海底地形變化程度合理選擇,一般一條側(cè)掃聲吶圖像應(yīng)在左右兩側(cè)各配合一條單波束測線進行建模。
[1]金翔龍. 海洋地球物理技術(shù)的發(fā)展[J]. 東華理工學(xué)院學(xué)報, 2004, 27(1): 6-13.
JIN Xianglong. The development of technique of marine geophysics[J]. Journal of East China Institute of Technology, 2004, 27(1): 6-13.
[2]陸秀平, 黃謨濤, 翟國君, 等.多波束測深數(shù)據(jù)處理關(guān)鍵技術(shù)研究進展與展望[J]. 海洋測繪, 2016, 36(4): 1-6.
LU Xiuping,HUANG Motao, ZHAI Guo jun, et al. Development and prospect of key technologies formultibeam echosounding data processing[J]. Hydrographic surveying and charting, 2016, 36(4): 1-6.
[3]王健, 邱文弦, 趙俐紅. 天然氣水合物發(fā)育的構(gòu)造背景分析[J]. 地質(zhì)科技情報, 2010, 29(2): 100-106.
WANG Jian, QIU Wenxian, ZHAO Lihong. Tectonic settings analysis of gas hydrate deposits development[J]. Geological science and technology information,2010, 29(2): 100-106.
[4]鄧希光, 吳廬山, 付少英,等. 南海北部天然氣水合物研究進展[J]. 海洋學(xué)研究, 2008, 26(2): 67-74.
DENG Xiguang, WU Lushan, FU Shaoying, et al. The research advances of natural gas hydrates in northern south China Sea[J]. Journal of marine science, 2008, 26(2): 67-74.[5]董慶亮, 歐陽永忠, 陳岳英,等. 側(cè)掃聲納和多波束測深系統(tǒng)組合探測海底目標(biāo)[J]. 海洋測繪, 2009, 29(5): 51-53.DONG Qinliang, OUYANG Yongzhong, CHEN Yueying, et al. Measuringbottom of sea target with side-scansonarandmultibeam soundingsystem[J]. Hydrographic surveying and charting,2009, 29(5): 51-53.
[6]HORNB K. Obtaining shape from shading information[M]. Massachusetts: MIT press, 1989.
[7]HORN B K. Shape from shading: a method for obtaining the shape of a smooth opaque object from one view[R]. US-MIT, 1970, 232.
[8]HORNBK.Height and gradient from shading[J]. International journal of computer vision, 1990, 5(1): 37-75.
[9]HORNBK,SZELISKIRS, YUILLEAL. Impossible shaded images[J]. IEEE transactions on pattern analysis and machine intelligence, 1993, 15(2): 166-170.
[10]ZHENG Qinfen, CHELLAPPA R.Estimation of illuminant direction, albedo, and shape from shading[J]. IEEE transactions on pattern analysis & machine intelligence, 1991: 13(7): 680-702.
[11]ZHANG Ruo, TSAI Pingsing, CRYER J E, et al. Shapefromshading: a survey[J]. Pattern analysis and machine intelligence, 1999, 21(8): 690-706.
[12]樊計昌, 劉明軍, 海燕,等. 計算尺度函數(shù)和小波函數(shù)中心頻率的 GUI 及其應(yīng)用[J]. 科技導(dǎo)報, 2008, 25(24): 36-39.
FAN Jichang, LIU Mingjun, HAI Yan, et al. GUI for computing center-frequency of the scale and waveletfunctionsandits applications[J].Science & technology review,2008, 25(24): 36-39.
[13]隋剛, 郝兵元, 彭林. 利用高程標(biāo)準(zhǔn)差表達地形起伏程度的數(shù)據(jù)分析[J]. 太原理工大學(xué)學(xué)報, 2010, 41(4): 381-384.
SUI Gang, HAO Bingyuan, PENG Lin. Data analysis of elevation standard deviation classifying relief degree of land surface [J]. Journal of Taiyuan University of Technology, 2010, 41(4): 381-384.
[14]劉曉萌, 常占強. 數(shù)字高程模型數(shù)據(jù)分布規(guī)律的研究[J]. 河北師范大學(xué)學(xué)報: 自然科學(xué)版, 2006, 30(4): 473-477.
LIU Xiaomeng, CHANG Zhanqiang. Research on the distribution law of DEM data [J]. Journal of Hebei Normal University: natural science edition,2006, 30(4): 473-477.
本文引用格式:
王愛學(xué), 趙建虎, 尚曉東, 等. 單波束水深約束的側(cè)掃聲吶圖像微地形反演[J]. 哈爾濱工程大學(xué)學(xué)報, 2017, 38(5): 739-745.
WANG Aixue,ZHAO Jianhu,SHANG Xiaodong,et al. Recovery of seabed 3D micro-topography from side-scan sonar image constrained by single-beam soundings [J]. Journal of Harbin Engineering University, 2017, 38(5): 739-745.
Recovery of seabed 3D micro-topography from side-scan sonar image constrained by single-beam soundings
WANG Aixue1,2,ZHAO Jianhu1,SHANG Xiaodong1,ZHANG Hongmei3
(1.School of geodesy and Geomatics, Wuhan University, Wuhan 430079, China; 2.Jiangxi Province Key Lab for Digital Land, East China Institute of Technology, Nanchang 330013, China; 3.School of Power and Mechanical Engineering, Wuhan 430072, China)
Without considering the impact of marine sediment, the gray intensity of a side-scan sonar image mainly depends on the changing seabed micro-topography, which can be used to compensate for the low-resolution deficiency of traditional bathymetry methods. In order to recover the micro-topography and transform it into a valid scale and absolute reference, we used the wavelet analysis method to extract high- and low-frequency terrains from single-beam soundings, in which the high-frequency terrains correspond with topography recovered from side-scan sonar images and the low-frequency terrains represent topography change. On this basis, we constructed a size restraint and reference correction model for side-scan sonar topography recovery. The test results show that the final recovered topography is consistent with the actual terrain in overall trend and local detail. Statistical results show that the internal and external root mean square error (RMSE) are both less than 15 cm, and the final recovery accuracy is the same as that for single-beam soundings.
side-scan sonar image; micro-topography recovery; shape from shading; constrain method; spectrum analysis of topography
2015-12-17.
日期:2017-04-26.
國家自然科學(xué)基金項目(41376109,41576107,41606114);江西數(shù)字國土重點實驗室開放基金項目(DLLJ201502);精密工程與工業(yè)測量國家測繪局重點實驗室開放基金項目(PF2015-7).
王愛學(xué)(1984-), 男, 講師, 博士.
王愛學(xué),E-mail: axwang@sgg.whu.edu.cn.
10.11990/jheu.201512057
P229.1
A
1006-7043(2017)05-0739-07
網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20170426.1031.010.html