-
立体匹配是一种从不同视角获取影像中寻找像点对应关系的技术，已有匹配算法可以分为局部匹配算法、全局匹配算法和半全局匹配算法三类[1]。全局匹配算法以全局能量最小化为目标，理论上可以得到全局最优解，但计算复杂度极高，目前还没有有效的解决办法。Hirschmüller提出的半全局优化算法同样以全局能量最小化为目标，平滑项以P1、P2两个参数作为视差变化为1和视差变化大于1的惩罚项，通过多个方向的置信度传播和匹配代价聚合实现全局能量最小化[2-4]。大量实验表明，该算法执行效率较高，已经在近景、航空、卫星等多个领域的影像立体匹配中得到应用[3-6]。半全局优化算法要求计算整幅影像视差范围内所有视差对应的匹配代价，设最大视差范围为z，影像长、宽分别为h、w，需要计算的匹配代价个数为h×w×z。对于航空或卫星立体影像来说，z对应于影像范围内的高度差。显然在地形起伏较大的山区高度差较大，而在平坦地区高度差较小。如果平坦地区存在高大建筑，则局部最大高差由最高的建筑决定。由于半全局匹配的匹配代价传播来自于相邻的像素，因此在已知局部搜索范围的前提下，只进行局部搜索范围对应的匹配代价即可实现置信度传播的目的，那些对应整幅影像搜索范围且不在目标像素搜索范围内的匹配代价计算则是不必要的。本文通过金字塔影像结构逐层匹配，为下一层影像提供搜索范围的初值。结合改进的平滑参数设置方法实现了高效的半全局优化。资源三号为三线阵立体测绘卫星，对应同一区域存在正视、前视和后视3个视角的影像，可以构成3个立体像对。多视匹配一般分为两种策略，一种是采用逐像对匹配加上左右匹配结果一致性检验，需要进行4次或6次独立的匹配，然后将匹配结果融合，这种匹配一般以核线立体像对为基础，不需要通过地面投影寻找对应像点，通过左右一致性检验可以检测出大部分匹配误差；另一种策略是通过物方高程变化形成的投影轨迹寻找对应匹配点[7]，由于多幅影像上的同名像点对应同一地面点，在正确的搜索高程处匹配代价最小，而且这时的匹配代价是所有匹配点的综合结果，可靠性较高，虽然在遮挡区域会形成一定的干扰，但由于遮挡可能只存在正视-前视和正视-后视中的一种，所以仍可能获得正确的匹配结果。同时由于半全局优化算法在存在匹配点的区域对平滑参数变化不敏感，而在不存在匹配点的区域如遮挡或大面积纹理缺失区域对平滑参数的变化比较敏感，利用这一点结合局部纹理分析和高程分析可以探测出误匹配点。本文以物方方案为研究对象，提出了适合资源三号三线阵影像数字表面模型(digital surface model，DSM)提取的技术方案，并以高分辨率的激光探测与量测(light detection and ranging,LiDAR)数据获取的高精度DSM和高精度数字正射影像图(digital orthophoto map,DOM)为控制和检验数据进行了实验论证。
-
本文提出的DSM自动提取方法以物方半全局优化为基本匹配策略，结合可疑区域检测与精化匹配方法，实现了资源三号三视影像DSM提取。DSM提取流程如图 1所示。
图 1 DSM提取流程图
Figure 1. Flowchart of DSM Generation
(1) 三线阵影像及有理多项式参数(rational polynomial coefficient,RPC)导入。根据文件名区分正视、前视、后视影像和对应RPC，RPC已经过控制点校正。
(2) 影像预处理与纹理分析。由于资源三号正视影像与前后视影像分辨率差异较大，因此需要预先将影像变换为近似的分辨率，本文采用了两种策略，一种是将正视影像分辨率降低一倍，这时只需在生成影像金字塔时处理即可，另一种是将前后视影像插值成与正视相同的分辨率。影像预处理的另一个关键步骤是确定影像的最小有效值，并将低于该值的灰度值统一设为最小有效值，这样可以实现这些低值区域的平滑过渡，否则会出现一些虚假纹理影响匹配。资源三号影像的一个突出问题是在灰度过低的区域如起伏较大的山区背光区域，容易出现较大面积无纹理或假纹理(表现为横条纹)，这些区域通过基于局部窗口方差计算检测出来，这些区域的匹配代价统一设为一个固定数值，通过匹配代价传播步骤可以实现平滑过渡。
(3) 生成影像金字塔。通过小波变换生成影像金字塔，存储在数据体的对角区域。这里的关键是金字塔的层数及设置的顶层搜索范围是否与地区内的高差对应。
(4) 基于CENSUS[8]的顶层匹配。CENSUS(一种影像局部结构描述算子)操作方法如下：将一个窗口内所有像素的灰度与中心点的灰度比较大小，大于中心灰度的赋值为1，否则赋值为0，变换后的窗口描述算子为由0和1组成的二值数字串。两个CENSUS变换数字串中对应位置数值不同(0、1或1、0) 的个数称为CENSUS距离，距离越小说明两个窗口越相似。对顶层影像做CENSUS变换并进行基于半全局优化的匹配，探测可疑区域，若出现并非纹理缺失导致的大面积可疑区域，且可疑区域的高程都接近最大或最小高程，则认为搜索范围设置过小。将搜索范围加大并再次对误差区域进行匹配，若误差区域不再表现为大面积区域，则顶层匹配结束，若区域减小且还有较大面积误差区域，则仍需加大搜索范围，直到不存在因搜索范围导致的可疑区域(前后两次匹配结果无明显差异)。匹配的结果是正视影像每个像点的高程。
(5) 基于CENSUS距离和互信息(mutual information，MI)进行逐层匹配、可疑区域检测与精化匹配。采用分块匹配，相邻两块影像有一定的重叠区域，对顶层以下的每一层影像逐层进行CENSUS变换，并利用上一层每个正视像点获取的高程计算前后视校正影像，计算MI查找表。当前层每个点的初始高程由上一层匹配的高程得到，每个点的搜索范围由该点邻域(邻域大小取决于局部高差)内的高程确定，并留有一定余量。进行基于MI和CENSUS的半全局优化匹配，并根据平滑参数的设置和多路径独立匹配结果和聚合后的匹配结果间的差异性等多种方法确定可疑区域。根据正视影像和前后视校正影像在误匹配区域的纹理特征、误匹配区域周围高程和误匹配区域的匹配测度等对误匹配区域进行分类，并对不同类型的可疑区域进行精化匹配，精化匹配时MI查找表计算不使用可疑区域的影像。
(6) 生成DSM。底层匹配完成后，根据正视影像像点、像点对应高程、RPC计算地面点，根据设置的DSM网格间距通过距离倒数加权法生成DSM。
-
物方半全局匹配与像方半全局匹配的主要区别在于前者通过地面高程确定匹配像点间的对应关系，后者通过视差直接确定。因此在确定对应像点的时候物方方法计算量要大一些，本文采用一种局部拟合方法，可以较快地实现正视像点和前后视像点间的对应。本文提出的算法会存储每个像点对应的高程，而不是每层都生成DSM，因此可以快速计算出像点对应的地面坐标，并通过前后视影像的RPC参数计算对应像点。由于需要在当前高程的基础上在一定范围内搜索匹配像点，因此需要计算每个正视像点对应不同高程在前后视影像上的对应点。半全局优化需要较大的存储，因此需要分块进行匹配。本文首先计算当前块影像四角的初始高程h0、高程最大向上偏移量hu、高程最大向下偏移量hd各自对应的前后视影像坐标h0(xb,yb)、hu(xb,yb)、hd(xb,yb)(b表示后视，前视同)。然后通过最小二乘建立正视平面坐标(xn,yn)，搜索高程相对于初始高程的偏移量ho和偏移高程对应前后视坐标(xo,yo)与初始高程对应坐标h0(xb,yb)偏移量间的对应关系。
对于卫星影像来说，地形起伏较大的区域可能只存在于局部，如果计算影像范围内最大最小高程区间内的所有匹配代价将会造成计算资源的极大浪费。因此本文只计算上层匹配结果确定的有效搜索区域内的匹配代价。这样会造成相邻像元可能只有部分高程的匹配代价是对应的，但这并不影响半全局优化的实施。半全局优化也同样存在首末两个视差偏移，一个像素无对应传播者的情况，设 dmin和dmax分别代表最小、最大搜索位置，可以将首末超出搜索范围的dmin－1和dmax+1赋值为dmin和dmax即可，这时计算最小匹配代价时会出现两个相同的匹配代价，但不影响算法的正确性。同样，在多个高程的传播方向上没有对应高程时，若首尾之外的高程距离首尾只有一个搜索步距，设Lmin和Lmax分别为dmin和dmax对应的传播后的匹配代价，则计算最小匹配代价时只计算Lmin+P1或Lmax+P1及L0min+P2中的最小值(L0min为传播来向的最小匹配代价)。超出搜索范围的匹配代价统一加P2。其余有对应值部分同原半全局算法实施方法相同。
由于物方方案的高程初值是浮点型的数值，因此需要做一次对齐，设搜索步距为Δh，根据邻域高程确定的最小高程为hmin，则hmin应变换为int(hmin/Δh)×Δh。另外预先对整幅影像做CENSUS变换会比实时进行CENSUS变换减少很多计算量。但由于投影坐标为浮点数，因此可能存在半个像素以内的坐标误差。为了减小误差，可以在对前后视影像进行CENSUS变换时，除计算整像素处的CENSUS，还要计算半像素处的CENSUS。该变换只在近似核线方向进行，对于资源三号影像来说就是行方向。采用9×7窗口计算CENSUS并用两个整数存储结果，这时前后视的CENSUS对应数据量是影像数据量(16 bit)的8倍，正视影像不计算半像素处的CENSUS为4倍。加上MI和CENSUS一起使用及匹配后高程坐标的抛物线拟合基本可以忽略微小的坐标误差引起的匹配误差。
-
16个传播方向的设置方法如图 2(a)所示，由于第0~7个方向与目标像素相邻，高程出现较大变化的可能性低于第8~15个方向，因此第8~15个方向的P1、P2较第0~7个方向的P1、P2应减小，本文将8~15个方向的P1、P2均乘以0.5。
图 2 16个传播方向与平滑参数调整示意图
Figure 2. 16 Propagation Direction and Smooth Parameter Adjustment
图 3 动态调整平滑参数对“尖刺”的抑制效果
Figure 3. Effect of Dynamic Smooth Parameter on Peak Inhibition
如图 2(b)所示，左边三行表示传播来向的搜索空间，右边5行表示被传播方向的搜索空间。黑色方格表示最小匹配代价对应高程，本文根据当前高程与传播来向最小匹配代价对应高程的距离，动态调整P2。距离按步距个数计算，设两者高程相距n个步距，则P2要乘以n，n越大越容易抑制较大的高程变化，本文将0～7个方向的n上限设置为5，第8～15个方向的n上限设置为3，可以做到既保留了合理的高程突变，又抑制了“尖刺”。
-
误匹配形成的原因主要包括第1节中提到的5个原因及搜索范围过小导致的误匹配区域，搜索范围过小主要包括：高差较大的山区及所占影像面积较小但高度较大的建筑两类。前者在金字塔顶层匹配时解决，后者在底层精化匹配时解决。
形成误匹配的原因多种多样，表现在匹配代价上就是匹配代价最小的点并不是正确的匹配点，半全局优化一个突出的特点是容易形成误匹配的区域对平滑参数变化的敏感度较高，而正确匹配的区域对平滑参数的变化敏感度较低。因此可以通过变化平滑参数的设置并找出匹配结果发生明显变化的点作为可疑点。本文采用两步匹配法确定可疑区域：第一次采用动态调整平滑参数，第二次采用固定的平滑参数。但第二次使用的传播和聚合的匹配代价是第一次传播和聚合后的结果，但由于经过了多个方向的聚合，因此需要将传播后的匹配代价L除以一个系数，本文将该系数设为聚合方向的二分之一。
半全局匹配的另一个特点是正确匹配的区域多个聚合方向各自的最小匹配代价对应高程多数与聚合后得到的高程一致，而误匹配区域往往呈现多个聚合方向各自的最小匹配代价对应高程多数与聚合后得到的高程出现较大差异。本文根据多个聚合方向各自的最小匹配代价对应高程与聚合后最小匹配代价对应高程间的距离相差超过一个步长的个数确定可疑区域，若个数超过总聚合方向的一半，则认为该点可疑。
将正视影像像点根据匹配高程投射到前后视影像上，可以得到与正视影像大小相同的校正影像。根据正视影像和两幅校正影像像点间的相似性决定可疑区域。
通过前面三种方法确定的有效点可能并不完全形成区域，而是密度较大的可疑点群，若要形成区域还要做进一步处理。方法是将误匹配点“裹挟”的没有被确定为可疑点的点在5×5邻域内向周围探测可疑点，若探测数目超过一半则认为该点可疑。然后通过连通域保留算法去掉面积较小的可疑区域。最后做一次膨胀-腐蚀操作，确定最终的可疑区域。
-
针对大面积纹理缺失导致的可疑区域，第一次匹配完后，将高程通过距离倒数加权法用16个方向探测到的有效匹配点插值，精化匹配时将所有搜索步距对应匹配代价设为相同值，在精化匹配时可以实现平滑过渡。
针对小区域纹理丰富、边缘清晰误匹配区域，若该区域沿原始影像行方向(近似核线方向)两侧高度变化较小，则认为是搜索范围过小导致的误匹配。将区域膨胀，计算该区域沿原始影像行方向最大跨度d。精化匹配时该区域搜索范围为d乘以一个系数，由于高大建筑等在影像上只占很小的区域，也可统一按最大搜索范围设置。
针对最小MI不在对角附近的误匹配，计算匹配代价时只采用CENSUS。
针对最小CENSUS过大区域(表现为窗口对应点与各自中心点的差值相反)，计算匹配代价时只采用MI。
-
实验采用山西太原地区高精度LiDAR(3 m分辨率)和DOM(初始分辨率为0.5 m，匹配控制点时将DOM降为2 m分辨率)为控制和检验数据进行了DSM提取和精度评价，实验数据为两景包含正视、前视、后视影像的资源三号L1级影像。首先在每景影像和DOM间通过对数极坐标匹配获取4个控制点，控制点高程在LiDAR上通过平面坐标双线性插值获取。通过控制点对原RPC的分子常数项和归一化经纬度对应的分子一次项共6个参数进行校正，从而实现匹配获取的DSM与LiDAR坐标基准的对应。DSM的精度通过一定误差范围内DSM与LiDAR对应点差值的均值、绝对值均值和均方根误差等指标来表示。由于第二幅影像有大片居民区，而LiDAR上的人工建筑已经滤除，因此将第二幅DSM的山区和平地分开进行精度统计。统计范围为-8 m～8 m内的误差，低于-8 m在误差图中用黑色表示，高于8 m用白色表示。
通过表 1可以看出，实验1中DSM1均值接近0，说明系统误差很小，相比之下实验2中DSM2存在较大的系统误差，误差直方图(图 4(a)、图 4(b))的峰值对比更加明显。虽然城市区域的高大建筑没有列入统计范围(黑色区域)，仍然会引起较大的统计误差。通过实验3、实验4中DSM2山地、平地分开统计结果来看，城市区域所在的平地系统误差明显大于山地。而山地同样存在一定的系统误差，这与图 5(b)整体偏淡蓝色对应，估计是定向参数导致的误差。图 5(a)中椭圆区域误差较大，该区域DSM要高于LiDAR。通过观察影像可以看出该区域位于山的阴面，有大面积的树木遮蔽。而山的另一侧树木很少，误差也较小(绿色和黄色)。实验5中DSM2椭圆区域去掉人工建筑后仍有一些区域误差较大，从影像上可以看出这些区域均位于纹理很弱的大区域农田，需要指出的是这些区域只用MI时要明显优于只用CENSUS，但只用MI整幅DSM统计精度要低于CENSUS，而仅用CENSUS统计精度低于MI+CENSUS。实验6、实验7将正视影像分辨率降低一倍后精度出现较为明显的下降，超出8 m统计范围的区域明显增多。见图 6~图 9。实验8~实验15将搜索范围余量(括号内为上下各留的余量)与对应精度进行对比，可以看出余量为2时精度有一定下降，大于4时精度变化不大。同时对动态调整平滑参数策略和固定平滑参数策略下的DSM精度进行对比，结果表明前者精度高于后者。图 10和图 11表明本文提出的可疑区域探测方法可以较准确地探测出光照很弱引起的假纹理区域和局部搜索范围不够(主要是高大建筑)引起的匹配误差。通过精化，前者实现了可疑区域的平滑过渡，后者的高大建筑也可匹配出来(由于LiDAR上没有这些地物，因此不做精度统计)。
表 1 不同实验区域/条件下的误差统计结果(实验8～15括号内的数字为上下匹配决定的搜索范围基础上加的余量)
Table 1. Error Statistics under Different Experimental Areas and Conditions
序号 实验区域/条件 均值/m 绝对值均值/m 均方根误差/m 1 影像1(DSM1) 0.063 2.237 2.842 2 影像2(DSM2) -1.440 2.291 2.892 3 DSM2山地 -1.080 2.219 2.807 4 DSM2平地 -1.907 2.396 3.008 5 DSM2无人工建筑 -1.801 2.005 2.424 6 DSM1半分辨率 0.240 3.096 3.738 7 DSM2半分辨率 -0.927 2.632 3.280 8 DSM1矩形区域动态参数(7) 0.388 2.113 2.683 9 DSM1矩形区域固定参数(7) 0.396 2.141 2.718 10 动态参数(2) 0.534 2.194 2.775 11 固定参数(2) 0.605 2.240 2.830 12 动态参数(4) 0.381 2.157 2.732 13 固定参数(4) 0.400 2.158 2.736 14 动态参数(15) 0.382 2.140 2.713 15 固定参数(15) 0.605 2.240 2.830 
图 4 不同实验条件下的归一化误差直方图
Figure 4. Error Histogram in Different Tests
图 5 DSM相对LiDAR偏差分布
Figure 5. Error Image Between DSM and LiDAR for Image 1 and Image 2
图 6 DSM2山地、平地及去掉人工筑的统计区域
Figure 6. Different Test Zone for DSM2
图 7 纹理缺失平坦区域DSM误差分布图
Figure 7. Error Image for Texture-less Areas
图 8 纹理缺失区域CENSUS和MI结果对比
Figure 8. Comparison Between CENSUS and MI for Texture-less Areas
图 9 树木遮蔽区域DSM与LiDAR对比
Figure 9. Comparison Between DSM and LiDAR for Tree Obstructed Areas
图 10 对高大建筑的误匹配探测与精化
Figure 10. Detection Error for High Buildings and Refinment
图 11 对虚假纹理导致的误匹配区域探测与平滑过渡
Figure 11. Detection Error for False Texture and Smooth Transition
-
通过本文提出的物方半全局优化匹配方法实现了较高精度的资源三号三视影像DSM的提取。通过动态调整平滑参数可以抑制不合理的高程突变，且整体精度较固定平滑参数高。通过上一层的匹配结果并加上少量余量确定当前层的搜索范围可以明显减小计算量，同时可以保证精度。通过可疑区域检测可以较准确的定位虚假纹理和局部搜索范围过小导致的误匹配并在精化匹配时得到明显改善。
ZY-3 DSM Generation Method Based on Semi-global Optimization
-
摘要: 提出了一种基于物方半全局优化的资源三号卫星影像数字表面模型（digital surface model，DSM）提取方法，将半全局优化策略和影像金字塔结合，利用上一层金字塔的匹配结果来动态确定下一层每个像点的搜索范围，打破了半全局匹配需要建立匹配代价立方体的局限；以互信息和影像局部结构描述算子（CENSUS）距离加权和作为匹配代价函数，结合可疑区域探测技术，实现了高精度的DSM提取。通过实验对影响DSM精度的因素进行了总结分析。Abstract: This paper details a DSM (digital surface model) generation method using ZY-3 images based on object semi-global optimization. This method avoids the limitations imposed when building a match cost cube in the standard method. The proposed method combines semi-global optimization and an image pyramid to dynamically determine the search space of every pixel in the next pyramid layer according to the match result of the previous layer. Combining outlier detection technology and using mutual information and CENSUS as cost function, it realizes high-precision DSM generation from ZY3 images. We analyzed the key factors affecting DSM accuracy through experiments.
-
Key words:
- semi-global optimization /
- image pyramid /
- mutual information /
- CENSUS /
- outlier /
- ZY-3 satellite /
- DSM
-
图 2 16个传播方向与平滑参数调整示意图
Figure 2. 16 Propagation Direction and Smooth Parameter Adjustment
图 3 动态调整平滑参数对“尖刺”的抑制效果
Figure 3. Effect of Dynamic Smooth Parameter on Peak Inhibition
图 5 DSM相对LiDAR偏差分布
Figure 5. Error Image Between DSM and LiDAR for Image 1 and Image 2
图 8 纹理缺失区域CENSUS和MI结果对比
Figure 8. Comparison Between CENSUS and MI for Texture-less Areas
图 9 树木遮蔽区域DSM与LiDAR对比
Figure 9. Comparison Between DSM and LiDAR for Tree Obstructed Areas
图 11 对虚假纹理导致的误匹配区域探测与平滑过渡
Figure 11. Detection Error for False Texture and Smooth Transition
表 1 不同实验区域/条件下的误差统计结果(实验8～15括号内的数字为上下匹配决定的搜索范围基础上加的余量)
Table 1. Error Statistics under Different Experimental Areas and Conditions
序号 实验区域/条件 均值/m 绝对值均值/m 均方根误差/m 1 影像1(DSM1) 0.063 2.237 2.842 2 影像2(DSM2) -1.440 2.291 2.892 3 DSM2山地 -1.080 2.219 2.807 4 DSM2平地 -1.907 2.396 3.008 5 DSM2无人工建筑 -1.801 2.005 2.424 6 DSM1半分辨率 0.240 3.096 3.738 7 DSM2半分辨率 -0.927 2.632 3.280 8 DSM1矩形区域动态参数(7) 0.388 2.113 2.683 9 DSM1矩形区域固定参数(7) 0.396 2.141 2.718 10 动态参数(2) 0.534 2.194 2.775 11 固定参数(2) 0.605 2.240 2.830 12 动态参数(4) 0.381 2.157 2.732 13 固定参数(4) 0.400 2.158 2.736 14 动态参数(15) 0.382 2.140 2.713 15 固定参数(15) 0.605 2.240 2.830 
下载: 导出CSV
-
[1] Scharstein D, Szeliski R. A Taxonomy and Evaluation of Dense Two-frame Stereo Correspondence Algorithms[J]. International Journal of Computer Vision, 2002, 47(1):7-42 http://cn.bing.com/academic/profile?id=2177253483&encoded=0&v=paper_preview&mkt=zh-cn [2] Hirschmüller H. Accurate and Efficient Stereo Processing by Semi-global Matching and Mutual Information[C]. IEEE Conference on Computer Vision and Patten Recognition, Beijing, 2005 [3] Hirschmüller, H. Stereo Processing by Semi-global Matching and Mutual Information[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(2):328-341 doi: 10.1109/TPAMI.2007.1166 [4] Hirschmüller H, Daniel Scharstein, Evaluation of Stereo Matching Costs on Images with Radiometric Differences[J], IEEE Transactions on Pattern Analysis and Machine Intelligence,2009,31(9):1582-1599 doi: 10.1109/TPAMI.2008.221 [5] Gehrke S, Morin K. Semi-global Matching:An Alternative to LiDAR for DSM Generation[C]. IEEE Conference on CVPR, Miami, 2009 [6] Pablo d'Angelo, Peter Reinartz. Semiglobal Matching Results on The ISPRS Stereo Matching Benchmark[C].The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Hannover, 2011 [7] 张力,张继贤.基于多基线影像匹配的高分辨率遥感影像DEM的自动生成,武汉大学学报·信息科学版,2008,33(9):943-946 http://ch.whu.edu.cn/CN/abstract/abstract1689.shtml Zhang Li, Zhang Jixian. Automatic DEM Generation from High-resolution Satellite Imagery Based on Multiple-baseline Image Matching[J]. Geomatics and Information Science of Wuhan University,2008,33(9):943-946 http://ch.whu.edu.cn/CN/abstract/abstract1689.shtml [8] Zabih R, Woodfill J. Non-parametric Local Transforms for Computing Visual Correspondence[C]. Proceedings of the Third European Conference-volume Ⅱ on Computer Vision, London, 1994 -
点击查看大图
计量
- 文章访问数: 1406
- HTML全文浏览量: 42
- PDF下载量: 572
- 被引次数: 0
