一種適于海浪測量的立體攝影物理模型?

王英霞 姜文正 喬方利 陳思宇

1)(中國海洋大學海洋與大氣學院,青島 266100)

2)(國家海洋局第一海洋研究所海洋環境與數值模擬研究室,海洋環境科學與數值模擬國家海洋局重點實驗室,青島 266061)

3)(青島海洋科學與技術國家實驗室區域海洋動力學與數值模擬功能實驗室,青島 266237)

一種適于海浪測量的立體攝影物理模型?

王英霞1)2)3)姜文正2)3)喬方利2)3)?陳思宇2)3)

1)(中國海洋大學海洋與大氣學院,青島 266100)

2)(國家海洋局第一海洋研究所海洋環境與數值模擬研究室,海洋環境科學與數值模擬國家海洋局重點實驗室,青島 266061)

3)(青島海洋科學與技術國家實驗室區域海洋動力學與數值模擬功能實驗室,青島 266237)

(2016年10月17日收到;2016年11月25日收到修改稿)

由于在海面上布設控制點較為困難,本文發展了一種無需海面控制點的立體攝影物理模型.該模型依次通過相機固有參數標定、相對定向和基于平均海面的絕對定向三個環節實現,文中給出了各個環節的控制方程和具體算法.利用室內實驗時拍攝的一個影像對,將共面方程和距離方程聯立構成閉合方程組,實現相機固有參數的標定;利用外海測量時拍攝的一個影像對,采用共面方程實現相對定向參數的標定;利用外海測量時拍攝的一個影像對序列,在左相機坐標系中通過影像匹配確定平均海面,建立特定的物方坐標系,并實現絕對定向.室內實驗檢驗了固有參數控制方程的收斂性并表明定標精度較高;外海實驗成功地構建了海面高度場,這表明本文提出的相對定向和絕對定向方法是可行的.該模型只需將兩個相機按一定方位架設在平臺上即可進行海浪測量,免去了繁瑣的海面控制點布設過程,降低了工作難度和成本,有利于立體攝影海浪測量技術的推廣和應用.

立體攝影,海浪,固有參數,控制點

1 引 言

由于在海氣相互作用、海洋上層混合和海洋工程方面的重要作用[1?3],海浪一直是物理海洋學研究的熱點.外海觀測獲取第一手資料是深入開展海浪理論及其應用研究的基礎,方向譜的觀測是現代海浪觀測的核心.立體攝影能夠直接測量波面高度場及其演化過程,是一種非常有效的方向譜觀測手段,它可以彌補目前海洋觀測設備的不足[4,5].得益于數碼相機和計算機的發展,近幾年立體攝影海浪測量技術獲得了顯著的提高,并逐漸成為一種實用化的海浪測量手段[6?19],如文獻[16—19]利用立體攝影測量的三維海浪數據重新審視了海浪的動力學特性,研究了海氣相互作用機制,得出了一些從浮標等單點測波數據中無法得到的結論.在此基礎上,2016年Bevgamasco等[20]開展了船載立體攝影海浪測量技術研究,在理論上提出了一種可行的船載立體攝影物理模型.

物理模型的構建是立體攝影研究的主要課題之一,受陸地近景攝影測量的影響,目前大部分立體攝影海浪測量系統是建立在直接線性變換理論[21]基礎上的,較為典型的有Wanek等[6,7]發展的全自動化三目立體攝影海浪觀測系統,Vries等[11,12]發展的立體攝影海浪測量系統等.直接線性變換理論定標時需要物方(海面)控制點,如Vries等的定標方法是將帶有全球定位系統(GPS)的小船駛入攝影區域,以GPS天線為控制點,并調整相機與GPS時間同步,通過GPS獲取控制點大地坐標系坐標.在海面上布設控制點是困難的,因此基于直接線性變換理論發展的測量系統定標過程非常繁瑣.近幾年,Benetazzo等[8,9]和姜文正等[22]各自提出了無需海面控制點的立體攝影海浪觀測技術,兩者都是基于平均海面建立物方坐標系并據此實現無海面控制點的系統定標.在具體方法上兩者又有些區別,比如Benetazzo采用理想的針孔相機模型,不考慮畸變參數;而姜文正采用一般非測量相機,考慮了徑向畸變但定標精度較低,穩定性較差.本文在文獻[22]的研究基礎上,構建了一種適于海浪測量的立體攝影物理模型,介紹了一種新的相機固有參數定標模型及其算法,改進并提高了基于平均海面的立體攝影系統定標方法的精度和穩定性.

2 物理模型與算法

2.1 相機固有參數的標定

2.1.1 相機成像模型

采用非測量相機,假定鏡頭的主光軸垂直感光傳感器平面,像素排列的行方向與列方向互相垂直,像素之間的間隔相等.事實上,這種假定對絕大多數相機而言是合理的.為了便于數學描述,本文定義兩個右手笛卡爾坐標系(如圖1所示),其一是以相機主點為坐標原點,x軸、y軸分別平行于像素排列的行方向和列方向,z軸沿主光軸指向物方,稱為相機坐標系;其二是以相機感光傳感器中心為坐標原點,u,v軸分別平行于x,y軸,稱為影像坐標系.

圖1 相機成像的物理模型Fig.1.The physical model of camera imaging.

設(xm,ym)是不考慮畸變時物點M(xM,yM,zM)的像點在相機坐標系中的坐標(如圖1所示),(um,vm)是點M的實際像點在影像坐標系中的坐標,則(xm,ym)與(um,vm)的關系可以表示為

其中x0,y0是影像坐標系原點在相機坐標系中的坐標,δxm,δym是畸變修正項,可近似為[23]

其中k1表示徑向畸變系數,

表示像點到像主點的距離.本文將(1)式中含有的3個參數x0,y0,k1和鏡頭焦距f稱為相機固有參數.

2.1.2 相機固有參數的控制方程

本文定義面向相機鏡頭時左側的相機為左相機,右側的相機為右相機.攝影測量時,若物點A(xA,yA,zA)在左、右相機感光傳感器上分別形成像點則a,a′和基線oo′共面(如圖2所示).設左相機坐標系o?xyz連續繞y,x,z軸旋轉βy,βx,βz角度后與右相機坐標系o′-x′y′z′對應坐標軸相互平行; 再設o-xyz連續繞y,x軸旋轉αy,αx角度后,z軸與基線重合,則共面方程可表示為

其中(xa,ya,?f)表示像點a在左相機坐標系的坐標,(ξa′,ηa′,ζa′)滿足關系式

其中(x′a,y′a,?f′)表示像點a′在右相機坐標系的坐標,R0表示坐標系o-xyz與o′-x′y′z′之間的轉換矩陣,

本文將(5)式中含有的參數βy,βx,βz,αx,αy稱為相對定向參數.

方程(3)中含有5個相對定向參數,8個相機固有參數.每一個同名點對可列出一個方程(3),但從文獻[24]可知,不可能僅憑增加方程(3)的個數求得上述13個參數,還需要其他限制條件.

如圖2所示,在左相機坐標系中A,B兩點距離可表示為

其中S表示A,B兩點間的距離,需通過測量得到;(xA,yA,zA)是直線oA,o′A的交點坐標,可表示為

其中(ξa′,ηa′,ζa′)由(4)式給定, 設基線的長度oo′為D,則

(3)式和(6)式中含有5個相對定向參數αx,αy,βx,βy,βz,8個相機固有參數x0,y0,f,k1,x′0,y′0,f′,k′1和基線長度D共14個參數.本文采用附有參數的條件平差方法確定這些參數.

2.1.3 控制方程的解算過程

對固有參數x0,y0,x′0,y′0做如下數學變換:

其中f0,f′0分別表示左右相機焦距的初始估計值.

圖2 立體攝影成像模型Fig.2.Stereo photography imaging model.

以參數

為平差參數,像點的影像坐標u,v,u′,v′為觀測值,將(3)式、(6)式線性化:

以上迭代過程的初始值可分兩步確定:第一步,f,f′取鏡頭焦距,θx,θy,θ′x,θ′y,k1,k′1取零,5個相對定向參數由共面方程的直接解確定[25];第二步,測量兩相機鏡頭中心點間距離D,將其作為已知量進行αx,αy,βx,βy,βz,f,f′,θx,θy,θ′x,θ′y,k1,k′1等13個參數迭代,迭代結果和D即為上述迭代過程的初始值.多次實驗表明采用該方法得到的初始值時迭代收斂.

2.2 相對定向

相機架到海上平臺上后,室內標定的相對定向參數不再適用,需要基于海浪測量時拍攝的影像對重新標定.本文采用(3)式作為解算相對定向參數的控制方程.以像點影像坐標u,v,u′,v′為觀測值,以參數為平差參數,將(3)式線性化:

其中F1(i)表示第i對同名點滿足的共面方程(3),n表示同名點對的個數,由于海面影像匹配的精度較低,本文建議n>300.這樣1的值為

2.3 絕對定向

海浪測量時每組影像由一序列等間隔的影像對組成,選取任一影像對,采用文獻[25]所述的影像匹配方法進行匹配得到300×300個近似均勻的同名點對的影像坐標.如圖2所示,在左相機坐標系中與同名點對a,a′相應的直線oA,o′A可表示為

其中(xa,ya,?f)是像點a的左相機坐標系坐標,(1)式給出了它與像點a的影像坐標之間的關系,(ξa′,ηa′,ζa′)由(4)式給定,是基線的三個分量.

4個方程3個未知數,A點的坐標的平差值可表示為

按上述方法求出所有同名點對對應的物點坐標,然后將這些點擬合成一個平面ax+by+z+c=0,擬合原則為各物點到平面距離的平方和最小:

其中N=90000表示同名點對的個數.這樣得到一個近似為平均海面的平面P.

用同樣的方法求出該影像組所有影像對對應的平面,然后進行疊加得到平均海面Pm.將平均海面Pm的方程表示為amx+bmy+z+cm=0的形式,則Pm與z軸交點坐標為的單位法線n為

其中ex,ey,ez表示左相機坐標系x,y,z軸方向的單位矢量,n的球面坐標系坐標(1,φ,θ)為:θ=

本文將物方坐標系O-XY Z定義為以左相機坐標系z軸與平均海面Pm的交點為坐標原點O,Z軸垂直于Pm平面并指向上方,X軸垂直于z軸與Z軸確定的平面并在面向左相機時指向右側(如圖3所示).顯然,將物方坐標系O-XY Z繞X軸旋轉θ,再繞Z軸旋轉0.5π?φ后,其與左相機坐標系對應的坐標軸平行,因此兩坐標系之間的轉換關系為

其中x,y,z表示左相機坐標系坐標,X,Y,Z表示物方坐標系坐標,而

至此,本文通過相機固有參數標定、相對定向和絕對定向建立了一個完整的無需海面控制點的立體攝影物理模型和算法.

圖3 相機坐標系與物方坐標系間的幾何關系Fig.3.Geometrical relationship between camera coordinate system and object coordinate.

3 實驗與討論

3.1 相機固有參數標定實驗

實驗在國家海洋局第一海洋研究所海洋環境與數值模擬研究室進行,時間為2015年11月3日.實驗利用兩臺SI-6600CL工業相機及KOWA LM8XC型8.5 mm鏡頭同時采集影像;相機分辨率為2210×3002,像素大小為3.5μm.如圖4所示,在實驗區的水平地面上鋪設沙脊,沙脊上有AA′,BB′,CC′,DD′,EE′,FF′,GG′7條線段.表1給出了線段的長度及端點的影像坐標,其中u,v表示左影像坐標系坐標,u′,v′表示右影像坐標系坐標,S表示線段的長度;影像坐標以像素大小為單位,線段長度以mm為單位.

將圖4左影像均勻分成10×10個分區,在每一個分區內提取一個特征點,并通過影像匹配找到其右影像中的同名點.這100個同名點對加上7條線段的14個端點共可列出114個(3)式,7條線段可列出7個(6)式.這121個方程即為求解14個參數的條件方程,表2給出了參數的初始值和最終的平差值結果.

圖4 固有參數標定時采用的影像對 (a)左相機拍攝的影像;(b)右相機拍攝的影像Fig.4.Image pair used for calibration of intrinsic parameters:(a)Image taken by the left camera;(b)image taken with the right camera.

表1 線段端點影像坐標(單位:像素)及線段長度(單位:mm)Table 1.Endpoint image coordinates(unit:pixel)and line length(unit:mm).

表2 14個參數的初始值及平差值Table 2.Initial value and adjustment value of 14 parameters.

為了檢驗14個參數標定的精度,將圖4左影像均勻分成30×30個分區,在每個分區內提取一個特征點,并通過與14個參數無關的影像匹配方法找到其在右影像上的同名點,由圖2可見,同名點應在核線上.將14個參數代入(3)式求出核線解析表達式,并計算同名點橫坐標對應的核線的縱坐標值.圖5(a)給出了各同名點縱坐標與相應核線縱坐標的差值.從圖中可見,除極少數的錯誤匹配點(差值大于1個像素的點)外,各點的差值均勻對稱地分布在線周圍.圖5(b)給出了差值的概率密度分布,其數學期望為?0.01像素,標準差為0.21像素;圖中的虛線表示相同標準差的正態分布曲線,可見差值的概率密度分布近似于正態分布,但比正態分布精度更高.顯然,這些差值是由偶然誤差引起的,因此可以認為所有同名點都在核線上.由于同名點和核線是由相互獨立的方法得到的且差值的期望值為?0.01像素,精度為0.21像素,所以14個參數的標定是準確的且精度較高.

圖5 (a)同名點縱坐標與相應核線縱坐標的差;(b)差值分布的概率密度,實線表示差值的實際概率密度分布,虛線表示相同標準差的正態分布Fig.5.(a)The difference between the ordinate of the same point and the ordinate of the corresponding epipolar line;(b)the probability density of the difference distribution,the solid line is the actual probability density distribution,and the dotted line is the normal distribution of the same standard deviation.

3.2 相對定向的實驗檢驗

外海實驗在廣東省博賀氣象-海洋綜合觀測平臺進行,時間為2015年11月11日.博賀平臺位于111?23′E 21?26′N, 距岸邊約7 km, 水深約16 m.相機安裝在風塔距海面約21 m處,左相機主光軸指向?40?(如圖6所示).拍攝前6 h內風速8—10 m/s,風向50?—70?,為向岸風.

實驗使用兩臺SI-6600CL型相機拍攝海面影像,其固有參數見表2,圖7為拍攝的影像對之一.將左影像均勻分成40×40個分區,在每一個分區內提取一個特征點,并通過影像匹配找到其同名點.由于海面影像的紋理信息少,影像匹配的精度和準確率較低,本文采用文獻[25]的識別和剔除方法,在1600個同名點對中選取最優匹配的407個點對計算相對定向參數,表3給出了由共面方程直接解得到的初始值和基于(12)式計算的平差值.

圖6 實驗地點和相機安裝位置Fig.6.Experimental site and camera installation location.

圖7 實驗拍攝的海面影像對 (a)左相機拍攝的影像;(b)右相機拍攝的影像Fig.7.Pair of sea images taken during the experiment:(a)Image taken by the left camera;(b)image taken by the right camera.

表3 相對定向參數的初始值和平差值(單位:rad)Table 3.Initial value and adjustment value of Relative orientation parameters(unit:rad).

3.3 絕對定向的實驗檢驗

如圖7所示,在左影像的白框內均勻提取300×300個特征點,依次采用基于核線的金字塔影像匹配,基于核線的最小二乘影像匹配,錯誤匹配點識別與剔除和補充匹配得到右影像上的同名點.由(16)式計算出對應的海面上的物點坐標,圖8給出了在左相機坐標系中這些物點的分布.

用平面ax+by+z+c=0按(17)式進行擬合,得到近似平均海面P的單位法線和截距.外海觀測時,相機連續采集256對間隔為0.8 s的影像序列,按上述方法確定出256個近似平均海面的單位法線和截距;所有單位法線矢量和的歸一化值即為平均海平面的單位法線,所有截距的平均值即為平均海平面的截距.計算單位法線的球坐標系坐標,然后根據(20)式和(21)式求出轉換矩陣,結果如下:這里矩陣T的單位為m.

至此左相機坐標系與物方坐標系之間的轉換關系完全確定,將所有左相機坐標系內的物點坐標轉化到物方坐標系中,稍做平滑后得到圖9所示的海面起伏的三維高程圖.

圖8 左相機坐標系中海面上物點的分布Fig.8.Distribution of object points on the sea surface in left camera coordinate.

圖9 海面起伏的高程圖Fig.9.Elevation map of the sea fluctuation.

由圖9可見,前后兩個波峰夾著波谷,這與圖7影像顯示的相同(細線位置是波谷);波峰與波谷間的波面位移差約2 m,這與實驗時8—10 m/s的風速引起的風浪大小相符,因此海面起伏被成功重構.這說明本文提出的無海面控制點的絕對定向方法是可行的,影像匹配也是準確的.由于所有影像匹配都是基于核線進行的,因此由相對定向參數計算的核線也是準確的,這從側面驗證了本文提出的相對定向方法是準確的.由圖9亦可見海面點的高程都是相對于平均海面的,所以本文提出的立體攝影物理模型剔除了潮汐對海浪測量的影響并且可以直接測量潮位變化.

4 結 論

本文發展了一種適合海浪測量的立體攝影物理模型,提出了一項利用共面方程和距離方程標定相機固有參數的方法,針對海浪這種特定的測量物,發展了穩定而精確的相對定向和絕對定向方法.實驗室內利用沙脊和沙脊上的7條線段計算出了精確的相機固有參數,這一方面表明了該方法的可行性,另一方面也表明與需要專業定標場的傳統方法相比,該方法更簡潔、方便.相對定向時,本文在1600個同名點對中選取最優匹配的407個同名點對進行相對定向,保證了相對定向的穩定性和精度;絕對定向時,本文發展了利用影像對序列確定平均海平面的方法,保證了確定平均海面的精度;這些方法在外海實驗中都得到了檢驗.與以往立體攝影海浪測量模型相比,本文提出的模型不需要在海面布設控制點,標定完固有參數的相機按一定方位架設到平臺上即可進行海浪測量,操作簡單、經濟;同時模型以平均海面為基準構建波面高度場,數據直觀且消除了潮位對攝影測量的影響,并能直接測量潮位變化.

[1]Melville W K 1996Annu.Rev.Fluid Mech.28 279

[2]Thorpe S A 1995Prog.Oceanogr.35 315

[3]Qiao F L,Yuan Y L,Yang Y Z,Zheng Q A,Xia C S,Ma J 2004Geophys.Res.Lett.31 293

[4]Krogstad H E,Barstow S F 1999Proceedings of 9th International Offshore and Polar Engineering ConferenceBrest,France,May 30–June 4,1999 p19

[5]Fedele F,Benetazzo A,Gallego G,Shih P C,Yezzi A,Barbariol F,Ardhuin F 2013Ocean Mode.70 103

[6]Wanek J M,Wu C H 2006Ocean Eng.33 723

[7]Bechle A J,Wu C H 2011Coastal.Eng.58 305

[8]Benetazzo A 2006Coastal.Eng.53 1013

[9]Benetazzo A,Fedele F,Gallego G,Shih P C,Yezzi A 2012Coastal.Eng.64 127

[10]Benetazzo A,Barbariol F,Bergamasco F,Torsello A,Carniel S,Sclavo M 2016Coastal.Eng.109 114

[11]Vries S,Hill D F,Schipper M,Stive M 2011Coastal.Eng.58 239

[12]Vries S,Hill D,Schipper M,Stive M 2009J.Coastal.Res.56 1484

[13]Gallego G,Yezzi A,Fedele F,Benetazzo A 2011IEEE Trans.Geosci.Remote Sens.49 4445

[14]Kosnik M V,Dulov V A 2011Meas.Sci.Technol.22 015504

[15]Liu P C 2013Int.J.Geosci.4 229

[16]Banner M L,Barthelemy X,Fedele F,Allis M,Benetazzo A,Dias F,Peirson W L 2014Phys.Rev.Lett.112 114502

[17]Benetazzo A,Barbariol F,Bergamasco F,Torsello A,Carniel S,Sclavo M 2015J.Phys.Oceanogr.45 2261

[18]Leckler F,Ardhuin F,Peureux C,Benetazzo A,Bergamasco F,Vladimir D 2015J.Phys.Oceanogr.45 2484

[19]Mironov A S,Yurovskaya M V,Dulov V A,Hauser D,Guérin C A 2012J.Geophys.Res.117 C00J35

[20]Bergamasco F,Benetazzo A,Barbariol A,Carniel S,Sclavo M 2016Comput.Geosci.95 105

[21]Abdel-Aziz Y I,Karara H M 1971Proceedings of the ASP/UI Symposium:Close-Range PhotogrammetryUrbana,USA,January 26–29,1971 p1

[22]Jiang W Z,Yuan Y L,Wang Y X 2012Acta Phys.Sin.61 119101(in Chinese)[姜文正,袁業立,王英霞 2012物理學報61 119101]

[23]Wang Z Z 1997Principle of Photogrammetry(Wuhan:Wuhan University Press)p235(in Chinese)[王之卓1997攝影測量原理(武漢:武漢大學出版社)第235頁]

[24]Zhang Z X,Zhang Q J 1997Digital photogrammetry(Wuhan:Wuhan University Press)p103(in Chinese)[張祖勛,張劍清 1997數字攝影測量學(武漢:武漢大學出版社)第103頁].

[25]Wang Y X 2015Ph.D.Dissertation(Qingdao:Ocean University of China)(in Chinese)[王英霞2015博士學位論文(青島:中國海洋大學)]

PACS:92.10.Hm,91.10.Lh,92.10.Kp DOI:10.7498/aps.66.059201

A stereo photographic physical model for ocean wave measurement?

Wang Ying-Xia1)2)3)Jiang Wen-Zheng2)3)Qiao Fang-Li2)3)?Chen Si-Yu2)3)

1)(College of Oceanic and Atmospheric Sciences,Ocean University of China,Qingdao 266100,China)
2)(Ocean Environment and Numerical Simulation Laboratory,The First Institute of Oceanography,SOA;Laboratory of Marine Science and Numerical Model in SOA,Qingdao 266061,China)
3)(Laboratory for Regional Oceanography and Numerical Modeling,Qingdao National Laboratory for Marine Science and Technology,Qingdao 266237,China)

17 October 2016;revised manuscript

25 November 2016)

Unlike the traditional models,a new stereo photographic model requiring no sea control points that are difficult to lay on a sea surface is developed in this paper.It is realized according to the order of camera intrinsic parameters calibration,relative orientation and absolute orientation based on average sea surface,and we give the governing equations and its algorithms.First of all,in the paper we present an imaging model that adopts non-measurement camera with considering only the radial distortion coefficients,and then give the method that calculates the camera intrinsic parameters by combining the co-planar equations and the distance equations to form a closed system of equations.For verifying the convergence of the governing equations,we lay four parallel sand ridges on a level ground in laboratory,measure seven lines of different heights and take a pair of the images by two cameras.The intrinsic parameters are successfully calculated by using the method,and further analysis shows that the results are highly precise.When observing the ocean waves in offshore sea,we need to recalculate the relative orientation parameters.At this time the governing equations are the co-planar equations,but the accuracy of the image matching is seriously related to the accuracy of the relative orientation parameters.However,as the mirror reflection of the sea surface,the matching accuracy of the sea image is often low,so we select 407 optimal matching points from among 1600 conjugate points and calculate the parameters.It needs a sequence of sea image pairs to obtain the absolute orientation parameters,256 pairs in the paper.For each pair,we select firstly 300×300 feature points on the left image approximately evenly,then find the conjugate points on the right image by image matching,calculate the coordinates of the corresponding sea surface points in the left camera coordinate system,and then fit the points into a plane;thus a plane sequence can be obtained.We establish a special object coordinate system on the average sea surface that is obtained by summing the plane sequence,and calculate the absolute orientation parameters by the geometric relationship between the two coordinate systems.And it is proved that the methods of relative orientation and absolute orientation proposed in this paper are feasible by reconstructing sea surface height field.The model greatly reduces the difficulty in calibrating the stereoscopic photography for ocean wave measurement,which is beneficial to its popularization and application.

stereoscopic photography,ocean wave,intrinsic parameters,control points

PACS:92.10.Hm,91.10.Lh,92.10.Kp

10.7498/aps.66.059201

?國家自然科學基金委員會-山東省人民政府聯合資助海洋科學研究中心項目(批準號:U1406404)、國家自然科學基金青年基金(批準號:41506041)和國家海洋局第一海洋研究所基本科研業務費專項基金(批準號:GY02-2012T01)資助的課題.

?通信作者.E-mail:qiao fl@fio.org.cn

*Project supported by the NSFC-Shandong Joint Fund for Marine Science Research Centers,China(Grant No.U1406404),the National Natural Science Foundation of China(Grant No.41506041),and the Basic Research Foundation of the First Institute of Oceanography,SOA,China(Grant No.GY02-2012T01).

?Corresponding author.E-mail:qiao fl@fio.org.cn