基于机器学习的变换以改进外太阳系图像压缩

《Journal of Astronomical Telescopes, Instruments, and Systems》：Machine-learned transformations of outer solar system images for improved compression

　　

卡卡西号成像科学子系统（ISS）传回了迄今所见分辨率最高的土星及其卫星图像，但由于在轨压缩器过度有损，其本应产出更多的数据。事后研究发现，若压缩算法能在飞行期间自定义和更新，数据可在不损害科学目标的前提下，比实际任务中压缩 3 到 30 倍。本文提出一种机器学习框架，自动选择适配给定数据集的压缩算法。我们展示如何利用现成的压缩器加上先前开发出的神经网络压缩器，有效压缩低空间熵数据，以土卫六大气雾霾密度的卡西尼图像为例。对于空间熵更高的图像，我们用神经网络建模并推广基于小波变换的编码方式，训练其最小化图像重建误差和变换系数的可压缩性。我们展示了该压缩器在压缩到每像素 2.6 比特时（相较于无损 JPEG 的每像素 4.36 比特）能够保留卡西尼土卫八图像中用于关键点匹配的细节，并且在压缩到每像素 1.75 比特以下时，解析卡西尼土星环中微妙的螺旋桨状结构的性能优于 db5 小波变换编码。随后，我们将该网络与已知的图像光流估计网络结合，以更激进地压缩土卫八飞掠的完整图像序列。最后，我们通过对抗训练与一个区分压缩误差和噪声的判别器网络进行训练，使压缩器学会压缩具有更复杂建模噪声分布的图像。1. 引言：外太阳系数据压缩很容易忽视我们从中远距离星际探测器（尤其是外太阳系任务探测器）中提取科学信息的能力。例如，卡西尼成像科学子系统（ISS）传回了大量高分辨率的土星及其卫星图像，使得我们能够详细研究土星环和土卫六大气的异常。但实际上，我们接收到的外太阳系任务数据的每一比特，都是一项工程上的壮举——它克服了对航天器可用电力、天线尺寸以及将数据传回地球的带宽的严苛限制。鉴于这些任务的巨大成本，科学家们面临巨大压力：既要获取大量数据，又要保证数据的高质量。因此，尽管有损压缩技术有望大幅提高数据回传量，许多科学团队却对它心存戒备。这种犹豫很大程度上源于事实：标准的"即买即用"压缩器会给数据引入伪影，某些伪影甚至会高于其压缩数据所用的仪器本身的噪声水平。例如，卡西尼 ISS 配备了一个硬件 JPEG 压缩器，研究发现它在所有可用配置下都会在相机图像中引入严重的块状伪影。这一缺陷直到它已经在轨运行后才被发现，科学与工程团队最终都决定不使用它，转而选择下载更少的未经压缩的原始数据。接续文献 1 中的讨论，我们主张采用一套可快速在飞行期间更新的、软件实现的、可调的压缩算法套件。每个科学团队都可以选择合适的压缩模型或自行构建，并选择团队认为能充分保留数据科学质量的相应参数。该协议最适合长期任务，可以让科学家先下载一批数据进行实验以表征噪声、开展反演并相应地微调压缩算法。压缩器随后可上传至航天器进行长期使用，并可根据科学团队判断周期性重复此过程，以适应对任务和仪器性能新获得的认识。需要强调的是，我们主张压缩器应当具有相对较低的复杂度，既是为了最小化运行成本，也是为了便于高效上传至远距离航天器。文献 1 中已证明，现成的压缩器可被微调以有效适配特定数据集，而不损害其科学完整性，我们将在第 2.1 节中回顾。这项工作由一组科学家和压缩专家团队协作，相互指导，分析数据并对各种压缩算法进行系统测试。本文中我们更进一步，将机器学习引入以自动化压缩器及其参数的选择。具体来说，我们用神经网络推广前一项研究中表现良好的压缩模型，并把科学反演的要素纳入训练过程和神经网络结构。我们设想，可能代价高昂的模型训练将在地面完成，使用早期回传数据，或者更早地基于早期任务回传数据进行数据模拟。只有最终模型参数（以及相关算法和流程）会被上传到航天器。我们的新贡献如下：• 我们将机器学习驱动的有损压缩直接关联到科学数据中底层关注信号的提取问题；• 我们展示了在某些情况下机器学习可以补充先前的算法，并在卡西尼图像提取土卫六分离雾霾（文献 2,3）上验证了性能；• 我们提出了实现机器学习变换编码的神经网络模型与算法，这是对文献 1 中算法的推广。我们的模型类似于变分自编码器，并以文献 4–8 类似的方式，使用率失真驱动的损失函数进行训练；• 我们在卡西尼任务的两个数据反演示例上展示了它的性能：1. 土卫八飞掠序列图像中的关键点匹配，以及 2. 测量土星环中的螺旋桨状结构；• 我们还讨论了对原算法的修改，以纳入机器学习估计连续图像之间的光流，以及使用与我们的压缩/解压网络对抗训练的判别器网络来区分压缩误差与噪声。我们强调，本文的工作并不旨在对所有现有压缩方法（无论是使用还是不使用机器学习的众多方法）进行全面比较。我们也不声称本文所述算法能胜过市面最好的压缩器。相反，我们关注的是科学家们在机器学习协助下评估有损数据压缩器对其反演的影响并加以修改的能力。我们希望展示一个框架的可行性，该框架可以向科学家提供一套压缩工具，并允许他们利用机器学习快速微调其中一种，以满足其特定需求。具体而言，我们只考虑那些在文献 1 中被科学家标记为对其反演效果较好的几种压缩方法的机器学习改进。我们花了大量篇幅讨论科学反演本身的结果，以及它在机器学习增强前后受到有损压缩使用的影响。此外，我们专注于对数据进行机器学习的变换，将其转化为更可压缩的格式，同时保持其科学完整性。这包括使用科学信息驱动的失真度量来指导训练过程，以及利用对反演方法和底层科学的了解来设计模型本身。我们不涉及对压缩过程后续步骤的重大机器学习修改，例如变换系数的量化和熵编码。我们的系统是数据驱动的，在更具针对性的数据集上训练时能获得更高的可压缩性，并允许在"纳入关于数据环境结构和噪声的先验知识"与"让网络自行学习这些信息"之间进行权衡。此外，我们的算法具有可移植性，我们从轻量级网络中获得了高性能，这些网络可以快速训练并轻松上传。在我们的一种方法（第 3 节）中，网络甚至可以作为一种在地面的后处理形式来应用（与文献 9 类似，后者利用机器学习修复有损压缩伪影）。这可以进一步减轻对航天器计算资源的需求。本文的部分结果已在文献 10 中发表。1.1. 论文组织结构本文组织结构如下：第 2 节回顾了卡西尼 ISS 图像压缩以及机器学习应用于数据压缩的先前工作。第 3 节举例说明现成压缩器如何借助神经网络增强，以有效压缩平滑图像，并以土卫六大气温霾密度的反演作为演示。第 4 节使用神经网络编码器/解码器模型推广基于小波的变换编码，产生低熵的变换数据分布，并在以下两方面展示其有效性：1. 土卫八图像中的关键点匹配，以及 2. 土星环中螺旋桨状结构的反演。第 5 节将光流估计纳入我们的框架以压缩土卫八的飞掠图像序列。最后，第 6 节将我们的压缩器与区分噪声和压缩伪影的网络进行对抗训练。这使我们既能针对复杂的噪声分布调整压缩器，也能标记可能的压缩伪影。2. 相关工作2.1. 先前研究：卡西尼 ISS 的变换编码卡西尼 ISS 搭载了广角和窄角相机各一台，每台配备 1024×1024 像素 CCD 探测器阵列，2004 至 2017 年在轨环绕土星期间共传回数十万张图像。数据以每像素 12 比特存储，可选择使用平方根查找表转换为 8 bpp 格式，从而对像素值的动态范围进行更粗略的量化。卡西尼图像相应地以数字量（DN）表示，12 比特图像的取值范围为 [0, 4095]，8 比特图像为 [0, 255]。ISS 搭载了一个硬件实现的 JPEG 有损压缩器，它对 8×8 像素块进行操作，执行离散余弦变换、量化和 Huffman 编码。但发现它在图像中产生不可接受的瑕疵，在任务后半段未被使用。它还配有一个无损 Huffman 压缩器，但受限于在高熵区域截断交替行的机制，它无法将高熵场景压缩超过两倍。一项事后研究（文献 1）探索了对一小批土星、其环系及其卫星土卫六与土卫八的未压缩图像使用更复杂的变换编码形式，以实现更高的压缩比。变换编码顾名思义，即对每幅图像施加一种变换（通常是线性的），将其转化为一组基元素的系数，希望这些系数的熵低于原始图像的像素。JPEG 压缩即为例，采用离散余弦变换。研究发现，几种常见的小波（文献 11）和小波原子（文献 12）变换都能对图像产生较为稀疏的表示。这两种变换都可以通过多分辨率滤波器组实现：首先对图像施加低通和高通滤波器，然后对前一次迭代中低通滤波器的输出迭代重复施加。作者通过清零可忽略的系数并施加熵编码，实现约 10 倍的图像压缩，同时保持较低的平均平方重建误差。一般而言，为给定数据集和科学需求寻找合适的变换极具挑战。这是因为很难找到一组基，它不仅能紧凑地表示数据，而且不会因压缩而在数据中引入系统性误差，从而可能损害科学结论。例如，JPEG 压缩器在卡西尼图像中产生的块状伪影就对应于离散余弦变换所用滤波器的形状，它源于压缩丢弃了某些系数。2.2. 利用机器学习增强压缩将机器学习应用于数据压缩的研究有着丰富的历史，特别是涉及神经网络的使用。

我们引导读者参考文献13以获取广泛的概述，其中指出基于神经网络的图像压缩早在80年代末就已受到关注。14,15 那个时期也见证了自编码器网络的首次出现，16,17 自编码器长期以来一直是一种流行的方法，旨在利用神经网络将数据映射到更低维的空间，随后再将其恢复。随着机器学习和深度学习的革命蓬勃发展，我们拥有了大量的新工具来应对数据压缩问题，尤其是针对图像和图像序列的压缩。一些研究使用循环神经网络来压缩和重构图像及其误差残差，18,19 而其他研究则设计了基于变分自编码器（VAEs）的压缩器并结合了侧信息，4,6,20,21 同时训练一个超先验模型来捕获数据统计信息以辅助压缩。5 我们将提出的一些模型类似于β-VAEs22，具有高斯潜随机变量，类似于文献23和24中讨论的模型。其他研究采用了自回归模型25–27来估计图像概率密度，或使用生成对抗模型（GANs）28–30来尝试隐式建模图像分布以生成图像。图像序列压缩与视频压缩密切相关，后者在文献31–33等工作中已通过机器学习方法加以探索。这些研究展示了机器学习可应用于视频压缩过程中的各个环节，包括追踪重叠图像间场景中物体的物理位移和运动，以及压缩光流（即连续帧之间像素的映射关系）。机器学习驱动的光流估计和压缩本身就是一个有趣的问题，文献34–38中提出了一系列基于神经网络的相关模型。我们在第4节中讨论的模型是变换编码的神经网络推广形式。这涉及训练模型学习如何将图像变换到一个可以通过量化和熵编码其变换系数来实现高效压缩的域中。在我们的情况下，我们将讨论离散小波变换（DWT）的神经网络推广形式。类似的方法（使用机器学习来推广或增强离散小波变换或其他类型的变换，例如离散余弦变换）已在文献39–44等工作中得到研究。对于我们提出的深空压缩框架而言，这是一个有利的选择，因为我们发现它使我们能够使用相对较小且可移植的模型实现良好的有损压缩性能，这些模型可以轻松地在任务期间上传到航天器上。值得注意的是，近期的一些研究方向试图利用基于Transformer的架构45–47，这些架构通常包含非常大量的参数。这些模型有潜力超越更简单的架构，但训练难度大且需要大量的内存和计算资源，因此我们不在本工作中考虑它们。

我们主要关注这些工具在多大程度上能够方便地融入外太阳系图像压缩任务中，例如卡西尼号任务的图像。为此，我们首先研究如何编码文献1.3中所述有损压缩算法产生的图像误差残差。使用现成模型增强土卫六雾霾数据的有损压缩

作为第一个测试，我们考虑压缩土卫六图像的问题，其大气雾霾在其边缘周围呈现微妙的特征。为了突破文献1中描述的小波压缩算法的极限——该算法用每个图像占主导地位的db5小波变换系数的一部分来近似图像——我们尝试仅用3%的小波变换系数来重构图像。我们选择了那些在我们的训练集图像上具有最高方差的系数。这一选择基于几个原因：首先，我们预计图像中方差最低的部分对于机器学习模型来说最容易预测。其次，如果我们的模型仅能学习缺失系数的平均值，我们的选择至少可以保证重构中的总体均方误差最小化。

如果我们将以64位精度编码所保留的3%小波系数，这对图像而言相当于总共1.92 bpp。在将所有选定高方差系数以外的系数置零后，我们执行了逆小波变换，并计算了相对于原始图像的误差残差。正如预期的那样，解压缩过程在土卫六边缘周围产生了波纹状伪影，对应于小波滤波器的形状。

然后，我们将文献19中的神经网络编码器应用于误差残差的编码。该网络如图1所示，对32×32像素的图像块进行操作，通过循环卷积层将每个图像块编码为2×2×32的二进制编码。我们对每个图像块应用该网络八次迭代，从而对误差残差实现总体1 bpp的编码，并使用土卫六边缘附近的图像块进行网络训练，该区域包含了关于雾霾的大部分信息。我们将解压后的残差与原始小波编码图像相加，并检查了由此产生的新误差残差。令人印象深刻的是，我们发现该网络不仅恢复了图像的边缘细节，而且仅通过保留1%的原始db5小波系数就实现了这一点，总压缩率仅为1.64 bpp。在图2中，我们展示了土卫六雾霾图像、从3%的小波系数恢复时的压缩伪影，以及使用1%的小波系数和误差残差压缩网络纠正后的图像。

图1：下载全尺寸图像——文献19中用于在小波压缩后对图像残差进行编/解码的神经网络单次迭代。该网络应用一系列卷积和循环卷积层来压缩和解压32×32的图像子区域。每次迭代压缩相对于前一次的误差。默认情况下，我们对该网络应用八次迭代，总共添加1 bit/像素。

图2：下载全尺寸图像——左：原始图像（土卫六雾霾），动态范围[0, 4095] DN。中：从存储3%的db5小波变换系数重构图像时的误差（压缩至1.92 bpp）。右：仅使用1%的db5系数并借助神经网络时的误差（压缩至1.64 bpp）。这消除了大部分压缩伪影，同时实现了更高的压缩率。两幅误差残差图像均以动态范围裁剪至[?10, 10] DN的区间显示。

该网络甚至能够恢复土卫六雾霾最微妙的特征之一——所谓的"分离式雾霾突起"。这是雾霾密度在边缘正上方出现的明显凹陷和隆起，可以通过绘制一条穿过土卫六边缘的图像线来描绘。在图3中，我们通过绘制使用多达10%的小波系数重构图像时的重建线，以及保留1%系数并借助神经网络的情况来说明这一点。我们看到只有后者能够恢复分离式雾霾突起并与原始图像匹配。

图3：下载全尺寸图像——土卫六雾霾原始图像与重构图像的像素值。在存储1%的db5小波系数后对图像残差应用该网络，我们恢复了红框中标注的特征性"分离式雾霾突起"。在没有神经网络辅助的情况下，即使存储十倍数量的系数也无法恢复分离式雾霾。

在土卫六雾霾测试之后，我们尝试将此方法应用于冰卫星土卫八（Tethys）的图像，这些图像具有比土卫六图像更高的空间熵。我们用其db5系数的1%来近似原始图像，并再次对32×32像素的残差图像块应用该网络。结果如图4所示。在此压缩级别下，我们可以看到虽然网络对重构图像做出了明显的改进，但它无法捕捉到我们在土卫八图像中期望的更精细的细节。因此，虽然我们的土卫六雾霾测试表明，对于相对平滑的景观，有能力利用现成模型有效补充有损小波变换编码以提高整体压缩率，但我们仍需额外的工具来捕捉更复杂的图像。为此，我们接下来试图开发一种适用于任意图像类别的机器学习框架，并有望将失真降低到图像固有噪声水平以内。

图4：下载全尺寸图像——（左）土卫八原始图像的一部分。（中）用其db5小波系数的1%主导系数近似图像后的同一部分。（右）对每个32×32像素图像块应用残差神经网络进行纠正，这改善了质量但仍无法捕捉到图像的完整分辨率。

4. 机器学习设置：损失函数与模型结构

如第3节所述，图1所示等现成模型可以改善某些数据集（如平滑的土卫六雾霾图像）在压缩率和图像科学内容失真之间的权衡关系。然而，同样的模型在其他数据集（如复杂的土卫八图像）中可能会无法捕捉到感兴趣的内容而无需一些修改。我们的目标是拥有适用于广泛数据类别的简单工具，以及以自动化方式将其微调至特定数据集的框架。

在典型的机器学习压缩框架中，给定相似分布的随机数据向量x∈?^N×1，我们寻求为编码器?和解码器??拟合模型，分别从各自的模型候选集??和???中选择，以最小化函数f?,??(x)（代表x的重构误差）和g?,??(x)（代表压缩成本）的线性组合。然后我们寻求执行优化方程（1）：

最小化 ?∈??, ??∈??? ??(f?,??(x)) + λ·??(g?,??(x))

其中λ是精心选择的常数，??(·)是我们对x分布上的期望算子。通常，我们的数据统计信息是未知的，我们依赖于基于随机梯度下降的算法，使用大量训练数据来执行方程（1）中的优化。

我们的目标是将我们对久经考验的特定数据压缩方法的理解作为先验知识加以利用。我们不希望从头开始使用神经网络构建压缩器，而是希望优化此前研究的"现成"压缩算法，重点在于产生相对较低复杂度和可移植的压缩器，以便上传到航天器进行低成本使用。为此，我们概述以下高层流程：

1. 构建编码和解码神经网络?和??，其网络拓扑推广了此前研究的编码和解码变换。
• 在编码器网络中包含额外的可训练变量（由向量ξ表示），代表变换系数?(x)的统计分布。

2. 在特定数据集上训练网络，最小化损失函数方程（2）：

C(x) ? ‖x ? ??°?(x)‖2 + λ·L(ξ|?(x))

其中L(ξ|?(x))是负对数似然，基于由ξ参数化的y ? ?(x)的建模分布P(·)计算：L(ξ|?(x)) ? ?(ne) ∑????? log P(y?|ξ)。

这里，y?是编码向量y的第i个系数，ne是编码系数的数量。在我们建模的分布中，我们出于简化目的隐式假设编码系数y?之间相互独立。

3. 对变换系数?(x)进行量化，并使用熵编码对量化后的系数进行编码。

通过这种成本函数的选择，方程（1）中的目标函数变成了受失真约束的码率优化的拉格朗日函数，类似于文献4–8、48中分析的内容。让我们讨论策略中的几个选择。通过显式建模由ξ参数化的潜变量分布，我们可以在训练过程中提示潜在的信号统计特性。这使我们能够融入关于数据方面的知识，甚至尝试在训练过程中检索信号统计特性本身。注意，虽然通常考虑均方误差函数作为f?,??(x)，但我们选择L2范数误差，因为其在误差变小时梯度不会消失，这使得它在低噪声环境下非常理想。我们选择负对数似然函数来表示压缩成本项g?,??(x)，其动机在于最小化其在真实数据分布Ptrue上的期望等价于最小化Kullback-Leibler（KL）散度DKL(Ptrue||Pξ)，其中Pξ(y) ? ∏????? P(y?|ξ)是我们的建模分布（这里，ne是编码系数的数量——即y的维度）。

KL散度表示使用为我们建模分布优化的编码，对从真实分布中抽取的系数 {????}??????=1 进行编码所需的额外比特数。当 ???? 与 ??true 匹配时，它会在 ?? 上达到最小值。

4.1. 神经网络设计

在此阶段，我们借鉴了文献1中关于压缩卡西尼号ISS图像已知方法的有效经验。研究发现，即使仅保留小比例的小波变换系数，通过Daubechies双正交小波实现的小波变换也能很好地描述图像。该变换通过对图像的每个维度施加5抽头低通和3抽头高通滤波器实现，产生四个子带系数（图5），然后迭代地对低通子带应用滤波器（图6）。

图5
第一级小波变换的泰坦雾霾图像。（左）原始图像，（中）对原始图像的每行应用低通和高通滤波器，产生新图像。每行现在包含低通系数和高通系数的连接。 （右）对每列应用低通和高通滤波器，产生四个不同的子带系数。

图6
完全变换的泰坦雾霾图像，通过对上一级最低子带系数迭代应用小波变换实现。

为了在编码器神经网络模型 ? 中捕捉这种变换的精髓，我们用四个卷积核替换小波变换滤波器；每个的输出代表四个子带之一（图7）。我们选择每个卷积核大小为8×8，以大致匹配Daubechies db5小波滤波器的大小，并在每一层引入绑定偏差和ReLU激活。由于我们允许在每一层使用不同的卷积核，我们期望该模型能包含比典型小波变换（在每级应用相同滤波器）更丰富的编码器集合。重要的是，我们不假设在每次卷积级之后进行降维。换言之，四个卷积层每个输出的系数数量等于输入尺寸的四分之一，?(??) 中最终编码系数数量等于原始图像??的像素数量（尽管在实践中，编码器会将许多这些系数驱动至零）。我们通过每个方向长度为2的卷积步长实现这一点。我们对每个卷积层的输入应用零填充，在输入的每一侧添加 ??=3 个零。为构建解码神经网络 ??，我们模仿编码网络拓扑，按反序排列层。对于解码网络，我们用转置卷积层替换卷积层，每个独立训练。

图7
神经网络编码器泛化的??级小波变换（由于空间限制，此处展示3级，尽管实践中通常选择??=5）。(a) 一个32×32像素图像通过四个8×8卷积层。 (b) 在每级，第一层的输出迭代通过另外四个卷积层。 (c) 最终级的第一层输出与所有级的第二、第三和第四层输出组合，向量化为1024个系数。

除了卷积核权重外，我们还包括一组可训练参数，用于编码器网络输出系数的概率分布。这使我们能够根据建模分布计算编码系数微分熵的代价。为简化起见，我们尝试为每个系数拟合独立的正态分布（见图8），需要为每个输出系数 ???? 训练均值 ???? 和标准差 ????。根据参数化为 ??={????,????} 的假设分布，我们编码系数 ??=[????] 的相关可压缩性代价为公式(3)：

??(??|??) ? ?∑? log ??(???|???,???) = ?∑? log(1/√(2?????) exp(?(???????)2/(2???2)))

然后我们训练 ?、?? 的权重和参数 ??，以最小化来自公式(2)的总体损失函数 ??(??)，同时受建模高斯分布的最小标准差约束，???≥??min，这有效地约束我们的编码器以防止损失函数发散。

图8
对于编码器输出的每个系数，训练一组概率分布参数以最好地表示数据。在此示例中，训练均值???和标准差???以最小化系数???的微分熵，假设为正态分布。

我们注意到，由于 ??(??|?(??)) 与上文讨论的 ??_KL(??true||????) 相关，训练过程将驱使参数 ?? 形成与变换系数的经验分布匹配的分布 ??(???|???,???)（同时驱使编码变换匹配参数化分布），本质上是在变换域学习数据统计。这种使用神经网络学习数据统计信息和变换的过程，与变分自编码器20（将输入数据向量映射到采样获得编码向量的概率分布）和基于流的生成模型21,49–51（使用正态化流方法寻找从简单随机变量到数据向量的变量替换，以学习数据集的潜在概率分布）技术类似。

备注1. 公式(3)中 ??(??|??) 的选择特别适用于 ?? 具有多元高斯分布 ??(????,????) 的情况。在这种情况下，如果我们将模型候选 ?? 和 ??? 限制为与 ?? 兼容维度的一致矩阵，则最小化 ??(??(??)) 将驱使训练参数趋向真实统计分布的各种特征。特别是，对协方差矩阵进行奇异值分解，????=??????????????，编码器矩阵和解码器矩阵将分别收敛于??????和????。类似地，???2将趋向????的奇异值，而???将满足??([???]??????=1)=????，其中??是训练的解码器。为详尽起见，我们在附录中具体展示。注意在这种情况下，变换系数???将成为独立的高斯随机变量，其各自分布确实是??(???,???2)，正如建模所示。这使我们能够独立地量化和熵编码它们，以压缩至给定均方误差失真下的最优比特率（如文献52中的反向注水论证所示）。

由于我们的意图是使网络小巧且可在有限数据上快速训练，我们选择32×32像素的图像块作为编码网络的输入图像。从1024×1024像素的卡西尼图像样本中，我们将每个图像分割为1024个非重叠的32×32子块。我们在编码器中包含五级卷积层（注意由于空间限制，图7展示了较小的三级编码器）。根据应用情况，我们修剪图像块以捕获可能包含感兴趣特征的块来构建训练集和测试集。此外，我们对所有图像块进行单次平移和重新缩放归一化，使其值落在0到1之间。

由于我们使用绑定偏差，每个卷积核总共贡献 8×8+1=65 个可训练参数。由于我们有五级小波级，每级有四个卷积核，编码器或解码器网络中的参数总数为 5×4×65=1300。此外还有 2×32×32=2048 个参数用于建模与编码器输出系数相关的分布的均值和标准差。因此，此模型总共有 (2×1300)+2048=4648 个可训练系数。

4.2. 选择系数量化级别

一旦训练了编码器和解码器网络，我们就有了将图像块??变换为与原始块像素数量??相同维度的向量空间的方法，但?(??)的每个条目都来自低熵分布。我们的实际压缩源于对????(??)的条目进行量化。由于我们为每个条目???拟合了各自的独立分布??(???,???2)，我们选择值δ?并执行量化映射公式(4)：

??(δ?)?: ??? ? ??? ? ?(???????)/δ??

这将变换系数映射到距离估计均值???的δ?倍数的整数倍。为恢复???的近似值，我们取相关区间的中点：公式(5)

?? ? ??? + δ?·(??(δ?)?(???) + 1/2)

我们在图9中展示。

图9
量化变换图像块的示例，蓝色表示原始系数???，橙色表示量化系数???(???)。

为简洁起见，我们写为公式(6)

?? ? [δ?,…,δ??]?

公式(7)

??_??(??) = [??(δ?)?(???) ? ??(δ??)??(????)]

我们期望每个量化系数接近零，并对它们的值应用算术熵编码器以获得最终编码。一旦我们使用编码器网络?变换图像的1024个子块中的每一个，我们将编码向量 {???}?=1,…,1024 连接成矩阵公式(8)

?(??) = [??? … ??????]

其中第??行包含每个图像子块的第??个编码系数???。定义??????_??(???)，相应的量化矩阵为公式(9)

??_??(?(??)) ? [??? … ??????]

然后我们可以选择将整个向量化矩阵传递给熵编码器，或选择分别编码每一行。后者通常是我们的首选方法，不仅因为每行应对应于相似分布的量化系数（由于我们的优化），而且因为它允许我们独立选择量化级别δ?。因此，我们定义相关熵公式(10)

??(??_??(?(??))) ? ∑?=1? ??(???)

其中???是公式(9)矩阵的第??行，??(???)表示???离散量化值分布的熵。对于用于优化向量??的图像块集合 ???{???,…,???}，我们可以考虑平均熵公式(11)

??_??(??) ? (1/??)∑?=1? ??(??_??(?(???)))

这是编码图像所需预期比特数的代理，随后公式(12)

β_??(??) ? (1/??)??_??(??)

即每像素的预期编码比特数。为选择量化级别，我们选择期望的每像素比特数??，并将??初始化为均匀细量化级别，????[δmin,…,δmin]?。然后我们对向量??执行梯度下降，以最小化每像素比特数：对于足够小的固定学习率dδ，我们计算公式(13)

??? = ????? ? dδ·?β_??(?????)

迭代直到差值 |β_??(??) ? ??| 在某个选定阈值内。在实践中，我们发现这种最速下降方法能够充分最小化目标压缩级别的量化误差。

4.3. 在土卫卫星忒瑞斯忒斯图像上的锚点匹配

为测试我们神经网络对基于小波的压缩方法的泛化能力，我们在土星冰卫星忒瑞斯忒斯的图像样本上训练了模型，其表面代表高空间熵场景。图10展示了一个这样的例子。作为基线，我们展示从10%的db5小波变换系数中以32位精度存储的此图像32×32像素块的重建——压缩到约3比特/像素——揭示了大量条纹伪影和误差残差中的大量信号结构（图11）。

图10
卡西尼ISS拍摄的忒瑞斯忒斯图像(N1477639092_1)。颜色对应于数字（DN）中的像素值。

图11
通过用db5小波变换主要10%的系数近似（压缩到约3 bpp）压缩和恢复忒瑞斯忒斯图像子块。我们看到明显的压缩伪影，并丢失了原始图像地形的某些特征。颜色对应于缩放的像素值（DN/256），使所有值落在0到1之间。

接下来，我们在来自另一组忒瑞斯忒斯图像的块集上训练了我们的神经网络编码器/解码器，并将系数量化到0.5到7 bpp之间的几个压缩级别。我们在图12中展示了几个压缩和重建图像及其相关的误差残差。我们看到图像重建的显著改善，即使压缩到1.36 bpp时，图像的重要特征仍可识别。

图12
使用神经网络泛化的小波变换编码器压缩和恢复忒瑞斯忒斯图像子块。颜色对应于缩放的像素值（DN/256）。第一行描绘同一图像压缩至1.36 bpp（左）、2.12 bpp（中）和2.95 bpp（右）。第二行描绘与原始图像相关的残差误差。

经网络优化的编码器比先前的db5小波变换更好地保留了图像细节，即使在更高的压缩比下，误差的动态范围也小了2到4倍。为了验证我们的压缩器，我们考察了它对Cassini探测器飞掠Tethys卫星序列图像进行相邻图像间控制点匹配能力的影响。该过程包括识别图像之间的相应地形特征并匹配相应的像素，以便将一幅图像扭曲变换到下一幅图像。我们选择这颗土星卫星是因为其明亮的冰质表面具有相对均匀的反照率，因此大部分可见亮度变化可归因于地形起伏而非区域成分异质性。我们重点关注两张Cassini窄角相机图像N1477638891和N1477638933（表1），因为它们非常相似，但观测和照明几何条件并不完全相同，而且这些图像在从探测器下载到地球时是无损压缩的。Tethys的平均直径为1062公里，在每幅图像中对应约694个像素的直径。每幅图像中的49.7°相位角将Tethys表面上的太阳直射点——即太阳位于正上方的位置——置于距圆盘中心向亮边缘方向约50°经度的位置。

表1 Cassini拍摄的Tethys立体控制点测试图像详情

| Cassini图像 | 航天器距离（公里） | 子航天器点（纬度，西经度） | 太阳直射点（纬度，西经度） | 相位角 | 最高分辨率（公里/像素） |
|---|---|---|---|---|---|
| N1477638891 | 254,923 | 22.16°，270.47° | ?22.87°，291.1° | 49.7° | 1.5237 |
| N1477638933 | 254,887 | 22.15°，270.50° | ?22.87°，291.2° | 49.7° | 1.5235 |

对于我们的控制数据库，我们首先使用JPL的VICAR" PICREG"监督式控制点匹配程序手动收集与图像N1477638891和N1477638933中均可见的清晰撞击坑相关联的控制点。这些无损压缩的全圆盘图像是在相位角条件下获取的，展示了表面对比地形特征的外观如何随太阳入射角和发射角的大范围变化而变化。围绕每个控制点的64×64像素框的匹配一致在所有位置上都优于0.7像素。然后我们使用自己的自动模式匹配程序（"Feature-Track"），尝试将参考图像（N1477638891）中的每个像素与测试图像（N1477638933）进行关联。

Feature-Track工具包?3??是一组FORTRAN计算机程序的集合，允许用户半自动地逐像素测量立体对字节格式数字图像的视差偏移量。这些程序产生用于立体摄影测量的位移场数字地图，以及相应的逐像素相关统计量，用以表征参考图像中在配对图像中被扫描的相关窗口内特征匹配的质量。

位移场提供了从立体对图像生成数字地形图的手段。在搜索开始时，图像被缩小到起始尺寸级别。随着搜索的进行，图像在每一级上以两倍因素逐步放大，直到在全尺寸级别（我们的测试中为1024×1024像素）找到最佳结果。在每个阶段，程序逐像素遍历参考图像，提取可变大小的相关窗口，然后在配对图像中的相应试验位置上逐步扫描该窗口，测量每个点窗口亮度模式的平均相关值，寻找最大窗口相关值。相关窗口的尺寸在用户定义的范围内可变（表2），Feature-Track具有根据窗口内的亮度分布自我调整窗口尺寸的能力。用户还可设置相关阈值，使Feature-Track能够在即使最佳找到的相关值低于下限的情况下，根据周围高相关基础上的偏移量插值控制点（样本、行）像素偏移量。

表2 Feature-Track运行的四阶段控制点匹配方案

| 阶段级别 | 精度（像素） | 图像尺寸缩减 | 相关窗口尺寸（最小，最大） | 最大扫描范围（Δ样本×Δ行） | 插值相关阈值（起始，结束） |
|---|---|---|---|---|---|
| 1 | 1/8 | 256×256 | 1024×1024 | 7, 31 | 12×12 | 0.10:0.85 |
| 2 | 1/8 | 1024×1024 | 1024×1024 | 7, 31 | 5×3 | 0.10:0.85 |
| 3 | 1 | 1024×1024 | 1024×1024 | 9, 17 | 7×7 | 0.10:0.85 |
| 4 | 1/8 | 1024×1024 | 1024×1024 | 7, 19 | 0.5×0.5 | 0.10:0.90 |

在其原始形式中，Feature-Track被设计用来匹配像素位移场，以关联配对图像，其中配对中的两幅图像都被数字重映射到平面投影中，所有区域以相同的水平比例显示，并裁剪为仅显示相互重叠的区域（参见参考文献55）。然而，在我们的研究中，我们需要以探测器下发的原始未修改形式表示行星图像。因此，我们的Cassini拍摄的Tethys图像没有任何重映射；它们的尺寸比例略有不同，Tethys在每幅图像中的位置也略有不同，并且图像包含卫星体形状和曲率造成的透视缩短和其他畸变效应。此外，我们没有进行任何常规的有系统性的恢复处理来消除仪器噪声和随机噪声效应。然而，由于图像包含大面积的暗天空，我们排除了距Tethys圆形边缘超过其半径5%的图像部分。

此前使用Feature-Track的经验建议采用四阶段控制点搜索策略（表2）。阶段1代表粗搜索，其中全尺寸（1024×1024像素）图像首先缩小到256×256像素，然后逐步以两倍放大，从参考图像到其配对图像逐步微调像素（样本、行）偏移量的阵列。阶段2输入阶段1的结果，以全尺寸精度搜索至1/8像素精度。在阶段3中，输入阶段2的结果，以更紧的、更小的相关窗口尺寸范围执行1像素精度的搜索。最后，在阶段4中，输入阶段3的结果，搜索重复至1/8像素精度，但允许相关窗口尺寸略大于阶段3。

一般而言，对于立体摄影测量，期望像素偏移不确定度达到亚像素精度。我们在无损压缩图像对之间的基线自动模式匹配测试中，相对于我们的控制数据库，得到平均控制点误差为0.17（标准差0.61像素）。对于任何给定的控制点，我们考虑相对于控制基准偏移量1像素或更大的误差为不可接受的。对于100个控制点的集合，我们定义当平均相关值加一个标准差超过1像素宽度时，即为压缩水平将图像保真度降低到超过可用范围的上限。在Feature-Track原始设计的应用中，我们通常期望平均控制点误差和相应的标准差会随压缩比的增加（即每像素位数的减少）而均匀增加。

然而，在我们使用Tethys压缩/解压图像对组进行测试时，我们发现可见圆盘上相关值存在显著的系统性变化，导致平均控制点误差及其标准差随压缩状态增加呈现出更复杂的趋势。其中一些趋势根据图像本身的基本特性是可以定性预测的。例如，由于图像包含像素级的随机噪声效应，我们期望少量的压缩将表现得像简单的低通滤波器，衰减像素级随机噪声，可以将图像之间的相关性改善到优于原始基线测量值。

同时，我们预先认识到全圆盘图像中的地形细节在垂直观测几何条件下（靠近圆盘中心）的分辨率将明显更高、失真更小，而靠近边缘或晨昏线处则较差。这在太阳直射点附近（也靠近照亮边缘）尤其如此，因为高太阳角度和Tethys的高表面反照率相结合，使地形阴影和地形坡度效应减弱，可能导致可匹配的小尺度表面细节很少。因此，解压图像在选定的近垂向区域表现最理想，而靠近卫星边缘时表现较差，因为在边缘处地表特征以更大的斜射角观测，地形定义在参考图像和测试图像之间被削弱且空间变化较大。Feature-Track可以通过将控制点的相关值赋予空值来拒绝差匹配或未匹配的控制点。这种方法避免了匹配那些过于平淡无特征的相关窗口——它们否则会被登记为100%相关。在我们对Feature-Track的改编中，一个意想不到的结果是，匹配特征的系统性模式识别出每个场景中适合立体摄影测量的部分和应当裁剪掉的部分?。

4.4. 结果

我们调查了两种情况：案例#1中参考图像和配对控制点图像在相同压缩水平下进行有损压缩，案例#2中参考图像无损压缩而仅配对控制点图像有损压缩。这两种情况遵循了Cassini成像实验在飞掠序列中设定的模式，其中一系列有损压缩图像周期性地穿插着一幅无损压缩图像。

基线控制点：我们的基线控制点搜索是通过将Feature-Track应用于原始未压缩的Cassini图像获得的（表2）。参考图像N1477638891及其100个人工选择的控制点位置如图13(a)所示。图13(b)显示，即使使用无损压缩图像，仍有31个控制点被拒绝，原因是极低的correlation值导致未匹配、相关窗口在图像边缘被截断，或匹配质量差（即相关加权像素偏移误差大于1.0像素）。

图13 Tethys两幅未压缩飞掠图像的基线控制点。(a) 参考图像N1477638891上的100个人工选择的控制点，使用PICREG手动与测试图像N1477638933上的对应点进行匹配。(b) 参考图像和测试图像之间强匹配、弱匹配（或未匹配）的控制点。(c) 图像之间对应控制点周围窗口的相关值直方图。(d) 对未压缩图像应用Feature-Track时其自动控制点预测的相关加权误差直方图。

正如预期的那样，几乎所有匹配良好的值都围绕卫星圆盘中心的垂线点分布。大部分未匹配点出现在照亮边缘附近和太阳直射点附近。控制点13、39、43和44均位于两幅图像的最顶部边缘附近，那里Tethys的圆盘被突然截断。最后，控制点14和17恰好位于晨昏线上，那里场景中大部分由低DN值像素组成，背景噪声比明亮区域的噪声更显著，并且掠射照明角度的微小变化可能导致场景的大幅变化。

在图13(d)中，我们展示了相关加权误差的直方图（两幅图像之间手动确定的控制点偏移量与使用Feature-Track自动计算的偏移量之间的不匹配）。直方图显示了一个强单峰分布，平均误差为ε=0.13像素，标准差为σε=0.24像素。这些值将作为我们在以下测试结果中用于比较的基线（见表3和表4，它们显示了经训练网络所实现的平均压缩位数/像素下的性能）。我们在图14中绘制了加权平均误差，并为案例1和案例2的数据点分别拟合了独立的经验趋势曲线。对于相应的标准差图，一条经验趋势曲线被拟合到案例1和案例2的合并数据上。

表3 小波网络：两幅图像均有损压缩（案例1）

| 压缩（位数/像素） | 平均控制点误差（像素） | 标准差（像素） | 平均+1标准差（像素） |
|---|---|---|---|
| 无损 | 8.00? | 0.13 | 0.24 | 0.37 |
| 有损 | 6.43 | 0.18 | 0.31 | 0.49 |
| 有损 | 5.04 | 0.22 | 0.51 | 0.73 |
| 有损 | 3.25 | 0.16 | 0.31 | 0.47 |
| 有损 | 2.09 | 0.16 | 0.18 | 0.34 |
| 有损 | 1.58 | 0.25 | 0.33 | 0.58 |
| 有损 | 1.21 | 4.36 | 2.73 | 7.09 |
| 有损 | 1.07 | 1.57 | 2.01 | 3.58 |

?无损未压缩的位数/像素

表4 小波网络：参考图像无损压缩（案例2）

| 压缩（位数/像素） | 平均控制点误差（像素） | 标准差（像素） | 平均+1标准差（像素） |
|---|---|---|---|
| 无损 | 8.00? | 0.13 | 0.24 | 0.37 |
| 有损 | 6.43 | 0.29 | 0.89 | 1.17 |
| 有损 | 5.04 | 0.25 | 0.64 | 0.89 |
| 有损 | 3.25 | 0.26 | 0.73 | 0.99 |
| 有损 | 2.09 | 0.24 | 0.70 | 0.94 |
| 有损 | 1.58 | 0.22 | 0.59 | 0.81 |
| 有损 | 1.21 | 0.32 | 0.60 | 0.92 |
| 有损 | 1.07 | 0.34 | 0.72 | 1.06 |

?无损未压缩的位数/像素

图14 表3和表4中小波网络测试的加权平均误差（ε）及相应标准差（σε）。如图所示，实心黑点代表案例1，即参考图像和配对控制点图像均进行有损压缩，而空心圆圈为案例2，即参考图像无损压缩。注意，对于此面板，ε以Log??标度绘制。曲线为对测试结果数据的简单经验函数拟合。左图：实线曲线为ε=1.305?1.24x+0.46x2?0.070x3+0.004x?+10?·exp(?12.0x)，其中x为以位数/像素（bpp）为单位的压缩级别，拟合自案例1结果。虚线曲线绘制了经验函数ε=0.105+0.029x+0.549·exp(?x)，拟合自案例2结果。右图：整体趋势的实线曲线参数为a=0.178，b=1.55，c=?0.22。

每种情况的最佳拟合曲线显示，当bpp值达到其范围内最小值时，加权平均误差ε呈现显著增加，正如我们预期的那样，这是由于小尺度表面细节的丢失以及压缩伪影的逐步发展所致。

每条曲线也展现出一个独特的??最小值，均约为2.6 bpp，该值小于最低压缩量（即6.43 bpp）对应的值。这与我们的假设一致：在某种程度下，压缩过程表现为一种低通滤波器，能够去除随机图像噪声。正如预期，Case 2中??随bpp的变化幅度相较于Case 1（其中两张图像均经过有损压缩）相对较小。标准差σ_ε随??增大而增大的总体趋势与压缩量增加时直方图变宽的现象大体一致。正如预期，在我们初步测试的压缩/解压缩图像对中，平均控制点误差随压缩比的增加而增加。在我们的测试中，当码率为5.04 bpp时，平均误差达到1.00，标准差为1.20像素。然而，需要注意的是，可见圆盘上相关性存在显著的系统性变化。特别是，我们考察的全盘图像在正视几何条件下表现明显优于靠近边缘或晨昏线附近的条件。因此，在实际应用中，合理预期解压缩图像在表面选定区域内将满足要求。图15展示了四张被认为适合控制点匹配的阋神星压缩全盘图像，它们实现的压缩码率低至1.58 bpp。

图15下载 | 全尺寸图像 | 使用我们的神经网络编码器（应用于32×32像素子块，圆盘总直径约为750像素）对阋神星全图像（N1477638933）进行压缩和恢复。上行显示压缩至不同码率（从左至右分别为5.04、3.25、2.09和1.58 bpp）的图像。下行显示与原始图像的差异，经对比度拉伸以显示（以黑色标示）有损图像中更亮的像素。图像经过遮罩处理，排除了距离阋神星中心超过盘径5%的天空区域。

4.5 土星环中的螺旋桨结构

在下一项测试中，我们训练网络来压缩和解压缩包含所谓"螺旋桨"结构的土星环图像。56 螺旋桨因其形状而得名（见图16中的示例），是土星环中由小天体引力引起的低对比度特征。我们期望网络能够在不直接针对此类特征进行训练的情况下保留这些特征。图17展示了一个来自12位图像的极微弱螺旋桨示例，以及使用我们的方法将图像压缩至1.53 bpp的结果。我们看到，尽管网络并未明确针对螺旋桨图像进行训练，但它仍然能够恢复螺旋桨。这表明我们的压缩器有能力在太阳系外行星任务过程中支持新的发现。我们还在图18中绘制了机器学习小波变换网络相对于对相同32×32子块应用的db5小波压缩方法的拟合模型参数的均方误差（MSE）。56 在0.6至1.75 bpp的压缩范围内，机器学习方法实现了比原始算法更低的MSE。

图16下载 | 全尺寸图像 | 卡西尼图像N1586628622的土星环上半部分，其中包含一个螺旋桨（黄色方框内）。

图17下载 | 全尺寸图像 | 来自卡西尼（N1586627590）的土星环图像子块示例，包含螺旋桨（左），以及使用我们的神经网络压缩器压缩至1.53 bpp的结果（右）。螺旋桨的两个叶瓣被网络成功恢复。颜色表示缩放后的像素值（DN / 4096），因此值在0到1之间。

图18下载 | 全尺寸图像 | 机器学习小波编码器和原始小波压缩算法的均方误差（MSE），通过对压缩图像拟合螺旋桨模型56并将模型参数与对应无压缩图像的模型参数进行比较来计算，基于一组包含螺旋桨的12位土星环图像（064RI_AZSCAN系列）。在约0.6 bpp至1.75 bpp之间，机器学习方法的MSE低于经典小波变换。在更高压缩码率下，我们可以看到机器学习方法在估计本征噪声水平时达到最低MSE。

图19下载 | 全尺寸图像 | 螺旋桨从标定图像N1586627074重新投影后，按照文献56中的模型进行拟合（绿色椭圆）。原始图像显示在顶部。后续各行分别对应机器学习小波压缩器的0.62、1.26和3.78 bpp（左列）以及文献1中的原始db5小波压缩算法（右列）。

我们可以从图19直观地看到MSE改善如何转化为螺旋桨恢复效果，该图展示了压缩对文献56中描述的几何拟合的影响。在此图中，我们比较了应用于32×32像素子块的机器学习小波压缩和文献1中的原始小波压缩算法。我们看到，在此特定示例中，螺旋桨的指示被机器学习压缩保留至0.62 bpp，而原始算法则无法做到这一点。此外，螺旋桨的形状被机器学习更好地建模至1.26 bpp的码率。在图20中，我们将螺旋桨拟合的若干参数作为码率的函数绘制，分别展示机器学习和原始小波算法的结果。我们可以看到，尽管在极低码率下螺旋桨恢复对机器学习小波压缩器没有返回值，但它比原始算法更接近参数的真实值。这很可能是由于原始算法在选择不重要小波时的偏差所致，也证明机器学习压缩器能够学习到感兴趣的潜在信号，并将更多重要比特用于重建该信号。

图20下载 | 全尺寸图像 | 根据文献56的模型对螺旋桨拟合的四个参数作为码率的函数。文献1中原始小波算法的恢复值以红色（带圆圈）显示，机器学习小波压缩器的恢复值以绿色（带方块）显示。从无压缩图像获得的值由虚线表示。随着码率增加，机器学习压缩器图像产生的恢复参数更接近对应的无压缩值。

5. 多图像压缩：飞越序列

如我们在4.3节中讨论的，高精度重建压缩图像的需求之一是在飞越序列（如阋神星飞越序列）的连续图像之间进行控制点匹配。此过程与视频压缩中的光流估计非常相似，后者会为连续视频帧估计像素到像素的映射。对我们的目的而言，两幅??×??图像之间的光流是一个??×??×2张量，其中(??,??,??)坐标表示(??,??)像素在两张图像之间在水平或垂直方向的偏移，取决于??=1还是2。为了同时压缩多张图像，我们实现并训练了文献14中的SpyNet光流估计网络。该网络如图21所示，是一种低复杂度模型，具有空间金字塔结构，在一定程度上模拟了手工完成的控制点匹配。该网络将一对连续图像缩小并应用卷积神经网络（图中的??0）来估计光流的粗略近似。然后它上采样光流以将第一张图像的较大尺度版本扭曲变换到第二张图像，将两者传入另一个卷积模型（??1）以估计残余光流，并迭代此过程。

图21下载 | 全尺寸图像 | 用于估计光流的来自文献37的空间金字塔网络。类似于控制点匹配例程，该网络缩小连续图像以粗略估计光流。然后它放大光流以将较大尺度的图像相互扭曲变换，根据误差更新流图，并迭代执行。

在实践中，给定来自阋神星飞越序列的两张连续图像，我们对第一张图像在两个维度上执行循环平移，以最大化其与第二张图像的相关性。这使它们的32×32子块的特征紧密对齐。然后我们训练了三个版本的波压缩网络：第一个用于压缩初始图像，第二个用于压缩光流图（每个子块为32×32×2的图），第三个用于压缩将第一子块通过解压缩光流扭曲变换到第二子块后的残差误差。注意，光流压缩器将光流图视为两个图像子块。图22展示了两张连续图像子块的放大图，说明了光流如何将第一图像子块扭曲变换为第二子块几乎完全相同的副本，捕捉到两个子块残差中出现的几乎所有感兴趣特征。图23说明了残差误差压缩器如何捕捉光流补偿后仍留在子块中的细微特征。在图24中，我们展示了使用光流和残差压缩器的全图效果，与单独使用神经网络小波压缩器压缩图像的效果进行比较。我们发现，这一过程获得了更好的图像重建效果，去除了部分单独压缩子块所产生的伪影，并且以每像素更少的比特数实现。

图22下载 | 全尺寸图像 | 左上：阋神星飞越序列的第一图像子块。左下：经全图平移以实现最大相关性后，飞越序列第二图像中对应的图像子块。右上：两个对齐图像子块之间的差异。右下：使用光流网络将第一图像子块扭曲变换到第二子块后的新误差残差。误差大幅降低，几乎没有剩余的信号结构，代价是需要编码光流的两个维度。颜色表示通过将DN减去图像均值再除以图像标准差得到的归一化值。

图23下载 | 全尺寸图像 | 左：光流校正后的初始残差误差（光流压缩至1.2 bpp）。中：估计并压缩/恢复的残差误差（额外1.2 bpp）。右：最终残差误差。

图24下载 | 全尺寸图像 | 左：阋神星全图像（12 bpp）。中：使用神经网络小波压缩器单独压缩图像后的误差残差（2.95 bpp）。右：使用光流从飞越序列中的前一张图像预测图像，并使用单独的神经网络压缩器压缩前一张图像、光流和残差后的误差（每幅图像平均2.69 bpp）。

5.1 光流网络结果

我们在4.4节中使用的同一组阋神星参考和测试图像上验证了此方法。四个光流网络测试的结果如表5和图25所示。表5显示了每张图像的压缩程度（以每像素比特数表示）——如上所述，参考图像通过专用小波网络进行压缩，而测试图像的码率则包括通过两个独立网络压缩光流和误差残差的比特。该表表明光流和误差残差被比参考图像更激进地压缩，如果此方法要应用于飞越序列中更长的图像串，这将是理想的情况。在我们测试过的控制点匹配方法中，此方法最有前景，能在显著压缩图像的同时保留对立体测光应用而言满意的空间和光度保真度。所有压缩级别之间的主要差异发生在被排除控制点的边缘和一段连续的可接受控制点区域（见图25）。

表5 光流网络结果

| 压缩方式 | 参考bpp : 测试bpp | 平均控制点误差（像素） | 标准差（像素） | 均值+1标准差（像素） |
|---|---|---|---|---|
| 无损 | 8.00:8.00a | 0.13 | 0.24 | 0.37 |
| 有损 | 3.27:2.40 | 0.17 | 0.10 | 0.27 |
| 有损 | 2.64:1.82 | 0.15 | 0.09 | 0.24 |
| 有损 | 2.24:1.33 | 0.20 | 0.12 | 0.32 |
| 有损 | 1.90:0.62 | 0.29 | 0.12 | 0.41 |

a 无损未压缩的每像素比特数。

图25下载 | 全尺寸图像 | 光流网络测试结果。上行显示相关性加权像素偏移误差的直方图。下行显示控制点位置以及它们是否在1像素误差内匹配。

6. 使用对抗网络区分阋神星图像中的压缩误差与噪声

最后，我们简要讨论了推广L2范数来度量重建误差的适当方式。该范数在加性高斯噪声情况下是理想的，但一般情况下噪声模型将取决于仪器标定。

为此，我们对土卫五的一组图像样本进行了内在误差的估计，具体方法是：对每张图像施加一个3×3中值滤波核，并测量其与原始像素值的差异。我们拟合了一个模型，描述平均像素噪声随像素值变化的函数关系，该函数本质上是非线性的。我们测试了一种对抗判别器网络的使用效果，该网络由卷积层组成，旨在根据原始图像区分我们的压缩器的重建误差与基于模型的内禀噪声。对抗网络的发展与推广归功于文献[28]。

对于每个图像块，判别器接收原始图像块以及压缩误差或模拟噪声，将它们拼接为一个32×32×2的张量。判别器网络结构简单：首先依次施加两个5×5的二维卷积核，每个卷积核的步长为2，偏置共享（tied biases）。每层的输入以填充数P=2进行零填充，每个卷积核之后接一个leaky ReLU激活函数，负斜率系数为0.3。第一个卷积层有64个滤波器（输出16×16×64的张量），第二个有128个滤波器（输出8×8×128的结果）。这一过程如图26所示。第二层的输出被展平后通过一个带有单一输出和sigmoid激活的全连接层，试图输出一个二值判断，即输入是真噪声还是压缩误差。具体来说，网络的可训练参数总数（包括两个卷积核、全连接层及其偏置）为：(5×5×64+1)+(5×5×128+1)+(8×8×128+1)=12,995。

图26
下载全尺寸图像

测试中使用的判别器网络。将图像块与压缩误差或建模噪声进行拼接。判别器网络施加两个卷积层，后接一个全连接层，以产生一个二值输出，指示输入是压缩误差还是噪声。

在训练过程中，判别器依次查看压缩误差和模拟噪声，对每一个产生二值判断，并训练以最小化其各判断的二值交叉熵的线性组合。而压缩网络则同时训练以最小化判别器关联准确率与压缩系数的微分熵的线性组合。这一过程如图27所示。从图28可以看到，当压缩到1.54 bpp时，该方法成功重建了图像的特征，且重建残差确实与基于我们建模的内禀噪声分布的随机高斯噪声有所偏离。该方法的一个优势在于，训练好的判别器网络可用于标记全图像中压缩误差是否过高，并调用更为保守的压缩策略。

图27
下载全尺寸图像

GAN训练过程：对于每个图像块，我们生成压缩残差误差和模拟噪声，分别输入判别器网络（各自与原始图像拼接）。判别器试图将每一个分类为压缩误差或噪声，目标是最小化其预测的二值交叉熵。同时，编码器/解码器网络训练以最小化判别器关联准确率与微分熵的线性组合。

图28
下载全尺寸图像

左：土卫五的原始图像。中：使用对抗判别器网络训练的神经网络小波变换编码器压缩并重建的版本，压缩至1.54 bpp。右：重建误差。

我们再次使用与第4.4节相同的参考图像和测试图像来验证该方法在共点匹配方面的效果，分别考虑两张图像均由小波网络压缩的情况（案例1）和参考图像无损压缩的情况（案例2）。我们的结果见表6和表7，表中列出了训练好的网络在不同平均压缩比特数/像素下的平均像素误差和标准差。

表6 对抗网络：两张图像均有损压缩（案例1）

| 压缩（比特/像素） | 平均共点误差（像素） | 标准差（像素） | 均值+1标准差（像素） |
|---|---|---|---|
| 无损 | 8.00a | 0.13 | 0.24 | 0.37 |
| 有损 | 8.39 | 0.19 | 0.44 | 0.63 |
| 有损 | 4.31 | 0.19 | 0.49 | 0.68 |
| 有损 | 3.24 | 0.18 | 0.46 | 0.64 |
| 有损 | 2.64 | 0.27 | 0.75 | 1.02 |
| 有损 | 2.24 | 0.29 | 0.85 | 1.14 |
| 有损 | 1.96 | 0.43 | 1.17 | 1.60 |
| 有损 | 1.55 | 0.45 | 1.16 | 1.61 |

a 无损未压缩的比特/像素。

表7 对抗网络：参考图像无损压缩（案例2）

| 压缩（比特/像素） | 平均共点误差（像素） | 标准差（像素） | 均值+1标准差（像素） |
|---|---|---|---|
| 无损 | 8.00a | 0.13 | 0.24 | 0.37 |
| 有损 | 8.39 | 0.15 | 0.22 | 0.37 |
| 有损 | 4.31 | 0.13 | 0.14 | 0.27 |
| 有损 | 3.24 | 0.19 | 0.46 | 0.65 |
| 有损 | 2.64 | 0.14 | 0.16 | 0.30 |
| 有损 | 2.24 | 0.17 | 0.37 | 0.54 |
| 有损 | 1.96 | 0.19 | 0.34 | 0.53 |
| 有损 | 1.55 | 0.15 | 0.22 | 0.37 |

a 无损未压缩的比特/像素。

案例1和案例2的最佳拟合趋势曲线（图29所示）的形状比第4.4节原始小波网络测试结果（图14）更为相似，两者之间的差异主要体现在幅度上。这里，ε的值在所有位置均小于1像素。曲线中的极小值比第4.4节的结果不那么明显，尽管每条曲线在最小的bpp值范围内都表现出预期的上升趋势。

图29
下载全尺寸图像

表6和表7中对抗网络测试的加权平均误差（ε）及对应标准差（σε）。案例1是参考图像和配对共点图像均有损压缩的情况，案例2是参考图像无损压缩的情况。曲线拟合测试数据。左：实线曲线是ε=0.105+0.0104x+1.67·exp(?x)，拟合案例1结果，其中x为每像素比特数（bpp）。虚线曲线是ε=0.126+0.0023x+0.157·exp(?x)，拟合案例2结果。右：总体趋势的实线曲线参数为a=?0.642，b=6.94，c=?6.44。

7. 结论与未来方向

在本工作中，我们着重展示了机器学习在补充、推广和优化经典压缩算法方面的多种应用。此外，我们的模型具有低复杂度和快速可训练的特性，即使使用有限的数据集也表现出色。通过将网络限制在32×32图像块上运行，我们能够在数小时内将其训练至收敛，通常只需要数千甚至数百个充分采样的图像块就能在新图像上实现高质量压缩。这证明了科学与压缩领域专家的协作在选取合适的模型和数据集以产生令人满意的压缩算法方面的重要性和有效性。

当然，还有许多方向值得探索，以促进这一研究的成熟。虽然将注意力限制在小图像子块上允许我们在有限的算力下快速训练网络，但我们可以通过增大网络规模、训练集和训练时间来获得更好的性能。挑战在于如何在土星等外太阳系任务（如卡西尼号）上以较低的成本上传和运行压缩器。

此外，在训练模型时选择要最小化的损失函数方面，还有大量内容值得探索。正如本文所展示的，对信号和仪器噪声的了解应当指导训练过程中所执行的优化。第6节中探讨的对抗判别器网络为一种仅由仪器噪声知识所表征的隐式损失函数提供了有前景的途径，如果我们能构建更精确的仪器噪声模型，它将有望产生更丰富的成果。我们还想测试利用判别器来指导如何向图像子区域分配压缩比特的极限。在理想情况下，从给定科学检索中的目标量到传感器所测量值之间的完整失真函数表征，应当能够同时转化为最优的网络拓扑结构和用于训练该网络的损失函数。

正如我们讨论的，在本研究中我们仅关注了受文献[1]发现启发而设计结构的机器学习变换，我们有理由相信这些结构将很好地适用于我们考虑的科学检索，并且我们能够直接将其与非机器学习方法进行比较。将这些结果与一组已建立的变换编码方法（包括第2.2节提到的机器学习架构，以及如Karhunen-Loève变换等经典方法）的性能进行基准测试将是有益的。此外，我们尚未探索压缩流水线中其他部分的优化。例如，我们在第4.2节中采用的量化方法简单且快速，但未来我们可能会探索经典向量量化方法（如网格编码量化）所带来的增益，或者探索机器学习量化方法[57]。同样，我们尚未对不同熵编码方法（包括经典的算术编码器和机器学习方法）针对量化系数的性能进行基准测试。这些话题留待未来工作。

8. 附录：备注1的证明

若**x**为高斯向量，且我们的编码器函数为单位矩阵**E**，令**y** = ?(**x**) := **E****x**，则**y**也将具有多元高斯统计特性，均值为μ**y** = **E**μ**x**，协方差为**C****y** = **E****C****x****E**?。计算我们损失函数L(ξ|**y**)的期望，根据公式(3)得到：

??**y**(L(ξ|**y**)) = n_e · Σ????? [log(√(2π)σ?) + (σ2**y**? + (μ**y**? ? μ?)2)/(2σ?2)]

其中μ**y**?和σ2**y**?分别为y?的均值和方差，n_e为**y**中的元素数量。当μ? = μ**y**?且σ? = σ**y**?（对所有i）时，该函数关于训练参数μ?和σ?取最小值。代入后我们得到：

min{μ?,σ?} ??**y**(L(ξ|**y**)) = (n_e/2)·log(2πe) + (1/2)·n_e·Σ????? log(σ2**y**?)

设**C****x** = **U****x****Σ****x****U****x**? 和 **C****y** = **U****y****Σ****y****U****y**? 分别为**C****x**和**C****y**的奇异值分解，令diag(**C****y**)、diag(**Σ****y**)和diag(**Σ****x**)分别为**C****y**、**Σ****y**和**Σ****x**对角线元素的向量。注意σ2**y**?是**C****y**的对角线元素，根据Schur-Horn定理，向量diag(**C****y**)被diag(**Σ****y**)所主化。此外，根据Poincaré分离定理（又称Cauchy交错定理），diag(**Σ****y**)被diag(**Σ****x**)所主化。现在，对于任意向量**s** = [s?]?????，定义函数f(**s**) := Σ????? log(s?)。由于对数函数在正实数上严格凹，f(**s**)是Schur凹函数，因此我们有f(diag(**C****y**)) ≥ f(diag(**Σ****x**))。

这表明当**C****y**的对角线元素等于**C****x**的奇异值时，??**y**(L(ξ|**y**))取最小值，这可以通过设置**E** = **U****x**?来实现。如前所述，平方参数σ?2将设置为σ2**y**?，而后者现在等于**C****x**的奇异值。

最后，注意在正交编码器和解码器矩阵**E**和**D**的空间上，根据公式(2)，我们完整损失函数C(**x**)的期望在**D** = **E**?时取最小值，从而使得??(‖**x** ? **D****E****x**‖2) = 0。根据上述讨论，我们得到**D** = **U****x**。此外，由于矢量化参数[μ?]?????将收敛到向量μ**y** = **U****x**?μ**x**，我们有**D**[μ?]????? = μ**x**。