1　引 言

中国遥感卫星辐射校正场敦煌场以敦煌市郊戈壁滩为场地，具有表面均匀、气候干燥、方向特性较好等得天独厚的特性，是国际上公认的有利于传感器在轨辐射定标的戈壁场^［1］。自2002年通过国家验收以来，敦煌场每年都为我国在轨运行遥感卫星星上仪器辐射校正提供支撑，是推动我国遥感技术及其应用定量化发展必不可少的基础设施之一^［2］。经多年建设发展，敦煌场地设施逐步完善，已成为国产对地观测遥感卫星（气象卫星、资源卫星、海洋卫星、环境减灾卫星、测绘地震及侦察卫星等）真实性检验与辐射定标的最重要场地^［3］，与此同时也为一些国外对地观测卫星服务，敦煌场对国内外在轨卫星辐射定标和真实性检验都具有重要意义。

地面光谱采样是获取像元尺度地表反射率的重要方式^［4-5］，其准确度会直接影响传感器定标结果和反射率产品真实性检验结果^［6］。采样点的位置分布会对地面光谱采样准确度产生影响，可以通过优化采样方案提高地面光谱采样准确度，令地面测量点尺度数据更高精度聚合到像元尺度^［7］。

已有像元尺度地面采样方法主要分为基于参考影像的采样方法与无需参考影像的采样方法两类。无需参考影像的采样方法是指不需要样区的任何先验数据，直接按照某种规律进行地面采样的方法，这类方法包括随机采样^［8］、系统采样^［9-11］、对角线采样^［12］等。在敦煌场获取陆地卫星像元尺度地表反射率时，中国资源卫星应用中心常用这类方法，在550 m测试区每隔50 m设置一个样方进行地面光谱测量^［13］。进行气象卫星像元尺度地面光谱采样时，国家卫星气象中心也常用这类方法，对11个均匀分布在10 km中心场的目标点进行观测^［14］。基于参考影像的采样方法是指需要样区历史影像作为先验数据，针对样区地表设计采样点布局的采样方法，这类方法中具有代表性的有分层采样^［15］、四叉树采样^［16］和模拟退火采样法^{［17，18］}。

当前采样方法虽多，但对采样方法准确度的研究开展较少。上述两类采样方法中哪种更适合敦煌场地、哪种能满足尺度转换和定标精度要求？敦煌场地现有针对陆地卫星采样方法需要对上百个点进行观测，是否有必要？针对气象卫星在10 km场地内仅设11个目标点是否足够？这些样区能否适当减少采样点数量，从而减小进出场地造成的破坏^［19］？若能实现对采样方法准确度的定量分析，这些问题就能够得到解决。

本文参考了其他地表参数和其他场地采样思想及陆地卫星、气象卫星现有定标方法，开展了不同采样方法、不同采样点数量对点面尺度转换影响的相关研究，利用高分辨率影像作为参考数据，定量评价不同采样方法、不同采样点数量引起的点面尺度转换误差，从而评价采样方法的准确度、确定不同像元尺度的最优采样方法。

2　场地与方法

2.1　敦煌场地

敦煌国家辐射校正场位于甘肃省敦煌市西约30 km处（40.04° N~40.28° N，94.17° E~94.5° E，海拔约为1200 m），为党河冲击扇叶戈壁，地势平坦，表层基本无植物生长，由多种砾石、砂及少量黏土组成，戈壁表面粒径分布比较均匀^［3，20］。且该校正场具有太阳辐射性强、光照充足、降水量少、蒸发强烈、能见度好等特点^［21］。

同步观测区面积往往根据卫星产品的分辨率设定，考虑临近像元效应通常将其设为3×3~5×5个像元大小。早期敦煌陆地卫星同步观测区面积常设为550 m×550 m^［13］、气象卫星同步观测区常设为10 km×10 km^{［14，22］}。随着卫星传感器技术的发展，卫星产品分辨率较早期已有很大提高，现有陆地卫星以30 m分辨率居多、气象卫星以1 km分辨率居多，为了降低采样难度并提高采样准确度，本实验组将观测区面积分别设为150 m×150 m（5×5个像元大小）和3 km×3 km（3×3个像元大小）。

分别划定两块150 m陆地卫星同步观测区与两块3 km气象卫星同步观测区，对陆地卫星与气象卫星像元尺度地表反射率采样方法开展研究。首先在敦煌原550 m陆地卫星同步观测区内划出150 m区域（即“150 m国家场”），在原10 km气象卫星同步观测区内划出3 km区域（即“3 km国家场”）；然后以高分一号2 m分辨率卫星影像作为参考数据、空间变异系数为评价指标，在整个敦煌校正场范围内分别择出一块150 m×150 m最均匀区（后文称“150 m新选场地”）和3 km×3 km最均匀区（后文称“3 km新选场地”）。上述150 m国家场、150 m新选场地、3 km国家场与3 km新选场地位置分布如图1所示，经纬度信息如表1所示。

图 1. 场地位置示意图

Fig. 1. Schematic diagram of site location

下载图片 查看所有图片

2.2　采样方法

在地面光谱采样工作中，测量人员常在场地内初步确定多个2 m×2 m的基本采样单元，在基本采样单元内用光谱仪进行多点测量。本实验组在进行采样方法研究时，分别对2 m基本采样单元在场地内的布设方案和2 m基本采样内光谱仪的测量方案开展研究。

确定样区内2 m基本采样单元的位置分布时，应使多个采样单元能更好地代表整个样区，在进行这部分采样方法研究时本实验组用到系统采样法与模拟退火采样法；确定光谱仪在2 m基本采样单元内的最佳测量方案时，要令多个测点的反射率数据更好地代表2 m采样单元反射率，对其开展研究时本实验组用到随机采样法与系统采样法。

2.2.1　系统采样

系统采样按照一定的空间间隔布设样点，令样点均匀分布在目标区域。2 m基本采样单元内，地面光谱测量点个数分别为1，5，9，16，25的系统采样方案如图2所示。像元面尺度内，基本采样单元个数为1，5，9，16，25的系统采样方案如图3所示。

图 2. 2 m基本采样单元内光谱测量点系统采样分布

Fig. 2. Sampling distribution of spectrometer system in 2 m basic sampling unit

下载图片 查看所有图片

图 3. 像元尺度内基本采样单元系统采样分布

Fig. 3. Sampling distribution of basic sampling unit system in pixel scale

下载图片 查看所有图片

2.2.2　模拟退火采样

1982年，Kirkpatrick等首次提出模拟退火算法用于解决大规模组合优化问题^［18］。该算法受到Metropolis等于1953年提出的Metropolis准则启发，基于物理中固体物质退火过程与一般组合优化问题之间的相似性，对固体退火过程的物理图像和统计性质进行模拟，从而求出组合优化问题的近似最优解^［23］。模拟退火算法在迭代过程中以一定概率接受较差解，可有效避免陷入局部极小并最终趋于全局最优，是解决组合优化问题的有效近似算法^［24］。

采样点布局优化问题属于组合优化问题，即求在整个样区具有最优代表性的采样点组合，可以利用模拟退火算法实现对各场地样点布局的优化^［25］。本实验组进行2 m面尺度到像元面尺度采样方法研究时用到模拟退火采样。将样区的多时相历史影像作为参考数据模拟真实地表，将每个像元位置视为一个可采样点，用模拟退火算法对参考影像进行计算，在采样点个数确定的情况下，求出令采样准确度最高的采样点布局^［17-18］。采样相对误差是评价采样方案准确度的重要指标，将模拟退火算法目标函数设为采样点在多景参考影像上的采样相对误差之和。单景影像上的采样相对误差ERS与目标函数fERS的表达式为

式中：mallpoints为参考影像所有像元均值；msample为采样点均值；N为参考影像数量；ERS,k为第k景参考影像的相对误差。

2.3　基于高分辨率影像的采样准确度分析方法

采样准确度是指采样结果与场地真值间的一致程度，相对误差与平均相对误差是评价准确度的常用指标^［26］，相对误差与平均相对误差值小，采样准确度高，反之采样准确度低。本实验组基于这两个指标用高分辨率影像模拟真实地表，对不同采样方案的准确度开展定量评估。优化采样是为了减小由采样点代表性有限引起的后续尺度转换过程的误差，由于敦煌场地较均匀，本实验组选择这类场地常用的尺度转换方法——算术平均法对采样数据进行升尺度转换^［27-28］。该方法用地面多个点测量的算术平均值作为像元尺度相对真值与待检验产品估算值进行对比，这种尺度转换方法的相对误差计算公式与采样相对误差的计算公式恰好一致，所以本实验中采样点的尺度转换相对误差相当于采样方案的相对误差。

通过高分辨率影像分析地面光谱采样准确度对影像幅宽和分辨率都有一定要求。地面光谱测量时，光谱仪视场直径约为0.44 m，而高分一号全色卫星影像数据分辨率为2 m，像元覆盖区域比光谱仪视场范围大得多，所以仅采用卫星影像进行不确定度分析所得结果缺乏可靠性；无人机（UAV）遥感技术可以快速获取区域的遥感图像，数据可以满足分辨率方面的要求，但无人机飞行高度有限，捕捉地面的范围有限^［29］，数据无法覆盖整个样区大小，所以不能仅通过无人机数据对整个像元尺度采样准确度进行分析。为了克服影像数据的限制从而实现对采样方法准确度的定量评估，可以将尺度转换分解为地面点尺度到2 m面尺度和2 m面尺度到像元面尺度两个步骤，地面点尺度到2 m面尺度的采样准确度分析需要分辨率远高于0.44 m、数据覆盖范围不小于2 m的影像，无人机数据可以充分满足该条件；而2 m面尺度到像元面尺度的采样准确度分析需要分辨率不低于2 m、覆盖范围达到千米级的影像数据，高分一号全色数据与该要求吻合。因此本实验组分别利用高分辨率无人机影像与卫星数据计算这两部分采样相对误差，衡量采样方法准确度。

2.3.1　地面点尺度到2 m面尺度采样准确度分析方法

地面点尺度指光谱仪捕捉地面点范围尺度，根据光谱仪视场角与测量时的高度推算，点尺度约为直径为0.44 m的圆。地面点尺度到2 m面尺度的采样准确度即地面测量准确度，对其进行定量计算需借助高分辨率无人机影像模拟真实地表，将2 m基本采样单元的无人机影像所有像元值均值作为该基本采样单元真值，在采样单元内选取多个地面光谱测量点并求其均值作为采样单元的地面测量值，计算测量值与真值相对误差，用于衡量点到2 m面尺度的采样准确度，过程如图4所示。

图 4. 地面测量准确度分析示意图

Fig. 4. Schematic diagram of ground measurement accuracy analysis

下载图片 查看所有图片

地面光谱测量过程在无人机影像上进行模拟，将模拟值作为地面光谱测量值：光谱仪视场角通常固定为25°，故根据光谱仪到地面距离可推算光谱仪测量范围，例如图5，光谱仪距地面1 m，推算出光谱仪测量范围是直径约为0.44 m的圆。在2 m范围无人机数据内选取像元作为地面测量点中心做圆，采用反距离加权法模拟出地面光谱测量的像元灰度值。设无人机数据分辨率为0.08 m，光谱仪测量范围直径为0.44 m，单点地面光谱测量模拟值计算范围如图6所示。将地面光谱测量的各个单点模拟值取均值作为2 m采样单元的地面光谱测量值。

图 5. 地面光谱测量示意图

Fig. 5. Schematic of ground spectrum measurement

下载图片 查看所有图片

图 6. 无人机数据模拟地面测量点尺度

Fig. 6. UAV data simulation ground measurement point scale

下载图片 查看所有图片

无人机数据有3个波段（RGB），根据上述步骤模拟计算得到的单个采样点在n波段的模拟值为VSV,n，将其取均值得到样点均值VSM,n，影像所有像元在n波段均值为VAM,n，此时单景影像上n波段的采样相对误差的表达式为

本实验组用多景无人机数据评估各个采样方案准确度，将多景影像在多个波段的平均采样相对误差（MSRE）作为采样方案的准确度衡量标准，采样方案的单波段平均相对误差MSREn与MSRE的表达式为

式中：ERE,n-p指第p景影像上第n波段的采样相对误差；P为影像数量。

2.3.2　2 m面尺度到像元面尺度采样准确度分析方法

像元面尺度分为陆地卫星像元尺度与气象卫星像元尺度，分别为150 m×150 m与3 km×3 km。2 m到像元面尺度的采样准确度借助2 m分辨率卫星数据定量评估。用高分辨率卫星影像模拟真实地表，在影像上进行采样模拟计算准确度，将影像上目标场地范围内所有2 m像元Digital Number（DN）值取平均作为场地真值；根据2 m采样单元坐标位置找出其在目标场地影像中对应的2 m像元位置，将所有2 m采样单元对应的像元DN值取平均作为目标场地观测值；求出观测值与真值的相对误差，作为2 m面尺度到目标场地尺度采样准确度评价指标。

进行这部分采样准确度计算时，本实验组采用成像日期在参考影像之后的多时相真值影像衡量不同采样方案准确度，将多景真值影像上的最大采样相对误差（EMRS）作为这部分采样准确度的评价指标。

3　结果分析

3.1　地面点尺度到2 m面尺度的采样准确度

为了分析不同采样方法的准确度，对比采样方法优劣并找到2 m采样单元内的最佳测量方案，本实验组将2020年9月18日在敦煌辐射校正场拍摄的11景高分辨率无人机数据裁剪成137幅2 m大小影像对不同采样方法进行分析，裁剪过程如图7所示。无人机数据含RGB三个波段，根据影像上靶标的大小估算其分辨率约为0.0035 m。

图 7. 无人机数据裁剪示意图

Fig. 7. Schematic diagram of UAV data clipping

下载图片 查看所有图片

将2 m基本采样单元内的光谱采样点个数设为1，5，9，16，25，分析随机采样与系统采样两种方法的准确度，比较两种方法优劣。在2 m的无人机影像上模拟随机采样过程中，需注意在距离2 m边界不小于0.22 m（光谱仪视场半径）的范围内随机选择视场中心位置，以保证光谱仪测量视场完全落在2 m采样单元内。计算得到采样方案在137幅无人机数据上的137组采样相对误差值（RGB三个波段，每个波段137个值），将各波段的多个相对误差取均值，得到采样方案在不同波段的采样平均相对误差，结果如图8所示。

图 8. 地面点尺度到2 m面尺度采样平均相对误差。（a）R波段；（b）G波段；（c）B波段；（d）RGB波段

Fig. 8. Sampling average relative error from ground point scale to 2 m surface scale. (a) R band; (b) G band; (c) B band; (d) RGB bands

下载图片 查看所有图片

地面点尺度到2 m面尺度采样相对误差结果表明，在采样点数相同情况下，系统采样准确度明显高于随机采样，在2 m基本采样单元内进行光谱测量应优先选用系统采样法。在RGB三个波段上，不同采样方法对采样相对误差的影响趋势大致相同，采样误差都随采样点个数增多而下降；同种采样方法在R波段的误差略大于G波段，G波段的误差略大于B波段。

系统样点个数为1时仅增设4个点至5个采样点，R波段的采样平均相对误差即可从1.183%降至0.563%、减小了0.62个百分点，G波段采样平均相对误差从1.041%降至0.507%、减小了0.534个百分点，B波段采样平均相对误差可从0.968%降至0.498%、减小了0.47个百分点；25个点较5个点的采样方案多设了20个样点，误差降低幅度却很小，R波段采样平均相对误差由0.563%降至0.266%、仅减小了0.337个百分点，G波段采样平均相对误差由0.507%降至0.201%、仅减小了0.306个百分点，B波段采样平均相对误差由0.498%降至0.197%、仅减小了0.301个百分点。

除了RGB波段各自的采样误差外，将采样方法在三个波段的平均相对误差取均值作为衡量采样方法准确度的另一个参考指标。5点系统采样较1点系统采样仅多4个点，三个波段平均相对误差均值由1.064%降至0.523%，减小了0.541个百分点；25点系统采样三个波段平均相对误差均值为0.208%，较5点系统采样多了20个样点，平均相对误差仅减小了0.315个百分点。

以上结果说明：样点个数为5的系统采样方案准确度较高，采样平均相对误差仅为0.523%，也可以保证较高的采样效率。

3.2　2 m面尺度到像元面尺度的采样准确度

分析像元面尺度内2 m基本采样单元的最佳布设方案时，将基本采样单元个数设为1，5，9，16，25，分别采用系统采样法和模拟退火采样法进行对比分析。

本实验组利用4块研究场地2013年—2019年的晴天高分一号2 m分辨率数据对2 m到像元面尺度的采样方法准确度进行分析，其中以2013年—2015年数据作为参考影像用于计算模拟退火采样法的样点布局，2016年—2019年的数据作为真值影像用于模拟真实地表采样、计算系统采样与模拟退火采样方案准确度，具体情况如表2所示。

3.2.1　2 m到150 m面尺度的采样准确度

计算各采样方案在150 m场地的多景真值影像上采样相对误差，将多景影像上的采样相对误差取均值与最大值，得到采样方案的平均相对误差与最大相对误差，对比采样方法准确度时以最大相对误差为主要参考指标，结果如图9所示。

图 9. 2 m到150 m尺度采样相对误差。（a）150 m新选场地；（b）150 m国家场

Fig. 9. Sampling relative error from 2 m to 150 m scale. (a) 150 m new site; (b) 150 m national site

下载图片 查看所有图片

从图9可以看出，两块150 m场地中国家场采样准确度明显优于新选场地，在进行陆地卫星验证时，建议在原国家场进行地面光谱采样。150 m国家场在采样点数相同情况下，系统采样准确度明显高于模拟退火采样，选用系统采样法较便捷且采样准确度高。当采样相对误差需控制在0.5%时，建议用5点系统采样法（最大相对误差为0.391%），与1点系统采样的0.677%相比仅多4个点就可以将采样最大相对误差降低0.286个百分点；当相对误差需控制在0.3%时，选用16个样点的系统采样法即可，采样最大相对误差可低至0.229%。

3.2.2　2 m到3 km面尺度的采样准确度

计算得到采样方案在3 km场地多景真值影像上采样相对误差，将这些值取均值和最大值作为评价采样方法准确度的指标，并以最大采样相对误差为主要评价指标，采样平均相对误差与最大相对误差结果如图10所示。

图 10. 2 m到3 km尺度采样相对误差。（a）3 km新选场地；（b）3 km国家场

Fig. 10. Sampling relative error from 2 m to 3 km scale. (a) 3 km new site; (b) 3 km national site

下载图片 查看所有图片

由3 km场地采样相对误差结果可以看出，新选3 km场地的采样稳定性与规律性明显好于3 km国家场，在进行气象卫星验证时，建议在新选3 km场地进行地面光谱采样。3 km新选场地在采样点数相同情况下，模拟退火采样准确度明显优于系统采样，建议选用模拟退火采样法。当采样相对误差需控制在1%以内时，选用5点模拟退火采样方法即可将最大采样相对误差控制在0.787%，较1点模拟退火采样法（最大采样不确定度为2.768%）多4个样点即可将最大相对误差降低1.981个百分点，样点个数为25时，模拟退火采样最大相对误差为0.238%，较5点增设20个点相对误差却仅降低0.549个百分点，此时选用5点的模拟退火采样法不仅准确度高还可以保证较高的采样效率。当采样相对误差需控制在0.5%时，选用16个点的模拟退火采样方法即可将最大采样相对误差控制在0.339%。当采样相对误差需控制在0.3%时，选用25个点的模拟退火采样方法即可将最大采样相对误差控制在0.238%。

4　讨论

针对敦煌场获取不同分辨率像元尺度地表反射率应如何采样的问题，基于高分辨率影像的采样准确度分析法对不同采样方法进行了定量比较。本实验组将采样准确度分析过程拆解为地面测量点尺度到2 m面尺度采样准确度分析与2 m面尺度到像元面尺度采样准确度分析两个步骤，通过不同尺度下的各采样方法的准确度分析结果做以下四方面讨论：

样区的均匀程度等条件对采样准确度有很大影响，在大小相同位置不同的样区采样结果有较大差异，获取高分辨率像元尺度地表反射率建议在150 m国家场、低分辨率像元尺度则建议在新选3 km场地。这是因为选择新选场地时仅采用了一个日期的参考影像，影像日期前后时段场地可能会发生变化，但150 m国家场所在区域常年被围栏保护不易遭到破坏，场地均匀性与稳定性都比较好，所以其采样准确度整体较150 m新选场地好；而3 km国家场与3 km新选场地都没有围栏保护，两块场地都比较容易发生变化，单景参考影像选择的场地有一定参考价值，所以3 km新选场地的采样准确度优于国家场。

地面测量点尺度到2 m面尺度采样，采用系统5点采样法可以保证较高的采样准确度和采样效率，采样平均相对误差为0.523%。

获取高分辨率陆地卫星像元尺度真值，建议在150 m国家场按5点系统采样方案进行2 m基本采样单元的布设测量地面反射率，此时采样准确度可以满足定标要求，最大采样相对误差为0.391%，最短采样路径仅约为238 m，工作量较小。若有特殊要求，需将2 m到150 m尺度的采样相对误差限制在0.3%以内，可以按照16点系统采样方案布设2 m基本采样单元，最大采样相对误差为0.229%。150 m场地采样方案样点分布如图9所示，5点系统采样点经纬度坐标如表3所示，16点系统采样如表4所示。

获取低分辨率气象卫星像元尺度真值，建议在3 km新选场地用5点模拟退火采样测量，该方案采样时间较短效率高且采样准确度较高满足相对误差低于1%的定标要求，能将2 m到3 km尺度转换相对误差控制在0.787%以内，样点的经纬度坐标如表5所示，遍历这些样点的最短路径长度约5.16 km，按人的慢行速度3 km/h、每个2 m采样单元的测量时长为3 min估算，约2 h可以完成全部测量工作。若需将采样相对误差控制在0.5%或0.3%以内，可以采用16点模拟退火采样法和25点模拟退火采样法，采样最大相对误差分别为0.339%和0.238%。3 km场地采样方案样点分布如图11所示，5点模拟退火采样点经纬度如表5所示，16点、25点模拟退火采样点具体位置如表6、7所示。

图 11. 2 m到3 km像元面尺度采样方案。（a）5点模拟退火采样点分布；（b）16点模拟退火采样点分布；（c）25点模拟退火采样点分布

Fig. 11. Sampling scheme from 2 m to 3 km. (a) Distribution of 5-point simulated annealing sampling points; (b) distribution of 16-point simulated annealing sampling points; (c) distribution of 25-point simulated annealing sampling points

下载图片 查看所有图片

5　结论

经上述分析讨论，考虑采样效率与采样准确度两种因素，敦煌场地获取陆地卫星像元尺度地表反射率时，在150 m国家场用5点系统采样方案标定2 m样方位置最佳；获取气象卫星像元尺度地表反射率时，在3 km场地2 m样方布设用5点模拟退火采样方案最佳；在2 m样方内建议用光谱仪进行5点系统采样测量。