一种基于DFT的土壤表面起伏特征分形数据生成算法的制作方法

本发明涉及一种土壤表面起伏特征分形数据生成算法，具体为一种基于dft的土壤表面起伏特征分形数据生成算法，属于地球物理勘探领域电磁波无损探测方面的数值仿真技术领域。

背景技术：

土壤是指地球表面的一层疏松的物质，土壤颗粒在自然条件下是聚集在一起以土壤结构的形式表现出来。土壤表面起伏特征的复杂程度可用分形维数进行衡量，分形在这个世界是普遍存在的，分型可以归纳出两个特征。第一，分型具有自相似性及分型，自身可以看成是有许多与自己相似的大小不一的部分组成；第二，分形具有无穷多的层次，及无论在分型的哪一个层次，总能看到有更精细的下一个层次存在，或者说分形的图形有无限的细节，可以不断的放大，永远都有结构。分形维数主要描述分形最主要的参量，简称分维。分维反映了复杂形体占有空间的有效性，它是复杂形体不规则性的量度。

真实土壤表面呈现凸凹不平的起伏特征，通常会引起反射电磁波相当程度的抖动。目前，常见的电磁场仿真物理模型，通常采用标准平面模拟复杂起伏的土壤表面。由于表面结构过于简单、理想，难以再现电磁波与真实土壤表面的相互作用过程，仿真结果与实测结果偏差较大，严重影响电磁波无损探测结果的分析和解释。

技术实现要素：

本发明的目的就在于为了解决上述问题而提供的通过构建土壤表面凸凹不平的起伏特征，形成拟真的土壤表面形态的一种基于dft的土壤表面起伏特征分形数据生成算法。

本发明通过以下技术方案来实现上述目的：一种基于dft的土壤表面起伏特征分形数据生成算法，包括以下步骤：

步骤1、生成大小为x×y的随机二维矩阵a；

步骤2、对矩阵a进行二维快速傅里叶变换，并将结果进行平移，使零频率处于二维矩阵的中心；

步骤3、设矩阵a的中心点坐标为(ic,jc)，对于每一个数据点a(i,j)，计算其与中心点的模l；即

步骤4、根据近似幂律关系对矩阵a进行变换；采用公式为：

步骤5、再次对矩阵a进行逆傅里叶变换，将所得结果归一化到[-k,k]的区间范围内，此时矩阵a即代表土壤表面的起伏程度，则该数据可直接用于土壤表面的建模，实现具有分形特征的逼真土壤表面形态模型。

优选的，所述步骤1中，生成二维矩阵a的x和y值皆为400，且生成的矩阵为随机矩阵。

优选的，所述步骤4中，对分形特征数据进行离散傅里叶变换后，其频率f与相应的频率系数mf近似呈幂律关系，即

其中，c为常数，不影响土壤表面起伏特征的衡量；d为分形维度，其值是介于1到3之间的实数。

优选的，所述步骤4中，频域数据经过变换公式处理后能量集中到低频区域。

优选的，所述步骤5中，对矩阵a进行逆傅里叶变换之前需要再次对a进行平移操作。

优选的，所述步骤5中，生成的土壤起伏数据不包含海平面信息，且实际建模时，需要与模型中土壤的基准高度相加，实际的土壤表面高度坐标为：

其中，afractal为模型中土壤起伏曲面所代表的数据，amax为模型中土壤起伏曲面所代表的数据最大值，amin为模型中土壤起伏曲面所代表的数据最小值，hmin为指定模型中该土壤表面高度的最小值，hmax为指定模型中该土壤表面高度的最大值。

本发明的有益效果是：该基于dft的土壤表面起伏特征分形数据生成算法设计合理，可实现真实土壤起伏表面建模，客观再现电磁波与真实土壤表面的相互作用过程，缩小仿真结果与实测结果之间的差异，实现电磁波无损探测结果的准确分析和解释。

附图说明

图1为本发明土壤起伏表面建模实例二维示意图；

图2为本发明土壤起伏表面建模结果三维可视化示意图。

图中，a、随机二维数据，b、随机数据的二维傅里叶变换，c、应用频域变换公式后的结果，d、最后生成的土壤起伏数据。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1～2，一种基于dft的土壤表面起伏特征分形数据生成算法，包括以下步骤：

步骤1、导入算法实现过程需要的模块，并生成大小为x×y的随机二维矩阵a，生成二维矩阵a的x和y值皆为400，且生成的矩阵为随机矩阵，生成符合以上所定义尺寸的随机矩阵a，其数值的取值范围为0到1；

编程方式为：

importnumpyasnp

fromscipyimportfftpack

surfacedims＝(400,400)

r＝np.random.randomstate()

a＝r.randn(surfacedims[0],surfacedims[1])

步骤2、对矩阵a进行二维快速傅里叶变换，并将结果进行平移，使零频率处于二维矩阵的中心；

编程方式为：

a＝fftpack.fftn(a)

a＝fftpack.fftshift(a)

步骤3、设矩阵a的中心点坐标为(ic,jc)，对于每一个数据点a(i,j)，计算其与中心点的模l；即

编程方式为：

步骤4、根据近似幂律关系对矩阵a进行变换，频域数据经过变换公式处理后能量集中到低频区域；采用公式为：

对分形特征数据进行离散傅里叶变换后，其频率f与相应的频率系数mf近似呈幂律关系，即

其中，c为常数，不影响土壤表面起伏特征的衡量；d为分形维度，其值是介于1到3之间的实数。

编程方式为：

a＝fftpack.ifftshift(a)

a＝np.real(fftpack.ifftn(a))

步骤5、再次对矩阵a进行逆傅里叶变换，对矩阵a进行逆傅里叶变换之前需要再次对a进行平移操作，将所得结果归一化到[-k,k]的区间范围内，此时矩阵a即代表土壤表面的起伏程度，则该数据可直接用于土壤表面的建模，实现具有分形特征的逼真土壤表面形态模型。

生成的土壤起伏数据不包含海平面信息，且实际建模时，需要与模型中土壤的基准高度相加，实际的土壤表面高度坐标为：

其中，afractal为模型中土壤起伏曲面所代表的数据，amax为模型中土壤起伏曲面所代表的数据最大值，amin为模型中土壤起伏曲面所代表的数据最小值，hmin为指定模型中该土壤表面高度的最大值，hmax为指定模型中该土壤表面高度的最小值。

编程方式为：

a_min＝a.min()

a_max＝a.max()

a＝a*(2*k/(a_max-a_min))-k+((2*k)/(a_max-a_min))*k

工作原理：在使用该基于dft的土壤表面起伏特征分形数据生成算法时，可实现真实土壤起伏表面建模，客观再现电磁波与真实土壤表面的相互作用过程，缩小仿真结果与实测结果之间的差异，实现电磁波无损探测结果的准确分析和解释。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。