ENVI5中“资源三号”立体像对DEM提取

1.概述    

      数字高程模型(Digital Elevation Model)，简称DEM，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。DEM除了包括地面高程信息外，还可以派生地貌特性，包括坡度、坡向等，还可以计算地形特征参数，包括山峰、山脊、平原、位面、河道和沟谷等。建立DEM的方法有多种，从数据源及采集方式主要有：根据航空或航天影像，通过摄影测量途径获；野外测量或者从现有地形图上采集高程点或者等高线，后通过内插生成DEM等方法。

    ENVI 的DEM自动提取模块能够简单、快速地从扫描或者数字航空影像、摆扫式或推扫式卫星传感器，如ALOS PRISM, ASTER, CARTOSAT-1, FORMOSAT-2, GeoEye-1, IKONOS, KOMPSAT-2, OrbView-3, QuickBird, WorldView-1/2、SPOT1~6，以及国产的资源三号、天绘卫星系列等创建DEM。独特的向导工具使得DEM提取更加简便、自动；提取完成后还可以使用DEM编辑工具对DEM进行局部编辑；DEM模块全面支持RPC模型参数，尽可能用最少的地面控制点达到有效精度；重叠区可以选择基于灰度或者特征自动寻找匹配点，减少工作量。

     资源三号搭载了四台光学相机，包括一台地面分辨率2.1m的正视全色TDI CCD相机、地面分辨率3.5m的前视和后视全色TDI CCD相机、一台地面分辨率5.8m的正视多光谱相机。其中前正后视全色相机，推扫成像形成三线阵立体像对，可用于DEM提取。经过测试发现，前视和正视，或后视和正视，可以组成立体像对进行DEM提取，效果较好，优于前视和后视提取的DEM。

2.详细操作步骤

      DEM Extraction工具分为9个步骤，总体上我们可分为6个步骤，如下图所示。

图1.1 DEM提取流程

 

第一步：输入立体像对

（1）单击File->Open，以正视和后视图像为例。打开"ZY3_01a_hsnnavp_*"和"ZY3_01a_hsnbavp_*"文件夹内的tif文件。

（2）在Toolbox中，双击/ Terrain/DEM Extraction/DEM Extraction Wizard: New。

（3）单击Select Stereo Image按钮，选择正视图像（NAVP）为Left Image，后视图像（BAVP）为Right Image。

注：1、ENVI一般能自动识别RPC文件，否则需要手动选择RPC文件。

2、左影像选择垂直观测的影像，或者观测角度小的影像。可以简单通过对比立体像对两幅影像的地面分辨率，分辨率高的当作左影像。

（4）单击Next按钮。

图2.1 Step 1 of 9输入立体像对

第二步：定义地面控制点

      提供三种定义地面控制点方式：不定义（No GCPs (relative DEM values only)）、交互式定义（Define GCPs Interactively）和读取控制点文件（Read GCPs From File）。由于缺少地面控制点数据，这里选择不定义控制点，单击Next按钮。

注：不定义地面控制点得到的DEM是相对高程，即以卫星地势面为基准。

图2.2 Step2 of 9定义地面控制点

第三步：定义连接点

（1）提供三种定义连接点方式：自动寻找、交互式手工定义和外部读取控制点文件。选择选择自动寻找，参数设置如下所示，

• 连接点数目（Number of Tie Points）：1000。
• 搜索窗口大小（Search Windows Size）：481。
• 移动窗口大小（Moving Window Size）：101。
• 最小相关系数（Minimum Correlation）：0.90。
• 平均高程（Region Elevation）：自动从影像读取。
• 是否检查连接点（Examine and Edit Tie Points）：Yes。

图2.3 Step 4 of 9自动选择Tie点参数设置

注：几个参数的设置说明如下说明：

• 搜索窗口大小（Search Windows Size）：大于或等于21的任意整数，并且必须比移动窗口大。该参数的值越大，找到匹配点的可能性也越大，但同时也要耗费更多的计算时间。大致确定搜索窗口大小的方法：在立体像对（带有粗略地理坐标）两个图像上找到一个同名点，量测两个图像上同名点的距离D（像素单位），搜索窗口最小设置为（D+1）*2。
• 移动窗口大小（Moving Window Size）：在搜索窗口中进行检查，寻找地形特征匹配的小区域。移动窗口大小必须是奇数。最小的移动窗口大小是5，即为5×5像素。使用较大的移动窗口将会获得更加可靠的匹配结果，但也需要更多的处理时间。移动窗口的大小跟图像空间分布率有关系，根据如下所列设置：
  • 大于等于10 米分辨率图像，设置值的范围是9~15。
  • 5~10 米分辨率图像，设置值的范围是11~21。
  • 1~5米分辨率图像，设置值的范围是15~41。
  • 小于1米分辨率图像, 设置值的范围是21~81 或者更高。

（2）单击Next按钮，进入Step 5/9步骤。

（3）本例子中，自动找到414个Tie点，Maximum Y Parallax（Y视差）=2.3394，小于10即可进行下一步处理。从点的分布来看，有部分地区没有Tie点，下面手动交互添加一些Tie点。

注：如果Maximum Y Parallax（Y视差）>10，单击Show Table，选择Sort Table By Error，误差大的点排在前面，逐个选择Tie点查看精度，将偏离较大的点进行微调或者直接删除。

图2.4 5 of 9自动选择Tie点及分布

1）单击"Params…"按钮修改两个窗口大小为481和101。

图 预测参数设置

2）在左边图像的无Tie点区域，定位一个点，单击Predict Right按钮，自动在右边图像上定位一个点，确认两个点是同名点，单击Add按钮添加。

注：一般预测的点都很准确，个别不准确可以由右边图像上的点来预测左边图像（Predict Left按钮）的同名点。

3）重复以上步骤，手动增加一些Tie点，尤其是地形起伏较大区域。

（4）当误差达到一定范围，如这里的3.4373（Maximum Y Parallax<10），单击Next按钮。

图2.5手动添加一些Tie点的分布

第四步：设定DEM提取参数

（1）Step 6 of 9，ENVI会生成核线图像，可以用于立体观测，选择一个路径分别输出，单击Next。

图2.6 Step 6 of 9 输出核共线图像

（2）在Step 7 of 9中，设置DEM输出投影参数。这里直接设置为UTM坐标，输出像元大小为10m（记得回车自动计算输出像元行列数，否则结果会带很多黑背景），单击Next。

注：这里输出的像元大小没有一个明确的规定，一般可选择 10或者15米输出。

（3）在Step 8 of 9中，设置DEM输出参数：

• 最小相关系数阈值（Minimum Correlation）：0.7。值的范围是0~1，这个参数是评价两个点匹配好坏。这个阈值越大，匹配精度越高，能得到的匹配点越少。一般设定在0.65~0.85之间。
• 背景值（Background Value）：-999。设定DEM的背景像素值。
• 外边界清理焊缝(Edge Trimming)：0。值的范围是0.0~0.6。设定输出DEM外边界清理焊缝宽度，用占整个DEM的百分比来标示。
• 移动窗口大小（Moving Window Size）：13x13。定义计算两图像相关性的范围大小，于来执行图像匹配，值越大越可靠，精确的匹配结果越少，计算量越大。
• 地形地貌（Terrain Relief）：Moderate。三个级别——Low、Moderate和High。Low用于覆盖区域地形平坦；Moderate用于大多数地形；High用于地形起伏较大的山区，地形地貌很明显区域。
• 地形细部（Terrain Detail）：Level 6。设置DEM地形细节等级。等级越高，生成的DEM越精细，处理时间越长。
• 数据输出类型（Output Data Type）：Integer。16bit的Integer和32bit的Floating Point。

（4）选择输出路径和文件名。

图2.7 Step7 of 9、8 of 9DEM输出投影与输出参数设置

第五步：输出DEM并检查结果

      在Step 9 of 9中，单击Load DEM Result to Display按钮加载生成DEM结果。

图2.8 DEM结果

第六步：编辑DEM

      生成的DEM有些地方如果有厚云等，对获得的地形有影响，可以通过手动编辑来更改这些区域的高程数据。

（1）在Step 9 of 9中，单击Load DEM Result to Display with Editing Tool按钮，打开DEM编辑工具。

（2）在Image 窗口中绘制需要编辑的区域，在Method中选择一种修改方法，单击Apply to Region of Interest或者Apply to Entire Band进行修改。几种修改方法说明如下表：

（3）单击Save按钮，保存为工程文件。

表 2.2 编辑DEM高程值的7种方法

方法	说明
Replace with value	用指定的值替换感兴趣区内的高程值，需要设定一个替代常量。
Replace with mean	用感兴趣区内原来的平均高程值替换整个感兴趣区内的高程值。
Smooth	对感兴趣区内做低通卷积滤波，需要设定一个卷积核，默认为3x3。
Median Filter	对感兴趣区内做中值卷积滤波，需要设定一个卷积核，默认为3x3。
Noise Removal	如果感兴趣区内原高程值大于其周围高程值的标准差，则用周围高程值的中值代替。
Triangulate	用三角内插算法对感兴趣区内的高程值重新插值。
Thin Plate Spline	用薄板样条插值算法对感兴趣区内的高程值重新插值。

   如下图是利用相同的流程和参数，使用其他资源三号数据提取得到的DEM结果。

图2.9资源三号立体像对提取的DEM结果

 

3.常见问题解答

3.1整体效果差

   问题描述：呈现大范围"三角形"状，地形不平滑

  解决方法：1、可能Tie点太少造成，增加Tie数量。

                   2、可能立体像对拍摄夹角太大，如选择前视和后视作为立体像对时。

3.2部分地区地形插值不平滑

   问题描述：如下图所示，主要产生的原因是地形起伏较大，而且是背着传感器拍摄方向，尤其是山体存在阴影情况。

    解决方法：增加Tie点。

图3.1 DEM结果局部不平滑

3.3处理速度慢

     设置地形地貌（Terrain Relief）和地形细部（Terrain Detail）级别太高，如设置Terrain Relief：High和Terrain Detail：Level 7。

    一般设置Terrain Relief：Moderate和Terrain Detail：Level 6即能满足大部分处理要求。

3.4内存溢出

    问题描述：提示错误类似Unable to allocate memory: to make array。

    解决方法：设置的分块大小太大导致。在ENVI Classic中，File – Preferences 面板中（如下图），建议设置Cache Size为内存大小的50%；Image Tile Size设置可参考：

• 内存 – 2G，Image Tile Size – 1M
• 内存 – 4G，Image Tile Size – 4~10M
• 内存 – 8G及以上，Image Tile Size – 50~100M

图3.2 ENVI Classic分块大小设置