一种适于海浪测量的立体摄影物理模型∗

王英霞 姜文正 乔方利 陈思宇

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