一维RMD是一种常用的分形曲线生成方法，它通过迭代的方式产生具有自相似性特征的曲线。计算时先将待插值的两点间连一条直线，把直线的中点向上或向下移动一个随机量，作为一个插值点，从而形成两条线段，以此来代替原线段；继续分别对两线段做上述操作，如此下去，就可以得到复杂的、满足自相似性的分形曲线^[13]。其通用插值迭代公式为：

式中，P_mid为插值点；P_i、P_i+1表示线段的两个端点; Δn为中点位移量，其大小与自相似性参数(或粗糙度因子)H、方向控制参数τ、偏移量控制参数δ以及迭代次数n有关，关系式为：

式中，H的取值范围为[0, 1]，H值越大，生成的分形曲线越光滑; H值越小，生成的分形曲线越粗糙；偏移量控制参数δ可以设置为一个固定值，也可以用随机函数或正态分布的高斯函数来代替。图 1给出了一维RMD的示意图。

图  1  一维RMD示意图

Figure 1.  Schematic of One-Dimensional RMD

从式(1)、式(2)可知，一维RMD主要通过移位种子点(线段中点)、移位方向、移位偏移量和迭代次数来控制生成的图形形态。在对海岸线进行分形插值时，由于传统的一维RMD通常将海岸线要素视为普通的几何线对象，以几何线上的点作为图形的特征点，移位种子点始终为线段的中点，移位方向随机选择正上方或正下方，偏移量是个固定值或随机数，使得插值得到的结果图形虽然能够模拟海岸线的分形特征，但无法顾及海岸线的弯曲结构特征，也难以满足空间数据的质量要求。下文以海岸线分形插值示例来具体说明。

如图 2所示，曲线a为比例尺1:3 000万的中国广西沿海某段的海岸线，采用一维RMD进行3次分形迭代插值，得到结果如图 2中曲线b所示。从图 2中可以看出，曲线b虽然在视觉形态上能够模拟海岸线的分形特征，但其弯曲的方向和坐标偏移程度较原始曲线a出现了较大差异，在同一个弯曲内，插值得到的曲线方向随机地与弯曲方向同向或相反，削弱了海岸线的弯曲结构特征。另外，由于是对曲线进行整体均衡分形插值，无法体现不同岸线地貌类型的插值结果差异，难以正确反映海岸线的真实地理特征。因此，本文基于海岸线的地理弯曲特征保持考虑，提出一种顾及海岸线地理弯曲特征约束的可控分形插值方法，以满足空间数据的质量要求。

图  2  一维RMD对海岸线进行分形插值

Figure 2.  Fractal Interpolation Based on One-Dimensional RMD for Coastline

曲线插值就是在不够精细的地方加入描述细节的点，并尽量反映曲线的精细特征^[13]。对于海岸线而言，其精细程度主要取决于地图比例尺，比例尺越大，细节描述越详细，越能够准确地反映实际地理状况。当需要对小比例尺海岸线数据进行分形插值来重建相应大比例尺数据时，如何在插值过程中同时保持海岸线的分形特征和弯曲特征，是当前分形插值面临的一个难题。另外，海岸线要素只在一定的尺度范围内才具有分形特征，这个范围就是无标度区，只有在无标度区范围内分形特征才具有稳定性^[15]。由于无标度区的计算方法不是本文的研究重点，在此不作过多阐述。不过，已有研究表明，常见地图比例尺海岸线要素都在无标度区范围内^[16]，本文主要针对这类常规比例尺海岸线要素进行分形插值研究，如果起止尺度海岸线要素超出了无标度区范围，则分形插值无意义，不在本文考虑范围之内。

本文的基本思想是：在无标度区范围内，首先根据不同海岸地貌类型所呈现的不同弯曲特征和分形特征，对海岸线进行地貌单元划分，将传统的整体分形插值变换为以海岸线地貌弯曲特征为划分单元的分段插值组合；其次，利用一维随机中点移位法对各划分单元分别进行插值，并结合各划分单元的弯曲特征对分形参量分别进行约束控制，以保持海岸线不同地貌单元的弯曲特征；最后，将各分段插值结果顺次连接起来得到新的插值曲线。

海岸在构造运动、海水动力、生物作用和气候因素等共同作用下形成了各种海岸地貌，用海岸线进行抽象表达时，呈现出不同的弯曲特征。海岸地貌根据其地貌的基本特征，可分为侵蚀地貌和堆积地貌两大类。侵蚀地貌是岩石海岸在波浪、潮流等不断侵蚀下所形成的各种地貌，通常地势险峻，其岸线破碎曲折、几何形态复杂，分形特征明显。堆积地貌是近岸物质在波浪、潮流和风的搬运下沉积形成的各种地貌，通常地势平坦、岸线平直，分形特征不明显。因此，不同地貌类型的海岸线可通过其呈现的几何弯曲形态和复杂度来识别，形态复杂、弯曲度大的海岸线通常表示海蚀地貌单元，形态平直、弯曲度小的海岸线通常表示海积地貌单元。

海岸线的弯曲特征是由所处区域的海岸地貌特征决定的。由上述海岸线地理弯曲特征分析可知，整条海岸线由海蚀地貌单元和海积地貌单元组合而成。为了实现不同类型地貌单元的差异化插值，首先需要将海岸线按照弯曲特征进行划分。弯曲识别有多种方法^[17-18]，但从海岸线的地理特征来看，海湾被定义为海岸带向陆地凹进的海域部分。为了与海湾的空间认知一致，本文采用文献[18]中地理曲线弯曲识别方法对海岸线中的海湾进行识别。如图 3所示，根据文献[18]的方法，对图 3中海岸线进行弯曲单元提取后的有序集为{C_{2, 7}, C_{8, 10}, C_{11, 18}, C_{18, 22}}。

图  3  海岸线弯曲特征提取

Figure 3.  Bending Characteristic Extraction for Coastline

在海岸线弯曲特征划分的基础上，可利用各弯曲的特征对一维RMD的主要分形参量分别进行约束控制。从式(2)中可以看出，一维RMD的插值过程主要通过移位种子点、移位方向、移位偏移量和分形次数来控制生成的图形形态。可以结合海岸线要素的弯曲单元特征，对这些参量分别进行约束控制，以实现几何层面的可控分形插值。

以某一弯曲单元为例来阐述这一分形控制过程。假设图 4中的曲线L_i为海岸线某一弯曲单元，A和B是L_i的两个端点，P是曲线L_i上距离线段AB最远的点，O是P点到AB的垂足。线段AB称为弯曲单元的弦，P点与弦AB形成夹角θ称为分形偏移角，PO称为弯曲的深度。从图 4中可以看出，在弯曲弦长相同的情况下，θ角越大，弯曲深度越深，表明弯曲越复杂，通常对应海蚀地貌；θ角越小，弯曲深度越浅，表明弯曲越平直，通常对应海积地貌。因此，海岸线的弯曲地貌单元特征可以通过分形偏移角θ来表达。

图  4  弯曲单元示意图

Figure 4.  Illustration of Bending Unit

分形参量的约束控制方法如下。

1) 在移位种子点上进行控制。一维RMD每次以线段的中点为种子点进行移位加点，插值得到的图形与弯曲单元的形状特征不相符。如图 4弯曲单元中的P点，可以看作是线段AB以O点作为移位种子点的一次分形插值，但O点并不是AB的中点。因此，可选择O点分弦AB的比例作为弯曲单元上每条线段的移位种子点所确定的定比分点比例，从而确定该移位种子点的坐标。设O点分线段AB的比例为λ，弯曲单元上相邻两个顶点的坐标分别为(x₁，y₁)和(x₂，y₂)，则该线段上移位种子点的坐标为$\left( {\frac{{{x_1} + \mathit{\lambda }{\mathit{x}_{\rm{2}}}}}{{1 + \mathit{\lambda }}}, \frac{{{y_1} + \mathit{\lambda }{\mathit{y}_{\rm{2}}}}}{{1 + \mathit{\lambda }}}} \right) $。

2) 在偏移量上进行控制。一维RMD的偏移量参数δ是个固定值或随机数，产生的偏移量不可控。由于分形插值的过程可以视为曲线化简的逆过程，在进行曲线化简时，当某一弯曲深度小于设定的阈值时，该弯曲被删除，否则被保留，即所有弯曲单元中弯曲深度的最小值可作为该尺度下线化简的临界阈值；反过来，在进行分形插值时，最小弯曲深度则可视为插值坐标的最大偏移量。因此，可以选择所有弯曲单元中弯曲深度的最小值作为分形插值的最大偏移量。

3) 在移位方向上进行控制。一维RMD的移位方向始终固定于y轴方向随机向上或向下，无法正确地反映海岸线变化的各向同性的特点。且当曲线上部分相邻的坐标点在横坐标值比较接近时，沿y方向随机移位插值很容易产生曲线自相交。本文将移位方向改为垂直于待插值线段且随机分布于线段的某一侧，从而减少曲线自相交的情况，且有利于保持曲线的分形特征和各向同性的特点。

4) 对分形次数进行控制。由于分形图形是通过迭代的方式产生的，迭代的次数决定了分形的粒度，迭代次数越多，分形粒度越细，对应比例尺越大。因此，可以利用迭代次数表征分形粒度。其计算公式推导如下。

假设对线段AB进行规则的一维RMD分形插值(如图 5所示)，且移位方向垂直于待插值的边。设原始线段AB的长度为l₀，经过第一次分形插值产生的移位点C与AB的夹角为θ(即分形偏移角)，且以后每次迭代产生的移位点与对应边所构成的夹角均为θ，假设第i次分形插值后任意相邻两个坐标点之间的距离为l_i，则对于第i+1次分形插值，其插值后的边长与前一次的边长满足如下关系：

根据规则分形的特点，经过第i+1次分形插值后，共产生了2ⁱ⁺¹条边，此时，该曲线的总长为：

对基于分形分析的曲线插值，曲线长度的变化不是线性的，而是随比例尺变化，即推广的Beckett公式^[19]：

式中，D为曲线的分维数，其值可以通过D-P分维估值法计算得到^[13]；L_M1和L_M2分别为插值前后的曲线长度；M₁和M₂分别为插值前后的曲线比例尺分母。对于曲线插值，满足M₁>M₂和L_M1 < L_M2。

图  5  一维RMD规则分形插值示意图

Figure 5.  Regular Fractal Interpolation Based on One-Dimensional RMD

假设曲线由比例尺1:M₁插值变换到1:M₂需要进行n次分形迭代插值，根据式(4)，该曲线进行n次分形插值后的总长度为：

由于L_M1=l₀，因此，

故

由于迭代插值次数为整数，故需要对式(6)进行取整运算，即：

从式(7)可以看出，在插值前后比例尺一定的情况下，不同类型地貌弯曲单元的分形插值次数不同(因θ不同)。这与海岸线各弯曲单元具有不同的分形特征相符，说明本文方法能够很好地顾及海岸线的地理弯曲特征和分形特征。

当海岸线各弯曲单元实施分形插值后，剩下一些长短不一的非弯曲线段部分。对于这部分线段，仍可以利用一维RMD进行分形插值，即移位种子点选取每条线段的中点，但移位方向选择垂直于线段随机分布于线段的一侧，移位偏移量设置为弯曲深度的最小值，插值次数设定为同一目标尺度下所有弯曲的插值次数的最小值。待所有弯曲部分和非弯曲部分实现分形插值后，按照曲线的方向将各分段的插值结果顺次连接起来，得到最终插值曲线。

为验证本文方法的有效性和合理性，选取中国海南岛三亚市一段长度约45 km、比例尺为1:200万的海岸线要素作为实验数据。插值目标比例尺设定为1:100万和1:50万，并设计对比实验，将本文方法与传统一维RMD和文献[13]的方法进行比较，以及与真实的1:100万和1:50万同要素海岸线数据作形态和精度的对比分析。首先，对海岸线进行弯曲单元划分处理。如图 6所示，该海岸线被划分为8个弯曲单元{C₁，C₂，C₃，C₄，C₅，C₆，C₇，C₈}。

图  6  海岸线分形插值预处理

Figure 6.  Pretreatment of Fractal Interpolation for Coastline

其次，对每一个弯曲单元，按照本文提出的参量控制方法计算相应的分形参量值。其中，利用式(7)计算每个弯曲单元的插值次数较复杂，需要先确定曲线的分维数D和各弯曲单元的分形偏移角θ。分维数D可以通过D-P分维估值法计算得到。对于分形偏移角θ，可以根据图 5中规则分形的方法计算其近似值，即先计算每个弯曲单元的弯曲深度，获取其中的最小值d_min作为各弯曲单元的最大移位偏移量，然后针对各弯曲单元，计算其构成的线段平均长度l作为规则分形的起始边长，由规则分形中边长与偏移量的关系可求$\mathit{\theta = }{\rm{ta}}{{\rm{n}}^{ - 1}}\left( {\frac{{2{d_{\min }}}}{l}} \right) $。

表 1给出了不同目标尺度下各弯曲单元分形插值次数的计算结果。从表 1中可以看出，分形次数n不仅与目标比例尺有关，还与弯曲单元的特征有关。在同一弯曲单元中(如弯曲单元C₁)，目标比例尺越大，分形插值次数越多；在同一目标比例尺下，弯曲单元的线段平均长度越长，分形插值次数越大。

接着，通过对插值过程中分形参量进行控制，得到目标尺度的分形插值结果如图 7所示。其中，图 7(a)曲线是比例尺为1:200万原始海岸线数据；图 7(b)、7(c)曲线是采用本文插值方法对图 7(a)进行目标比例尺分别为1:100万和1:50万的分形插值结果；图 7(d)曲线是采用传统一维RMD对图 7(a)曲线进行3次分形迭代插值的结果；图 7(e)、7(f)曲线是采用文献[13]的方法对图 7(a)曲线进行目标比例尺分别为1:100万和1:50万的分形插值结果。从图 7(b)、7(c)可以看出，目标比例尺越大，插值得到的细节越多。与采用传统一维RMD的结果图相比，本文方法插值得到的结果显然更符合真实海岸线要素的形态特征。与文献[13]的插值方法结果图 7(e)、7(f)相比，由于本文算法的插值次数不仅能根据目标尺度的差异进行自适应调整，而且还能兼顾海岸线各段弯曲特征的差异进行自适应调整，从而使得结果图中各弯曲单元具有不全相同的插值细节，能更好地区分相邻尺度间插值结果的差异，也更符合海岸线要素的随机分形特征。另外，由于在插值前对海岸线进行了弯曲特征预划分处理，使得各分段插值结果不影响海岸线的整体弯曲层次结构，很好地保持了海岸线的地理弯曲特征。相反，由于文献[13]方法是对海岸线进行整体均衡分形插值，各弯曲的插值次数均等，使得部分不需要插值或少插值的弯曲单元也进行了相同的迭代次数插值(即过度插值)，导致目标尺度一定范围内的插值结果相同(如1:80万至1:120万区间范围内的插值结果与1:100万插值次数相同)，造成目标尺度的区分度较低。而且，因为存在过度插值，间接导致了海岸线出现弯曲特征弱化和自相交现象的产生，如图 7(e)、7(f)所示。因此，本文插值方法对海岸线的弯曲形态保持较文献[13]方法更优，也更符合真实的海岸线要素随机分形特征。

另外，为了检验本文的插值方法是否符合空间数据的质量要求，进一步将插值结果分别与真实的1:100万和1:50万同要素海岸线数据作形态和精度的对比分析。精度评估采用缓冲区方法，即对真实的同比例尺海岸线要素作一定精度半径的缓冲区分析，并检查插值结果是否落在缓冲区范围内。由于纸质地图图面精度要求优于0.1 mm，因此，可将其作为精度评价的缓冲区半径。在1:100万和1:50万的地图中，0.1 mm对应实际距离分别为100 m和50 m，分别对1:100万和1:50万的真实海岸线数据作半径为100 m和50 m的缓冲区分析，并将其与本文的插值结果和文献[13]的插值结果进行叠加显示，得到结果如图 8所示。图 8(a)和图 8(b)的曲线比例尺分别为1:100万和1:50万，其中，天蓝色的曲线为对应尺度下真实的海岸线数据，灰色带状区域为真实海岸线缓冲区，黑色曲线和红色曲线分别为本文算法和文献[13]算法的插值结果。从图 8(a)中可以看出，在曲线形态上，本文方法的插值结果与真实同比例尺海岸线要素较接近，比文献[13]的方法插值结果形态好；在插值精度上，本文方法的插值点绝大部分都落在缓冲区范围内，也比文献[13]方法的插值结果精度好，基本上能够满足该比例尺下空间数据的精度要求。但从图 8(b)中可以看出，虽然在形态上本文方法的插值结果仍优于文献[13]的插值结果，但与真实的1:50万海岸线要素相比，无论是本文方法还是文献[13]方法，插值结果与真实的海岸线要素在部分弯曲段出现了较大偏差。这说明当目标尺度跨度较大时(此时的目标尺度跨度为4倍)，两种插值方法的结果精度都难以满足空间数据的质量要求。

图  7  不同分形插值方法结果对比

Figure 7.  Comparison of Interpolation Results of Different Fractal Interpolation Methods

图  8  基于1:200万的插值结果与真实数据对比

Figure 8.  Comparison of Fractal Interpolation Results Based on 1:2 000 000 Scale with Real Data

为了进一步验证本文方法在较小尺度跨度范围内的插值结果质量的可靠性，继续对1:100万真实海岸线进行可控分形插值实验, 得到目标尺度1:50万的曲线数据，如图 9所示。图 9中同样对1:50万的真实海岸线数据作了半径为50 m的缓冲区，黑色曲线和红色曲线分别为本文方法和文献[13]方法的插值结果。从图 9中可以看出，本文的插值方法结果与真实的1:50万海岸线数据接近，插值坐标基本上都在缓冲区范围内，说明坐标精度能够满足当前尺度下的精度要求。多次插值实验表明，当插值目标比例尺与待插值数据原始比例尺跨度不超过3倍时，插值结果精度基本上能够满足空间数据的质量要求。另外，无论从曲线形态上还是精度上看，本文的插值方法结果都比文献[13]的插值方法结果要好。

图  9  基于1:100万的插值结果与真实数据对比

Figure 9.  Comparison of Fractal Interpolation Results Based on 1:1 000 000 Scale with Real Data

最后，从海岸线的分形特征保持来看，通过对插值结果海岸线进行定量的分维数计算可知，插值前海岸线的分维数为1.034，使用本文方法所得1:100万和1:50万的结果分维数分别为1.065和1.083，都与插值前的海岸线分维数相近，说明本文方法都能保持海岸线的分形特征。

本文针对海岸线分形插值缺乏地理弯曲特征约束和分形过程不可控的局限，提出了一种可控的分形插值方法。该方法的特点是根据不同岸线地貌类型所呈现的不同弯曲特征对海岸线进行地貌单元划分，将传统的整体分形插值变换为以岸线地貌弯曲特征为划分单元的差异化分段插值组合；进一步，利用一维RMD对各划分单元分别进行分形插值，并结合各划分单元的弯曲特征对分形参量分别进行约束控制，以实现海岸线不同地貌单元的弯曲特征约束。实际的海岸线要素插值实验表明，本文提出的可控分形插值方法计算简单，能够很好地保持海岸线要素的地理弯曲特征和分形特征，且当目标比例尺与原始比例尺跨度不超过3倍时，插值结果数据基本上能够满足目标尺度下空间数据的质量要求, 可为矢量曲线数据连续多尺度表达、在线渐进式传输和低分辨率矢量曲线数据进行高分辨率重建等方面提供技术支撑。相对已有的分形插值方法，本文提出的改进插值方法虽然能够减少插值后曲线的自相交情况发生，但当目标尺度跨度较大时，仍无法完全避免自相交的产生，下一步将继续对该方法进行改进和完善。