多重分形数值计算技术及其应用

在上面的讨论中，针对康托尔集上的质量分布为例，引入了多重分形的定义。根据这个定义，原则上讲，也可以计算出多重分形谱来。只要计算每个观测尺度下α的数值，以及对应的f（α）就足够了。但是，这种计算往往会带来很大的误差，这主要是因为在此计算中，实际上是用单个盒子的质量分布来推测不同质量分布随着尺度而变的斜率这一整体信息。

下面介绍的方法不是从局部出发，而是从整体出发，计算出在不同观测精度下的Renyi信息维χ。其中，引入了一个新的参数q。在不同的参数q下，能够得到Renyi维随着观测精度变化的速率χ（q）。计算χ（q）可以按照最小二乘法拟合出log（S（q））对log（b（t））的斜率。之后，根据这个速度χ（q）和多重分形谱f（α）构成了勒让德变换，对这一性质，就能估算出多重分形谱f（α）。

（1）Renyi信息维

首先要为多重分形定义一个新的量，叫做Renyi信息维度。对于一个多重分形场来说，当给定一个观测精度t，就得到了每一个小盒子b（t）中的质量分布mi，其中要求 ，那么，对于这些mi整体，我们可以计算：

Renyi信息维就是这个量当t趋于无穷大时的极限，即定义为：

其中q是一个参数，可以取遍所有实数。n（t）是在t精度下所有小盒子的个数。这个Renyi信息维与分形维的定义十分相似，它计算的是广义的信息量

随着尺度增长的速度，当q=0的时候，求和式的结果就等于n（t），于是这个维度χ（0）就退化成了多重分形的支撑集的分形维。当q=1的时候，这个量计算的是Shannon信息熵随尺度变化的速度。在实际的计算中，通过为待计算的集合划分不同精度的盒子，然后可以很容易地计算出随log（b（t））的变化曲线，在对数坐标下，这条线变成了一条直线，于是可以计算出该直线的斜率，即为χ（q），它反映的是多重分形场的整体信息，因为它对整个场的质量分布进行求和。q的作用相当于一种权重，如果q＞0，则那些质量大的小盒子就会对χ的贡献越大，从而χ就反映出来那些质量密集区域随着尺度变化信息量的增长速度。反过来若q＜0，则mi越小的区域会对整体的χ贡献越大，于是χ反映的就是质量小区域的信息量增长速度。

（2）Legendre变换

由于χ（q）和f（α）存在着天然的联系。下面做一定的推演。首先，在给定精度t下，根据χt（q）的定义：

其中n（t）表示t精度下小盒子的个数。在这些盒子中，有很多盒子具有相同的质量，即对于每个mi都有N（mi）个小盒子，这就意味着上式中的和式可以写为：

也就是说，可以把对不同小盒子i的求和转换成对不同的盒子质量mi的求和。而根据前面给出的α的定义，有：

同样的道理，根据f（α）的定义：

于是，代到和式中：

相当于把求和下标又变成了α。当t趋于无穷大的时候，这个和式b（t）会远远小于1。于是求和号中的最小的那一项会起到最主要的贡献，以至于可以用下面的近似式：

这样，把求和号去掉了，相应地变成了使得等式中最小的那个α所对应项的值。把这些结果代入到Renyi维数的定义中就有：

这样，函数χ（q）和f（α）就构成了一个勒让德（Legendre）变换对。按照Legendre变换的性质，也可以根据函数χ（q）而反过来求出f（α）。(www.xing528.com)

于是，便可以通过全局的Renyi维度计算出多重分形谱f（α）。

在具体的计算中，可以针对一个α值，跑遍所有不同的q值，从而得到使得αq-χ（q）最小的一个q，从而计算出一个f（α），但是这样计算起来过于烦琐。事实上，只要注意到，由于要让αq-χ（q）最小，可以得到一个α关于最优的q的响应函数的α（q）。再把最优解代回去，又会得到f（α（q）），也就是说无论α还是f，实际上都是q的函数。这样，再让q跑遍所有的可能值，就会得到一组α（q）和f（q）的点对，把这一组点对画在α-f平面上就得到了想要的多重分形谱。而根据χ（q）的定义，实际上能够具体地写出α（q）和f（q）。由于最优的q要让αq-χ（q）最小。对q求导，得到 α=χ'（q），而根据定义，有：

这样，就得到了α（q）的表达式：

其中，。这样，可以求出f（q）的表达式：

如此，就可以画出f（q）和α（q）的关系图。

以上为相关的数学推导。接下来不妨把多重分形谱的计算看作一个函数multifractal，它的输入是一个n*n矩阵M，代表一个二维的场，矩阵的每个元素代表对应场点的质量。当然至于更高维度的场也可以做类似的推广。那么multifractal计算需要如下几个步骤：

①初始化，计算出n*n的空间总共能分成多少个方格盒子。一般地，设最小的盒子边长为2，最大的盒子边长为2int（log2n）；

②对所有不同边长的盒子循环：

a.对每个盒子，计算盒子内的总质量mij；

b.将质量归一化，pij=mij/Σmij。

③对不同的q进行循环（事先给q一个预设的取值范围）：

a.再次对所有边长的盒子（分辨率）循环，计算在该分辨率下的A=Σμilog mi和B=μi log μi

b.根据不同分辨率下的A，B，对盒子大小b（t）在双对数坐标下做回归，计算出来α（q）和χ（q）

④输出α（q）和χ（q）的坐标组合，并画图。