在2016年之前,全世界最大的射电望远镜是位于波多黎各的阿雷西博射电望远镜(Arecibo radio telescope,Arecibo),于1963年11月1日落成,是美国利用天然海滩岩上的溶蚀洼地地形作为台址建成的,望远镜反射面口径305 m。2016年9月25日,中国利用贵州大型峰丛洼地——大窝凼的地形,建造了全球最大的单口径射电望远镜——500米口径球面射电望远镜(Five-hundred-meter Aperture Spherical radio Telescope,FAST),被称为“中国天眼”。由于FAST望远镜的接收面积比Arecibo射电望远镜提高了2.8倍,望远镜的灵敏度得到了显著的提升^[1]。当射电天文学家在追求更大的望远镜接收面积时,工程结构和力学支撑条件的限制使得全可动望远镜的口径只能到100 m左右; 若要追求更大的接收面积,射电望远镜只有“平躺”在地面、依附地形,以牺牲观察天区换取更大的接收面积。所以,形式上Arecibo和FAST都属于“平躺”式射电望远镜或“仰卧”式射电望远镜。Arecibo和FAST这2个超大的单口径射电望远镜都是利用圆形口径在300 m以上的洼地作为台址^[2]。FAST望远镜选址自1994年开始,至2007年7月获得立项批复,经历了13 a的艰难历程。FAST选址包括洼地地形条件、工程地质基础、自然灾害、无线电环境、社会经济条件等方面的分析。洼地的地形条件分析又包括收集各种地形图、卫星图像等资料,解读、量算、数字化、分析并建立洼地属性数据库,反复形态比较、绘制各种各样的地形剖面、拟合不同口径反射面的填挖方量,6个馈源支撑塔的地面分布优化等。找到大小和形态适宜的洼地,通过填挖方拟合和望远镜结构位置优化,确定FAST开挖中心、开挖标高,推荐6塔分布。堪舆台址,“给FAST找个家”并不容易^[3⇓-5]。开展省级范围FAST型大射电望远镜地形选址,包括应用30 m分辨率的数字高程模型(digital elevation model,DEM)搜寻大型洼地,研究洼地地形参数的适宜性。应用精细的洼地地形数据开展填挖方拟合和结构位置优化是选址第一阶段重要的基础工作,也是评价台址适宜性的关键^[6]。

FAST和Arecibo望远镜虽然都选址于峰丛岩溶洼地,但所处局地的海拔、相对高差和基岩背景等都不同(图1和表1)。FAST 望远镜台址大窝凼洼地大而深,东西宽990 m,南北长 1 050 m,深 360 m,最低垭口至洼地底部高差 140 m。Arecibo 望远镜台址洼地小而浅,东西长 460 m,南北宽 440 m,深 120 m,最低垭口至洼地底部高差70 m。大窝凼洼地的地质构造稳定,在小型断层构造交汇的背景下,SN向断层构造附近灰岩岩体破裂,两侧近水平的致密中厚层灰岩向外倾斜。若以东侧相距 1.2 km 的水淹凼洼地底作为水力基准面,高差达 104 m,水力梯度达到 86.6‰,这给洼地的溶蚀、侵蚀提供了较好的水动力条件。Arecibo 望远镜台址则多地震,珊瑚礁岩或海滩岩基岩松软,年降水量 1 000 mm,风化层 (土层)厚,植被生长茂盛。洼地底部与最近的低高程点(南偏东伏流出口) 相距 550 m,高差 50 m,水力梯度达到90.91‰。较大的水力梯度和有利的地质构造易于形成大型岩溶洼地。

FAST 望远镜的反射面是截取 600 m 直径球面上500 m 开口的球冠面,通过微小变形成抛物面,从而实现来源于宇宙射电信号的聚焦。通过平均分布于 600 m 圆周上的6个支撑塔悬挂馈源仓,在仓索的驱动下,使馈源仓精确移动到聚焦,实现宇宙射电信号的接收。FAST 望远镜反射面、馈源运移面和馈源支撑塔的几何参数如图2。相对于球冠面(反射面)底,500 m 口径的球冠深134.17 m,40°观察天顶角馈源运移的最高点176.38 m,分布于 600 m 圆周上的6个馈源支撑塔塔高 271.00 m^[8]。在理想状况下,最适宜 500 m 口径的 FAST 型大射电望远镜的洼地地形是圆形、 直径在 600 m 左右的洼地,底部平缓,坡面越高越陡; 洼地分水岭相对于洼地底高出270 m; 最低垭口相对于洼地底高出176 m,洼地底部134 m高度范围地形越圆、 越接近球面越好。这样不仅可以大大节约工程建设中的岩土开挖成本、减少反射面支撑和馈源支撑塔的建设,还使得望远镜借助洼地地形的阻挡,最大程度地减少外来电磁波的干扰。

DEM是一定区域位置(X, Y)及其相对应高程(H)属性的数据集。平面位置用格网表示,格网大小即是像素的大小,亦即平面分辨率。高度用高程值表示,以置信值为高程精度。单个 DEM 格网单元(像元)只能表示这个像元位置的高程值,所以 DEM 栅格数据也可以理解成具有三维空间属性(X,Y,Z)的连续点集。只有以连续的线、面分析 DEM 数据时,才能认识地形的高程在空间的分布差异和特征,认识区域地貌形态。通过模型或模式对 DEM 统计计算,提取高程点、等高线、分水岭、沟谷线、坡度、坡向、 坡长、坡度变化率、流域等地形参数信息。建立或生成 DEM 的方法有地面测量、 地形图等高线数字化、航空或航天立体测绘、激光高程扫描、声波-超声波(水下)等多种方法。

数字表面模型(digital surface model,DSM)是指包含了地表建筑物、桥梁和森林等高度的表面高程模型。随着航空、航天测绘技术的发展,全球测绘能力大大提升,但固体地表成片、成规模的地物如建筑物、桥梁、森林、冰川、常年积雪的高度无法确定并去除,所以获得的高程信息(DSM)中包含了 DEM 和 DEM 之上地物高度的信息。目前,网络共享的 DEM 都是应用地表类型修正 DSM, 从而获得反映固体地球地面高程的数据。

数字地形模型(digital terrain model,DTM)是基于 DEM 得到的单项数字地形参数,虽然也有空间差异和相同的格网,但没有了高程的属性,而是基于 DEM 的地形形态参数,如坡度、坡向、粗糙度等。

航空和卫星立体测绘、InSAR 技术测量可生产大区域的DEM或DSM,通常把DEM和DSM都统称为DEM。网站共享常用DEM的名称、生成技术、精度和覆盖范围等信息如表2。表中所列数字地形数据的水平分辨率,ALOS AW3D30_V2.2数据为12.5 m、SRTM3V00 3 数据的为 90 m,其他均为 30 m。NASADEM_HGTV001、SRTM3V003 的数据格式分别为 hgt 和 img,其他的均为 tif 格式。

图3(a)—(e)是表2 中所列 5 种数字表面数据 FAST 台址区域高程值拉伸后的灰度影像、派生的等值线图。十分明显的是 SRTM3_V003 的 DEM 灰度影像图的分辨能力不足,洼地地形反映不清晰,有明显的马赛克; NASADEM_HG 的 DEM 灰度影像图中亮点噪声较多,噪声反映的是局地地形直上直下,与实际不符; ALOS_AW3D30 的DEM 灰度影像图中亮点噪声和暗点噪声有所减少,影像也显得更加平滑; ALOS_PALSAR 的 DEM 灰度影像图有一些亮点噪声; ASTER_GDEMV3 的 DEM 灰度影像图虽然比 ALOS_AW3D30及ALOS_PALSAR 的在高亮值区范围要大些,但在地形的连续性方面更为合理,对低高程(洼地)的表达更明显。应用 ArcMap10.8 的表面插值功能,对各样例数据采用 20 m 等间距插值生成等高线,结果如图3(f)—(j)。 ALOS_PALSAR,ALOS_AW3D30及NASADEM_HG 的 DEM 等值线在较陡的区域有粘连,ASTER_GDEMV3 的 DEM 等值线显得相对平滑而少粘连。

通过上述 5 种 DEM 数据的 DEM 灰度影像和表面插值等值线的比较,在大范围的区域,开展岩溶洼地的搜寻以 ASTER_GDEMV3 的 DEM 数据为宜。

图4 是应用 DEM 坡度值的单色图像与地形高程的晕渲彩色图,色彩纹理和清晰度与 DEM 的分辨率、质量直接相关,分辨率越高晕渲彩色图反映的地形信息越精细、越清晰,效果越好。

图4(a)—(c)虽然分辨率都是30 m,但清晰度和效果差别较大。其中,图4(c)是最新发布数据NASADEM_HG 的晕渲彩色图,相对于图4(a)和(b),图面反映的地形更清晰,颜色反差和地形效果更好,只是局部噪声点比较大。相比较而言,开展岩溶洼地的搜寻以ASTER_GDEMV3 的 DEM 数据为宜。在实际选用 DEM 数据时,不仅要分析数据附带的精度和质量评价资料,更要选择地形变化剧烈、高差大的区域检测和对比分析数据的真实情况,选用适宜的数据源。

开展 FAST 型大射电望远镜选址洼地地形搜寻包括应用30 m辨率的 ASTER_GDEMV3 数据开展大区域的洼地搜寻、洼地形态参数提取,应用大比例尺地形图等高线矢量化的洼地填挖拟合2个方面的研究。ArcGIS10.8 的软件处理平台对具有统一投影坐标系统的点、线、面的矢量数据和栅格数据可开展管理、处理、分析和显示等。ArcMap Model Builder 可有序组织ArcMap10.8 附带的处理工具,实现多个处理过程的一体化、“一键完成”。因此,通过构建洼地地形参数提取模块、填挖拟合模块实现洼地地形参数“一键一个要素”的提取和分析。

DEM 的本质是高程数值的矩阵,洼地地形参数中最低点、山峰点、垭口点、分水线、洼地等深线的地貌部位形态、DEM 高程和处理提取方法如表3。

在 ArcMap10.8 的软件处理平台下,通过栅格和矢量的组合处理即可实现上述基本地形参数的提取。地形参数提取模块的主要环节和步骤如表4。这些参数提取模块已在贵州全省范围的 DEM 数据洼地搜寻处理分析中实际应用。由于存在 DEM 栅格的分辨率大小和提取采用的背景参数设置的差别,提取结果应该在 DEM 栅格图像背景上筛选判断和确认。

图5 是大窝凼周区洼地最低点、山峰点、分水线、洼地等深线图。分析发现: ①大窝凼洼地不管是洼地底高程、底部等值线,还是大窝凼底部地形态势,都与1∶1万地形图形态十分接近,只是局部小规模形态被“平滑”掉了; ②每个完整的分水线封闭区都有洼地最低点、洼地等深线分布于地形高程低值区; ③山峰点都较好地分布在地形高程高值区的分水线上或附近。由于最低点、山峰点是基于像素的计算,是某个或某几个像素的中心点, 而分水线虽然是基于像素计算得到,但经过后期矢量线的锯齿平滑等处理,并不都经过像素的中心,所以大多数山峰点只能在分水线附近。分水线高程采样点高值点及高程鞍部点及高程图分水线上的局部低点一定是垭口点,而局部高点可以是山峰点也可以是洼地一侧宽平山梁上的一个分水点,分水线上的垭口低值点与高值点(山峰或高分水点)一定是相间分布。图5只是提取洼地主要的山峰点及其高程,而不能把分水线上的高值点全部提取。图6 是采用30 m间距点提取的分水线高程值,在结果比较、判断分析、分类后得到的大窝凼洼地垭口(低值)点与山峰 (高值)点分布图。图7 是大窝凼洼地 30 m 间距点分水线高程值的展开图,可十分清晰地看到垭口(低值)点与山峰(高值)点的分布特征,垭口与山峰的形态及高程一目了然。参考图6,可分析垭口与山峰的分布方位。

搜寻到了一个深度和口径基本合适的洼地,还需对洼地的地形开展三维的分析评价。平面形态、剖面形态千差万别的洼地,只有通过拟合分析才能显示口径 500 m、球径 600 m 球冠反射面与地形的“交-切-割”关系。工程建设需要开挖土石方少,同时反射面下方空间小,从而既可以减少开挖成本、降低反射面下部支撑结构建设,又能降低造价,同时也可以提升望远镜的性能。研究应用洼地区 1∶1 万地形图数字化等高线,经插值和重采样,生成了1 m分辨率的 DEM。借助 ArcGIS10.8 的三维空间分析工具,进行填挖方量拟合,平衡开挖土石方量和反射面下方空间,寻找最优球冠安置的平面位置和高度位置。

ArcMap Model Builder 的填挖方量拟合模型包括下列几个主要步骤:

1)创建一个与 DEM 数据范围相同、分辨率相同的 0 值高程栅格文件;

2)设定球冠最低点三维位置(X₀,Y₀,Z₀)、球冠直径(D₀)、球冠开口 (D₁),应用球冠参数和几何位置计算并将球冠区域的栅格 0 值改为球冠面高程,得到球冠面 DEM;

3)应用球冠区的多边形矢量裁剪得到地形 DEM 栅格表面和球冠面 DEM 栅格表面文件;

4)应用[填挖方]工具,分别输入填挖之前的地形 DEM 栅格表面和填挖之后的球冠面 DEM 栅格表面文件,计算得到开挖区 DEM 差值图像和填方区 DEM 差值图像,同时得到开挖(填方)总方量、开挖(填方)区的垂直投影面积、开挖(填方)的表面面积。通过多次拟合计算,即可回归分析得到填挖方优选的球冠位置,并据此综合分析得到经济最优的开挖方案。

在 ArcScene 的环境下,可以应用地形、球冠面、馈源塔、填挖区域、剖面线等具有三维信息的表面栅格或 TIN 进行组织和参数设置,显示各要素间的三维关系。

图8(a)是某洼填挖方量拟合得到的结果图,蓝色为挖方区,红色为填方区域。应用开挖区 DEM 差值图像和填方区 DEM 差值图像,进行等值线处理,即可得到开挖深度等值线图和填方深度等值线图,如图8(b)。

图8显示了A-A'和B-B'这2个剖面的位置,图9 则显示了洼地的地形-球冠叠合剖面图,应用剖面显示开方和填方间的“交-切-割”关系。图10是某洼地 DEM 分析数据在 ArcScene 的环境下,三维显示地形表面、 开挖面、反射面、地形剖面、支撑塔的空间关系。可以交互变化观察视角和位置,直观分析各要素间的空间关系。

通过研究FAST型大射电望远镜选址中的洼地搜寻、地形分析、填挖拟合和三维分析,实现了基于 DEM 数据的大范围、快速、准确的处理分析。ArcGIS10.8软件平台不仅可以管理和分析各类空间数据,而且可以一体化快速分析提取洼地形态参数、拟合开挖、剖面叠合分析和三维显示。解决了省级区域大型岩溶洼地地形搜寻中的关键技术。

1)DEM 数据选用要根据工作需要、工作区域大小和研究条件来选定,适宜即可; 全省范围的大型岩溶洼地的搜寻宜应用 30 m 分辨率的 DEM; 洼地填挖方量拟合、剖面叠合、三维空间交互分析宜采用 1 m 分辨率由 1∶1 万地形图等高线生成的 DEM。

2)应用 ArcGIS 强大的空间数据管理、处理和分析功能,通过构建洼地地形参数提取模块,实现贵州全省洼地地形提取、识别和洼地地形参数“一键一个要素”,达到快速、准确。大量的空间图形计算分析需有优良的图形处理硬件环境,执行模块处理时要减少并行处理,以免经常报错。

3)填挖方量拟合模型以栅格为基础,受到栅格大小和高程精度的限制, 开挖区和填方区边缘区锯齿效应明显,生成等值线时要进行甄别处理。

4)由于岩溶洼地地形类型多样、地形高差变化较大,DEM 数据也存在局部质量和精度问题,参数提取、空间分析结果要交叉验证和检查,以减少虚假信息干扰。