论文推荐|曹彬才：高分辨率卫星立体双介质浅水水深测量方法

Original 2016-10-08 曹彬才，邱振戈 测绘学报

本文提出了一种基于高分辨率卫星多光谱立体像对的浅水水深测量方法。该方法利用多光谱近红波段消除太阳耀斑，以有理函数模型（RFM）构建测区原始DEM，通过水陆边界内插获得水面高程，并采用双介质近似折射改正模型消除目标点垂直坐标偏移。试验表明，本文的模型方法和处理流程在水面平静、底质纹理丰富的浅海岛礁水深反演中能取得优于20%的相对测深精度，可为浅水水深测量提供新手段。

浅海水深探测方法主要有声呐测深、SAR地形反演、LiDAR测深、多光谱/高光谱水深反演以及双介质摄影水深测量。声呐测深以测量船为运载平台，通过声呐测深仪向海底发射声波并接受回波获得海底地形^[]，SAR工作的微波波段虽然无法直接穿透水体，但可以通过探测海水表面粗糙度的变化间接反演水深^[]，LiDAR利用蓝绿激光发射接收设备，通过大功率、窄脉冲激光探测海底深度^[]。声呐虽然精度高，但受制于浅海区水深，即使是多波束系统也无法获得大幅宽，这种方法花费高、效率低，且有一定危险性，并不适于大面积浅海测绘；SAR依据流经水下地形的潮流对海洋表面波的水动力调制反演水下地形，要求风速3~8m/s，潮流>0.5m/s，极大地限制了该技术的实用性^[]；LiDAR可以获得密集的、高精度测深点，如SHOALS-3000最大探测深度50m，水平精度2.5m，测深精度25cm^[]，一旦突破了高应用成本的制约，将成为浅海测深的主流方法；目前，利用卫星多光谱影像反演浅海水深最为普遍，基本原理是通过少量的已知水深点建立水下辐射强度与水深值的函数关系，与前3种方法相比具有经济、直观的优势，缺点是需要水深控制点，测深精度受大气、水质及底质类型等影响^[-]，当观测条件不佳时，仅可作为补充手段，解决信息有无的难题^[]；高光谱水深反演的本质与多光谱方法类似，其光谱分辨率更强，有利于识别反射信号微弱的水体信息，在底质类型差异较大区域的水深反演中更有优势^[-]，但高光谱数据通常几何分辨率低、数据量大，需要进行波段配准、数据降维等预处理。

双介质摄影测量技术并非新生事物，从20世纪80年代起开始得到关注，基本原理是利用物方空间和像方空间处在两种不同介质中拍摄的图像确定被摄目标几何特性^[-]。文献[13—14]研究了双介质摄影测量的一般构象关系和基本公式，文献[]阐述了双介质相对定向的流程，文献[16—17]系统分析了波浪对双介质测深精度的影响，以上研究均采用模拟数据，并未给出双介质立体测水深的实际案例。文献[18—20]推导了近似折射改正公式，使用航空立体像对在低流速、清洁的卵石河床条件下获得了高质量的水下地形信息。

目前，双介质摄影测量在隧道监测^[]、水底测图^[]、考古^[]及工业测量^[]方面都有成功应用的案例，但以下难点制约了双介质水深测绘的精度：①水面位置必须已知或可以建模；②太阳耀斑和白浪会影响水下特征匹配；③水面不平和立体摄影延时造成的表面变化影响双介质几何结构等。近年来，IKONOS、WorldView-2等高分辨率敏捷卫星的出现为双介质浅海测绘提供了新思路，高精度无控立体定位可以克服多光谱水深反演依赖控制点的缺点，快速重访周期下可方便选择理想的大气、风速、潮流等拍摄条件。

本文探讨了目前尚未见诸报道的高分辨率卫星立体影像浅海双介质水深测绘。首先推导了双介质折射改正公式，定量分析近似改正公式的误差大小，随后通过耀斑消除、立体模型构建、水下点匹配、粗差控制及折射改正获得水下DEM，最后利用海图数据和LiDAR实测数据进行精度评估。试验结果表明，利用风浪小、水底纹理清晰的高分辨率卫星影像进行双介质测深，可以获得优于20%的相对测深精度，具备在一定条件下代替多光谱水深反演的潜力。

1 立体像对双介质测深基本原理

摄影测量学根据两个已知的摄影站点和两条已知的摄影方向线，交会出构成这两条光线的待定地面点的三维坐标。双介质摄影测量的原理与之类似，差别在于光线经过了空气/水两种不同的媒介会产生折射作用，如果仍然按照直线传播的几何结构建模解算，目标点三维坐标高于实际位置，从而造成明显的系统误差。如所示，由于光线的折射作用，水下目标点P沿着折线PP₁S₁、PP₂S₂被摄站位于S₁、S₂处的立体影像记录到。双介质摄影测量就是探测同名光线交点(由影像匹配获得观察点A)并恢复目标点真实三维坐标(点P)的过程。

图 1 双介质摄影基本几何关系Fig. 1 Basic two-medium geometric relationships

光线从一种介质进入另一种介质时产生的折射现象可以用折射定律概括，即

(1)

式中，r₁、r₂分别为左、右摄站到水下点的入射角；i₁、i₂为对应的折射角；n为水的折射率，通常取1.340，该值在较大的温度、盐度变化范围内误差小于1%^[]。

2 双介质折射改正与测深误差分析2.1 近似折射改正模型

为了简化折射模型，假设水质均一，折射率为常数，水面为平面。如所示，坐标系统假设沿飞行方向为X轴，构成右手坐标系，将立体像对重叠区分为航迹方向、航迹垂直平分线及其他位置3种情况展开讨论。

2.1.1 航迹方向

如所示，同名折射点P₁、P₂的连线与航迹方向平行，航高H。假设观察点A不在S₁、S₂的垂直平分面上，即D₁≠D₂，A点与S₁、S₂的连线确定了入射角r₁、r₂，由式(1)可以计算出射角i₁、i₂。在比h更深的任意高程面h_C处，x_C₁=h_C·tani₁，x_C₂=h_C·tani₂，P₁P_C₁和P₂P_C₂的交点即为目标点的真实位置P。

由于折射点P₁、P₂位于航迹方向，故Y_P=Y_A。设P₁P₂的距离为k，则

(2)

(3)

由几何关系和折射定律可得

(4)

(5)

(6)

(7)

真实位置的X大小关系为

(8)

将式(3)代入式(8)可得

(9)

式(9)表示了航迹方向目标点真实位置和观测位置的在平面上的定量关系，带入相关参数可以定量计算水平偏差ΔX，计算结果(具体见下文)显示航空/航天情况时该水平偏差都可以忽略。

2.1.2 航迹垂直平分线方向

如所示，当同名折射点P₁、P₂位于中航迹垂直平分线方向时，θ₁=θ₂，D₁=D₂，因此i₁=i₂，r₁=r₂，直线L₁、L₂必然交于一点，即X_A=X_P、Y_A=Y_P。

图 2 航迹垂直平分线方向几何关系Fig. 2 Geometry of two-medium photogrammetry in the perpendicular bisector of the in-track direction

因此有

(10)

将式(4)—式(7)代入可以解得真实水深值h。

2.1.3 其他位置

为了简化问题，忽略入射角差异造成的水平位置差(下文的定量分析解释其合理性)，即所有的点都遵循所示的几何关系，只是此时i₁≠i₂，k表示P₁P₂长度，有

(11)

因此

(12)

此处

(13)

将式(4)—式(7)和式(13)代入式(12)可以解得真实水深h。实际应用中并不一定采用这种同时考虑了左右片入射角度的双侧改正模型，文献[]指出双介质折射改正精度与其几何结构密切相关，越靠近地底点越容易受到水面不平的影响，因此常用入射角度适中的单张影像进行近似折射改正。仅考虑的左侧部分，式(1)变形得

(14)

式(3)变为

(15)

式(14)代入式(15)得

(16)

中L₁可用目标点A和摄影中心计算

(17)

简单的几何关系可得

(18)

式(18)代入式(16)得到A点真实水深。

2.2 测深误差分析

2.2.1 近似折射改正带来的水平误差

采用式(9)定量分析近似折射改正模型的水平误差。示例了航高1500m，基线长度500m条件下航线方向的不同点的观测位置与真实位置的水平方向差异，水深值显示了1m、15m和30m 3种情况。30m是可见光影像能观测到海底特征的极限深度，对应的最大水平偏差仅为0.04 m。此外，由于卫星影像具有高轨道、窄视场的特征，代入相关参数后计算结果表明航天情况的水平偏差比航空要小5个数量级。由于海洋测绘的误差容忍度大，因此对航空/航天立体数据而言，近似折射改正造成的水平偏差可以忽略。

图 3 航迹方向真实位置和观测位置的x误差Fig. 3 Difference in x values of observed and true positions for the in-track case

2.2.2 平均入射角带来的垂直误差

高分辨率遥感卫星普遍采用线阵CCD推扫成像，即影像中不同行相同列像素对应的光线入射角(地底角)相同，对相同行不同列像素而言，入射角略有差异但相差不大。示意了CCD沿飞行方向前后摆扫的几何关系，假设轨道高H，地面幅宽W，地面中心处入射角α₀已知，则

图 4 线阵CCD推扫成像几何关系Fig. 4 The geometric relationship of linear pushbroom CCD

(19)

例如WorldView-2的轨道高度700km，幅宽17.4km，设α₀=30°，则α₁≈29.9°，假设水深h_A=30m，代入式(16)可得由于入射角差异造成的高程差Δh=h₁－h₀≈45.048m－42.034m=0.014m，实际上只要卫星在一定角度范围内前后或侧方摆扫，同一行像素对应的入射角度差异都不大，由此造成的折射改正高程差可以忽略，即整景卫星影像的所有像元可以使用相同的平均入射角计算水深改正值。

2.2.3 其他误差

其他误差包括立体成像模型误差、同名点误匹配、波浪起伏、拍摄延时误差等因素。成像模型方面，可利用地面控制点(GCPs)改善三维定位精度，如WorldView-2无控定位精度为平面4.1m CE90%^[]，使用3个以上GCPs可获得优于1.5m CE90%；误匹配方面，太阳耀斑和不同的摄影角度影响了水下纹理的真实性，越深处纹理越匮乏，特征点可靠性降低，像点坐标的误匹配会导致真实点P的观测三维坐标A的直接计算错误，并且折射改正会放大该垂直误差，例如水深17m时，1个像素的误匹配最大会导致3.8m的测深误差；此外，波浪起伏导致的法线偏移，立体像对拍摄延时都会对测深造成较大影响，其中75%的测深误差来自水面波浪^[]。

3 高分辨率卫星影像双介质测深流程

如所示，运用高分辨率卫星立体像对进行双介质水深测量主要包括辐射处理、几何模型构建、水陆分离、水面高程计算、折射改正等步骤，由此得到摄影时刻的瞬时水深，再经过潮汐改正以获得符合制图需求的水深信息。

图 5 高分辨率卫星影像双介质测深处理流程Fig. 5 High resolution satellite imagery based two-medium photogrammetry bathymetry procedure

3.1 辐射处理

除大气校正、影像增强外，太阳耀斑是双介质水深测绘中需要考虑的重要问题。受海浪表面状态、太阳位置、拍摄角度等因素的影响，水面可能直接将太阳光线反射到传感器中。对分辨率优于10m的影像而言，波浪特征能得到清晰反映，太阳耀斑常表现为零散的亮斑或沿波浪边缘的白色条带，不利于立体像对水下同名特征点的提取。

采用Hedley算法^[]进行耀斑消除，该算法基于两点假设进行设计：①水体对近红外(NIR)波段有强烈的吸收性，浅水区(如数十厘米)足够使得近红波段的辐射强度趋近于零；②水体辐射强度仅仅与像点的地理位置有关，并且可见光波段(VIS)的最亮/最暗值基本对应NIR波段的最亮/最暗值，可以用回归直线表示VIS/NIR对应关系。Hedley算法使用深水区的影像子集计算VIS/NIR回归参数，改正后的i波段为

(20)

式中，b_i表示回归直线的斜率；L(NIR)为近红外波段灰度值; L_min(NIR)为目标子集内近红外波段灰度最小值；L_i(VIS)和L_i(VIS)′分别为纠正前后可见光波段灰度值。

3.2 几何处理

为了隐藏核心设计参数，影像供应商通常采用与传感器成像几何无关的、非严格的有理函数模型(RFM)用于高分辨率遥感卫星影像测绘处理^[-]。RFM将像点坐标(r_n, c_n)描述为以相应点地面坐标(X_n, Y_n, Z_n)为自变量的多项式比值，公式为

(21)

式中，a_ijk、b_ijk、c_ijk、d_ijk表示有理函数参数(RPCs)；n表示点个数；i+j+k定义了模型的次数，通常不大于3次，即i+j+k≤3。为了避免解算中参数级差过大，需要将原始物、像点坐标进行平移和缩放处理，使得标准化后的取值位于(－1.0，+1.0)之间。RFM立体定位时，将式(21)变形为

(22)

式中，F(X_n, Y_n, Z_n)、G(X_n, Y_n, Z_n)分别表示式(21)中r_n、c_n对应的等式右侧，将式(22)泰勒级数展开得误差方程

(23)

迭代解算中需要的初始值可以使用标准化平移参数的平均值，也可以用RFM一次项求解值。通过RFM立体解算陆地及水底匹配点三维坐标，建立原始水陆DEM。

3.3 水陆分离

原始水陆DEM构建以后，可通过水陆分离、边缘检测获得水陆分界像元，在原始水陆DEM中内插并取平均求得水面高程，以作为水下折射改正的起算位置。水陆分离使用NIR波段(在水中几乎无反射，图像强度与陆地有明显差异)，主要方法有区域生长法、阈值分割法、水体指数法等。阈值分割和水体指数都可能造成碎斑，而识别填充碎斑比较费时，因此以人工指定的深水区种子点进行区域生长的方法在水陆分离中更加有效，具体的区域生长门限可以由最佳分割阈值和经验确定。

4 试验验证4.1 试验影像

WorldView-2单程即可立体采集分辨率为2m的4波段多光谱影像，确保图像质量的连续性和一致性^[]。由于全色影像记录了从可见光波段到近红外波段的总辐射量，水下特征在全色影像中远不如蓝、绿单波段影像中清晰，因此双介质测深不宜使用全色影像。本文选取了4组采集时间间隔小、水质清洁、水下纹理丰富的WorldView-2多光谱立体像对(StereoOR2A级)验证卫星双介质立体测绘的可行性，试验区地理位置如所示，影像详细信息列于。

图 6 包含蓝、绿、红、近红波段并显示为真彩色的WorldView-2多光谱试验影像地理位置Fig. 6 Location of study area and shown on WorldView-2 true color

表 1 试验影像详细信息Tab. 1 Specific information about experimental imagery

立体像对		时间	平均入射角/(°)	影像质量	验证数据
永兴岛	左	2011-08-12(03:27:51)	1.5	严重耀斑	海图Navionics 采集时间、精度未知
右	2011-08-12(03:29:03)	31.0	无耀斑
甘泉岛	左	2014-04-02(03:33:31)	31.8	无耀斑	SHOALS-3000 LiDAR数据， 2013-1-19，XY2.5m，Z25cm
右	2014-04-02(03:32:15)	4.9	严重耀斑
黄岩岛	左	2012-10-12(03:16:21)	34.3	无耀斑	无
右	2012-10-12(03:14:50)	12.8	轻度耀斑
永暑礁	左	2014-05-31(02:59:08)	33.6	无耀斑	无
右	2014-05-31(03:00:40)	18.3	无耀斑

4.2 耀斑消除

以永兴岛和黄岩岛为例，显示了、中白框部分的耀斑污染和消除效果。原始影像中随波浪起伏的亮斑基本消失，水下纹理特征更加清晰，耀斑消除前的白线处光谱曲线上下波动十分剧烈，耀斑消除后光谱曲线趋于平稳。

图 7 太阳耀斑消除效果图Fig. 7 The effect of sunglint correction

4.3 水下点匹配与粗差剔除

如果水下特征在单波段图像上清晰可见，特征点匹配与常规的图像匹配并无差别，本文采用SIFT算子提取同名点。需要注意两点：①由于光在水体中的衰减，越深的地方纹理越匮乏；②深水区(可见光透水范围内的相对概念)与浅水区同名点的视差差别较大。

为了在深水区提取出特征点，可设置较大的SIFT最邻近/次邻近比率阈值，但会增加浅水纹理丰富区的误匹配数目，需要细致地剔除粗差。本文采用RANSAC法与局部视差比较相结合的粗差剔除方案：首先RANSAC整体粗差剔除，由于深水区与浅水区同名点视差差别大，因此RANSAC阈值不能太小，否则所有深水区或浅水区特征点可能被当作粗差剔除；随后，由于地形的连续性(除开断崖情况)，局部范围内视差(或高程)不会发生突变，可以凭此设置阈值剔除残余粗差。

以永兴岛和甘泉岛为例，如、所示，当SIFT比率阈值为0.5时，提取出特征点多数位于陆地部分、浅水区或深水区特征明显的位置；当阈值增大时，部分深水区也提取出大量的、位置准确的同名点。粗差剔除如所示，过小的RANSAC阈值将深水区正确同名点剔除，阈值2.0时保留了深水区特征点。

图 8 匹配及粗差剔除Fig. 8 Image matching and Gross error elimination

图 9 水陆分界与边缘提取Fig. 9 Water land separation and edge detection

4.4 折射改正与精度验证

经过影像匹配后得到同名点后，代入RFM立体定位模型获得三维坐标，建立起原始水陆DEM。随后对多光谱影像的近红波段进行区域生长和边缘检测()，获得水陆边缘，在原始DEM中内插水陆边缘位置高程，取平均作为海水面高程。

水下点折射改正时，选择几何结构更好的接近30°的入射角()，获得折射改正DEM。定量精度评价方面，由于缺乏黄岩岛、永暑礁的实测水深，仅对永兴岛、甘泉岛进行精度分析。永兴岛水深验证数据来自公开海图Navionics，该区域内有0.5m等高线的声呐图层，海图数据的平面精度较差，选择数个明显特征点(如海岛外围曲线拐点)完成与卫星影像的配准，随后均匀选择水深验证点(所示)；甘泉岛的LiDAR数据精度较高，内插中匹配点位置的水深值作为检查。

折射改正后DEM与原始影像叠加显示如所示，从目视效果看，地形起伏符合灰度、纹理的变化趋势，与二维视觉判断结果相吻合。定量分析结果如、和所示。中的ME、MAE、RMSE、REP、C²分别表示平均误差(mean error)、平均绝对误差(mean abs error)、均方根误差、相对误差百分比(relative error percent)和相关系数平方，其中C²在统计学中表示总离差平方和中可以由回归平方解释的比例，数据介于0~1之间，越接近1表示水深测量的精度越高，各评价参数计算公式如下

图 10 折射改正后的三维地形Fig. 10 3D view after refraction correction

图 11 永兴岛测深精度分析Fig. 11 Two-media Photogrammetry bathymetry accuracy for Yongxing Island

图 12 甘泉岛测深精度分析Fig. 12 Two-media Photogrammetry bathymetry accuracy for Ganquan Island

(24)

(25)

(26)

(27)

(28)

式(24)—式(28)中，x_i表示双介质测量水深值；y_i表示水深验证值；n为点总数。

由可知，折射改正后永兴岛的所有评价指标都大幅提升，得到了20%的相对测深精度，C²从0.61提升为0.89。左下方是与LiDAR数据同时获取的甘泉岛高分辨率航拍影像，结合LiDAR实测水深可将岛屿外划分为具有明显分界的浅水沙滩(水深小于1m)和深水礁盘(水深约1~20m)两部分，上方是WorldView-2立体像对，影像1水面平静，影像2有明显波浪，并且浅水沙滩部分波浪起伏强于深水礁盘。将折射改正后的检查点叠加在卫星影像上，可以明显看出大误差基本位于浅水沙滩，如右下侧所示，红色点相对误差大于50%，绿色点小于50%。由于波浪大造成了较大的测深误差，因此精度分析时不考虑浅水沙滩部分。列出了深水礁盘处检查点(绿色点)的精度情况，折射改正后可获得16%的相对精度。并且由和可以看出：折射改正前的估计水深值整体小于实测水深，折射改正让水深绝对值整体变大，这与前文推导的几何结构关系相符。此外，列出了试验区多光谱水深反演(双波段比值模型，具体见文献[])的结果，可以看出双介质方法取得了和多光谱反演相当的精度，甚至在REP和C²指标上更优。

表 2 双介质与多光谱测深精度比较Tab. 2 Accuracy comparison of two-media and multispectral methodm

m
方法	影像	折射改正	ME	MAE	RMSE	REP(%)	C²
双介质	永兴岛	否	4.25	4.42	3.25	37	0.61
是	0.53	1.95	2.44	20	0.89
甘泉岛	否	3.33	3.87	3.55	38	0.46
是	-0.04	1.51	1.90	16	0.88
多光谱	永兴岛		0.34	1.88	2.35	21	0.87
甘泉岛		0.21	1.44	1.76	18	0.86

4.5 试验误差分析

永兴岛与甘泉岛总体取得了优于20%的相对测深精度，如果影像上无明显的水陆分界(本文永暑礁、黄岩岛成像时水面淹没整个礁盘)，需要人工判断水面高程，可能造成测深误差；甘泉岛深水区部分匹配点可靠性差，难以被当作粗差剔除，也引入了部分误差；波浪方面，甘泉岛测深误差在浅海沙滩部分明显大于水面平静的深水礁盘区域，有理由推测双介质技术在风平浪静的情况下(如黄岩岛类的潟湖)能获得更佳的精度；验证数据方面，永兴岛用电子海图作为检查数据，与影像的配准、无潮位数据等可能带来误差，相比而言甘泉岛的精度评价更加客观、可靠。

5 建议与结论

理论上只要水下纹理清晰，能够组成立体像对的数据都可以用作双介质测深。为了控制影像质量，在卫星影像数据源选择时建议考虑以下方面：

(1)连续多光谱立体影像获取能力。全色影像包含的水下信息很少，无法用作双介质测绘；多数遥感卫星搭载了多光谱甚至高光谱传感器，但并不是所有卫星都有敏捷的立体成像能力，例如TH-1、ZY-3等国产卫星采用三线阵全色影像获得地理信息，但其多光谱数据为一次推扫成像，虽然可以通过异轨成像获得多光谱立体，但时间间隔久、水面位置变化大，不利于双介质测深。类似于WorldView-2等敏捷商业卫星则满足连续获取多光谱影像的能力。

(2)适中的立体拍摄角度。卫星影像的拍摄角度以适中为宜，过大或过小容易造成严重太阳耀斑。最佳的方案是多线阵推扫式的立体多光谱成像，但目前没有此类专门考虑海洋测绘需求的测绘卫星。

通过本文的理论分析和试验讨论，可得出如下结论：双介质摄影测量中，近似折射改正模型造成的水平差异可以忽略，即所有目标点都能当作光线相交的情形来近似折射改正；利用风浪小、纹理清晰的高分辨率卫星影像可以获得优于20%的相对测深精度。中国南海岛屿水下可见度高，随着卫星重访周期、敏捷性的逐步提高，选择风浪小的拍摄条件或存档影像资料不再遥不可及，本文提出的利用遥感卫星双介质测量浅海水深的方法在判断地貌走势、描述地貌特征、建立三维海底景观上具有一定应用前景。

全文内容请点击下方“阅读原文”，建议开启微信优化阅读。

曹彬才^1,2, 邱振戈¹, 朱述龙², 涂辛茹¹, 曹芳³, 曹斌¹