CN109509251A - 多因素三维土石混合体生成方法 - Google Patents

多因素三维土石混合体生成方法 Download PDF

Info

Publication number
CN109509251A
CN109509251A CN201811325980.7A CN201811325980A CN109509251A CN 109509251 A CN109509251 A CN 109509251A CN 201811325980 A CN201811325980 A CN 201811325980A CN 109509251 A CN109509251 A CN 109509251A
Authority
CN
China
Prior art keywords
block stone
stone
block
soil
dimensional
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811325980.7A
Other languages
English (en)
Other versions
CN109509251B (zh
Inventor
杜立兵
刘新荣
蔡白洁
何伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing University
Sichuan College of Architectural Technology
Original Assignee
Chongqing University
Sichuan College of Architectural Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing University, Sichuan College of Architectural Technology filed Critical Chongqing University
Priority to CN201811325980.7A priority Critical patent/CN109509251B/zh
Publication of CN109509251A publication Critical patent/CN109509251A/zh
Application granted granted Critical
Publication of CN109509251B publication Critical patent/CN109509251B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明涉及一种多因素三维土石混合体生成方法，属于形态仿真和数值计算领域。该方法包括步骤：S1：土石混合体信息获取；S2：块石形态库构建，对块石进行分类，并构建块石数字形态库；S3：块石随机生成；S4：多因素定向块体分割，基于等效椭球堆叠，对地层进行多因素定向块体分割。S5：多因素土石混合体生成，采用S3在块体内部随机生成块石，生成同时满足含石量、级配、倾角、长轴比、块石形态的多因素土石混合体模型。S6：将步骤S2‑S5整合到一个软件，提供土石混合体断面图像和块体三维形态信息，就可生成符合多因素土石混合体模型。解决无法生成同时满足含石量、级配、倾角、长轴比，以及块石形态的多因素土石混合体生成问题。

Description

多因素三维土石混合体生成方法
技术领域
本发明属于形态仿真和数值计算领域，涉及多因素三维土石混合体生成方法。
背景技术
土石混合体广泛存在自然界和工程建设中，其性质介于土和碎裂岩体之间，土石混合体是由具有一定工程尺度、高强度的块石，分布于低强度的细粒基质土体中，构成极端不均匀(inhomogeneous)的松散岩土介质系统。土石的阈值为工程尺寸(隧道为隧道洞径，边坡为边坡的坡高)的5％，土、石间极端的物理性质差异和土石混合体空间组成的非均匀特性，使得实际土石混合岩土工程难以采用常规的连续、均质的理论方法和数值方法进行计算和设计。特别现实中大体量的土石混合体岩土工程，由于块石尺寸巨大，土石混合体的物理力学性质通常无法通过常规室内试验直接获得。同时由于现场试验的高成本、长周期，和场地限制等特性，现场试验往往难以开展。从而土石混合体的物理力学参数难以获取，进而影响实际土石混合体岩土工程的安全性、稳定性、经济性。为此，当前土石混合体的研究多采用室内试验结合数值模拟的方式进行，故构建合理的土石混合体模型是研究土石混合体岩土工程问题的先决条件。块石是土石混合体的重要组成部分，块石的形态和空间组成又称为土石混合体的细观结构，精细化土石混合体包括块石的含石量、级配、长轴比、倾角、块石形态等多种因素，这些因素决定着土石混合体的物理力学性质，精细化土石混合体建模是土石混合体研究中不可或缺的部分。
故为了深入准确的研究土石混合体岩土工程问题，亟待开发能够同时反映块石的含石量，块石级配、长轴比、倾角规律，以及块石形态等多因素的精细化建模方法和软件，辅助获取土石混合体的物理力学性质，同时辅助设计验证土石混合体工程的稳定性、安全性，经济性。
现有土石混合体生成方法较少，分析研究中多将土石混合体考虑为均匀介质，或采用数字图像、CT扫描生成真实混合体模型，但数字图像需要出露土石混合体地层，且仅能生成二维土石混合体，CT扫描方法受尺寸大小的限制且成本较高。虽也有部分学者采用随机生成的方法生成随机土石混合体数值模型，但仅能生成多边形块石，形态单一，且只能考虑含石量和级配，无法同时全方位反映块石的含石量、级配、长轴比、倾角、块石形态等多种因素。故生成的土石混合体模型难以准确的反映复杂土石混合体工程变形和受力特性。
1真实土石混合体模型成本高，且模型尺寸受限
真实混合体模型来源于对CT图像或真实土石混合体照片进行数字图像处理得到，由于CT扫描成本高昂，且受尺寸大小限制，同时真实照片需出露土石混合体地层，其所需人工多，操作难度大(或基本不可能)，同时时间成本、经济成本高昂，故真实土石混合体模型难以获得，也难以用于实际应用。
2.随机土石混合体模型块石形态为多边形，无法反映真实块石
当前土石混合体模型多位随机土石混合体模型，但其中块石为多边形或常规圆形、矩形等，难以反映真实的块石表面形态。
3.随机土石混合体无法同时满足实际块石的含石量，块石级配、长轴比、倾角，以及空间分布规律
现有随机土石混合体模型无法同时满足块石的含石量，块石级配、长轴比、倾角，以及空间分布等规律，这将极大限制了生成土石混合体的应用范围，也将造成其与真实土石混合体模型的巨大差异。
4.随机土石混合体生成方法复杂，效率低
且当前方法，采用块石不断加入单一空间的循环迭代算法，随着块体数量增加，块石放置难度和算法执行效率将不断降低，且难以生成高含石量的土石混合体。
发明内容
有鉴于此，本发明的目的在于提供一种多因素三维土石混合体生成方法，是一种能同时满足含石量、级配、长轴比、倾角、块石形态，贴近真实土石混合体的便捷、高效土石混合体生成方法。
为达到上述目的，本发明提供如下技术方案：
多因素三维土石混合体生成方法，该方法包括以下步骤：
S1：块石信息获取；
S2：块石形态库构建；
S3：三维块石随机生成；
S4：多因素块体定向分割；
S5：多因素土石混合体生成；
S6：将步骤S1-S5整合到一个软件，使得基于真实块体形态信息和土石混合体信息生成细观的土石混合体结构模型，且进一步开发，能够生成直接用于计算的土石混合体细观计算模型。
进一步，所述步骤S1为：首先，对土石混合体开挖断面进行数码成像，并对断面图像进行二值化、滤波和图像分割方法，编程快速提取土石混合体的整体信息，包括含石量、块石粒径和长轴比，以及块石的二维形态轮廓和平面分布；随后，对土石混合体试样通过断层CT扫描，获取各个断面土石混合体图像，并通过三维重构，批量得到土石混合体试样内块石的三维形态轮廓；最后，对中、小型块石进行逐一3D扫描获取块石的三维形态轮廓；
具体包括以下内容：
①断面的成像及数字图像处理；首先对土石混合体断面进行高清成像，其次对断面图像采用图像二值化、滤波和图像分割方法，编程快速提取所有块石的二维形态轮廓；最后，根据块石轮廓编程计算块石面积、惯性矩阵、中心、长轴比和倾角几何信息，最终输出土石混合体含石量，块石的粒径、倾角、长轴比分布规律；
②试样的断层CT扫描；首先对中、小型土石混合体试样进行CT断层扫描；其次对获得的断层图片，进行三维重构，建立土石混合体试样的三维重构模型；最后，从三维重构模型批量获取块石的三维轮廓信息；
③块石的3D扫描；首先确定土石的分界阈值，若设阈值为1cm，则表示1cm以上为块石，以下为基质土；其次，对土石混合体进行筛分，将不同粒径下的块石区分，并得到块石粒径级配曲线，最后，对不同粒径下的块石进行除尘去泥，并逐一进行3D扫描获取不同粒径下的块石3维轮廓信息。
进一步，所述步骤S2为：首先，对获取的块石二维形态轮廓、三维形态轮廓根据粒径、长轴比和块石成分进行分类，初步建立真实块石形态库；其次，对所有分类的二维、三维块石进行平移、旋转和缩放块石标准化处理，将块石转化为面积或体积等于1，中心位于原点、主轴与x、y、z轴平行的标准化块石；随后，在球坐标系向下，对标准化三维块石轮廓进行球谐函数转换，得到三维块石轮廓的球谐形态列向量，进而得到三维形态库的球谐函数形态矩阵；最后，在极坐标系下，对标准化二维块石轮廓进行傅立叶变换，得到二维块石轮廓的形态列向量，进而得到二维块石形态库的傅立叶形态矩阵；
具体包括以下内容：
①真实块石形态库构建；对所有获得的二维、三维块石形态轮廓，根据粒径、长轴比和成分因素进行分类，初步建立真实块石库；即将块石按粒径界限10cm、20cm、50cm分为4类，长轴比界限1、1.5、2分为4类，则最终构建了4×4＝16个块石库；
②块石形态标准化，各个块石库内的块石轮廓对进行平移、旋转和缩放标准化处理，输出块石中心与原点重合，长短轴与x、y轴重合，且块石面积、体积等于单位面积、体积的标准块石轮廓，便于后续数字化处理；
③三维球谐函数形态矩阵构建；在球坐标系下，基于球谐函数将单个三维标准块石轮廓转换为形态列向量，首先，如式(1)所示，球坐标系下，三维标准块石轮廓表面节点V(x,y,z)到块石中心点O(0,0,0)距离的球谐函数级数表征为
其中θ∈[0,π]，为表面节点相对中心点的倾向和倾角，通过节点V(x,y,z)进行求解，为形态特征指数，表征整个三维块石的形态轮廓，m∈[-n,n]，因此总共有(n+1)2项，含有(n+1)2个形态指数同时为n阶球谐函数，其如式(2)所示
其中为勒让德函数，通过MATLAB内置的Legendre函数直接求解；研究表明n＝15时，球谐函数精确表征块石形态，此时含有256个形态指数而通常块石3D扫描下，标准块石表面轮廓节点数>2000，对所有节点建立式(1)所示方程，并通过最小二乘法求解形态指数求解的为复数最终形态列向量如式(3)所示
Sαi含有512个变量，对形态库的所有块石进行球谐函数变换，最终构建单个块石库的三维球谐函数形态矩阵SMIJ＝[Sα1 Sα2 … SαN]，I,J分别为形态库粒径、长轴比的编号；
④二维傅立叶形态矩阵构建；在极坐标下，基于傅立叶转换为级数将单个二维块石轮廓转换为形态列向量，首先，极坐标系下，对二维土石混合体标准石轮廓节点V(x,y)到块石中心点O(0,0)距离r(θ)如式(1)所示，
其中An、Bn代表特定形态的傅里叶形态指数，具有唯一性；根据研究An、Bn分别选取前10项精确描述颗粒形态，而二维块石轮廓节点数量>200，则对所有节点建立式(4)所示方程，并通过最小二乘法求解形态指数An、Bn，形成块石的形态列向量：
Pα1＝[A1…An,B1…Bn]' (5)
其次，对整个块石形态库的所有块石形态进行傅里叶级数展开，生成形态库矩阵，即对N个颗粒进行形态描述，得到形态向量矩阵：
PM＝[Pα1 Pα2 … PαN] (6)。
进一步，所述步骤S3为：首先，从构建的块石形态库中，直接随机选择真实形态的标准块石轮廓作为输出；其次，对真实块石形态库的球谐函数形态矩阵进行PCA主成分分析，再基于形态库的真实块石特征，通过主成分向量和随机场进行随机块石生成，最后，若没有三维块石形态矩阵，仅含有二维块石矩阵，则通过二维块石形态库求解形态特征向量，并基于傅立叶逆变换和多变量高斯随机场生成随机三维块石轮廓；
具体包括以下内容：
①块石形态库块石随机选取；若三维块石形态库中块石含有量足够大，即远远大于20，则根据块石编号，直接从块石形态库中随机选择三维块石块石；
②球谐函数块石随机生成；若三维块石形态库中块石含有量不足，则通过已有三维块石库的块石形态特征，通过球谐函数随机生成三维块石；步骤如下：
首先，对三维球谐函数形态矩阵SM(Mij)进行PCA主成分分析；其中SM(Mij)含有512×n个元素，n为形态库的块石个数；过程如下：
1)计算SMIJ(Mij)行的均质u(ui)，u尺寸为512×1，其中
2)计算SM(Mij)中心矩阵
SMcen＝SM-huT (8)
h为尺寸512×1，所有元素为1的列向量；
2)计算SM(Mij)自相关矩阵SMcon，并计算相关矩阵矩阵SMcon的特征向量P和特征值λ
其中特征向量Pi对应的特征值λi越大，特征向量Pi越能体现整个协方差矩阵SMcon的特征，从而选取特征值大的Pi作为分析主成分；根据研究表明，将特征向量Pi根据特征值λi的大小从小到大排列，其前50项的累计贡献为97.2％；故仅需要特征向量P前50项的特征向量Pi和对应特征值λi来表征整个形态矩阵，进而忽略掉其他非关键的特征向量，将形态矩阵的维度由512降为50；得到主成分向量PC50＝[P1 P2 P3 … P50]；
其次，根据主成分向量PC50，特征值λi随机生成新的随机块石；根据式(9)生成随机块石的球谐形态指数尺寸512×1；
其中，uT为块石形态矩阵的行均质列向量，yi为服从标准正太分布的随机数，为第i个特征值的平方根，Pi为对于特征值λi的特征向量；
③高斯随机场块石生成；当没有三维块石形态矩阵，仅含有二维傅立叶块石形态矩阵时，根据二维块石形态库矩阵的标准形态向量；球坐标系下块石轮廓节点V(x,y,z)到块石中心O(0,0,0)的距离为以为变量的二维随机场，通常含有均质、方差，以及概率密度函数；基于中心极限定理，大样本统计分布收敛于高斯分布，认为距离在平面内服从二维高斯正太分布，θ∈[0,π]，从而采用高斯随机场构建三维随机块石；具体流程为
首先，拟定三维随机块石轮廓含有2562个节点，傅立叶形态指数取前128项；将二维形态矩阵PM的均质列向量作为该形态库的标准块石形态参数，对标准二维块石形态指数进行傅立叶逆变换，得到关于角度的一维自相关函数
同理假定角度θ∈[0,π]的一维自相关函数与相同，则θ∈[0,π]的自相关函数为
且认为θ∈[0,π]，相对独立，联合自相关函数假定为
其次，随机生成块石含有2562个节点，各节点存在唯一的平面坐标，则通过节点与节点之间的相位差Δθ,带代入式(13)中，生成节点间的自相关系数；最终建立尺寸2562×2562的自相关矩阵C；
最后，对自相关矩阵C计算特征值λ，特征向量P，并由公式(15)随机生成块石；
其中{R}为所有节点的半径，λi为自相关矩阵[ρ]的特征值，{Φi}为对应的特征向量P，ξi为均质为0，方差为1的正太分布随机数，{μ}为{C}均值。
进一步，所述步骤S4为：首先，块石的堆叠椭球生成，并统计堆叠椭球的长轴比、倾角、级配统计规律；其次，基于开发的椭球快速波前堆叠算法，按照长轴比、倾角、级配规律将生成的堆叠椭球进行高质量、低孔隙的逐一相切堆叠；也能够采用颗粒流软件中球组成椭球clump进行堆叠，生成高质量椭球堆叠；其次，椭球多因素多边形生成，将堆叠椭球进行球分解(颗粒流clump分解)，由分解圆球生成Voronoi加权多面体，最后根据球与椭球的从属关系合并多面体，定多面体分割，生成的多边体将继承大部分的椭球特性，如长轴比、倾角，以及块石粒径，进而生成定多因素的多面体分割；
具体包括以下内容：
①块石的堆叠椭球生成，对获取的单个块石轮廓通过边框六面体——等效椭球——堆叠椭球的步骤进行求取，其中各部分的意义如下：
1.边框六面体，完全包裹块石的最小六面体，用于求解块石的倾角、长短轴、长轴比；
2.等效椭球，内切于边框六面体，且长轴比与六面体相等的椭球，其倾角、长短轴、长轴比与块石形同，用于等效块石；
3.堆叠椭球，基于含石量缩放的等效椭球，其倾角、长短轴、长轴比与块石形同，用于地层的椭球堆叠生成特定含石量的土石混合体；
生成各个块石的堆叠椭球后，统计其倾角、长轴比、级配分布规律，用于后续堆叠；
②堆叠椭球的高密度相切堆叠，首先根据统计规律，随机生成单个椭球，其次，将单个椭球放入堆叠空间，进行逐一相切堆叠；最后，堆叠方法采用椭球波前层层堆叠法进行堆叠，或采用颗粒流软件中椭球clump挤压堆叠；
③多因素块体分割，首先，将堆叠后的椭球进行球体分割，其次，根据加权Voronoi多面体，采用分割后的球体生成初始多面体，最后，根据椭球与分解球的从属关系，将来源于于同一椭球的多面体进行合并，生成分割块体；由于块体其将继承大部分的椭球粒径、倾角；长轴比特性，故生成了定多因素块体分割。
进一步，所述步骤S5为：首先，多边形分类，将多因素多边形根据块石类别(例如长轴比、粒径)进行分类，其次，多边形内块石生成，根据不同的类别，采用步骤S2，在多边形内部生成满足含石量虚拟随机块石；最后，检查块石是否与多边形边界重叠，如不重叠则接受，否则重新生成；最终生成多因素土石混合体模型；
具体包括以下内容：
①多因素块体分类；分割后的块体具有特定属性，因此需要对应特定的块石形态库，故首先对分割块体进行分类，为更好的选择形态库，分类的级别与块石形态库相同；
②多边形内部虚拟块石生成，根据块体分类，对特定块体由分类选择合适块石形态库，采用步骤S2生成虚拟块石；
③虚拟块石重叠性，检查生成的虚拟块石是否与块体边界相交，如不相交则接受，相交则重新生成块石，最终生成多因素三维土石混合体模型。
本发明的有益效果在于：
(1)开发满足真实块石形态特征，且高效的真实块石随机生成方法，以在土石混合体中生成真实形态块石；
(2)开发同时满足真实块石级配、长轴比、倾角，空间分布规律，且高效的块体分割方法，用以保证土石混合体中块石可同时满足级配、长轴比、倾角，空间分布；
(3)开发满足真实块石形态，块石含石量，块石级配、长轴比、倾角，空间分布规律，且高效的多因素三维土石混合体生成方法，扩大随机土石混合体模型应用范围，以更好的贴合实际土石混合体中的块石分布；
(4)开发执行效率高，循环迭代次数少的高效多因素三维土石混合体生成方法，扩大随机土石混合体模型应用效率范围，以更好的服务工程建设；
附图说明
为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：
图1为本发明关系图；
图2为块石数码拍照图；
图3为块石形态轮廓获取图；
图4为块石的堆叠椭球及几何特性求解示意图；
图5为块石的尺寸统计分布图；
图6为块石的长轴比统计分布图；
图7为块石的尺寸统计分布图；
图8为土石混合体断面块石分类图；
图9为标准化块石示意图；
图10为4类块石形态库图；
图11随机块石真实示意图；
图12堆叠椭球堆叠示意图；
图13基于椭球堆叠的定向块体切割示意图；
图14基于椭球堆叠的块体定向切割效果；
图15土石混合体生成示意图；
图16最终土石混合体模型；
图17土石混合体渐进网格划分。
具体实施方式
下面将结合附图，对本发明的优选实施例进行详细的描述。
如图1所示，本发明分为4个部分：
(1)块石信息获取
首先，对土石混合体开挖断面进行数码成像，并对断面图像进行二值化、滤波、图像分割等方法，编程快速提取土石混合体的整体信息(含石量、块石粒径、长轴比等)，以及块石的二维形态轮廓和平面分布；随后，对土石混合体试样通过断层CT扫描，获取各个断面土石混合体图像，并通过三维重构，批量得到土石混合体试样内块石的三维形态轮廓。最后，对中、小型块石进行逐一3D扫描获取块石的三维形态轮廓；
包括以下内容：
①断面的成像及数字图像处理。首先对土石混合体断面进行高清成像，其次对断面图像采用图像二值化、滤波、图像分割等方法，编程快速提取所有块石的二维形态轮廓。最后，根据块石轮廓编程计算块石面积、惯性矩阵、中心、长轴比、倾角等几何信息，最终输出土石混合体含石量，块石的粒径、倾角、长轴比分布规律。
②试样的断层CT扫描。首先对中、小型土石混合体试样进行CT断层扫描。其次对获得的断层图片，进行三维重构，建立土石混合体试样的三维重构模型。最后，从三维重构模型批量获取块石的三维轮廓信息。
③块石的3D扫描。首先确定土石的分界阈值(如假定1cm，即1cm以上为块石，以下为基质土)，其次，对土石混合体进行筛分，将不同粒径下的块石区分，并得到块石粒径级配曲线，最后，对不同粒径下的块石进行除尘去泥，并逐一进行3D扫描获取不同粒径下的块石3维轮廓信息。
(2)块石形态库构建
首先，对获取的块石二维形态轮廓、三维形态轮廓根据粒径、长轴比、块石成分等进行分类，初步建立真实块石形态库。其次，对所有分类的二维、三维块石进行平移、旋转、缩放等块石标准化处理，将块石转化为面积或体积等于1，中心位于原点、主轴与x、y、z轴平行的标准化块石。随后，在球坐标系向下，对标准化三维块石轮廓进行球谐函数转换，得到三维块石轮廓的球谐形态列向量，进而得到三维形态库的球谐函数形态矩阵。最后，在极坐标系下，对标准化二维块石轮廓进行傅立叶变换，得到二维块石轮廓的形态列向量，进而得到二维块石形态库的傅立叶形态矩阵。
包括以下内容：
①真实块石形态库构建。对所有获得的二维、三维块石形态轮廓，根据粒径、长轴比、成分等因素进行分类，初步建立真实块石库。例如将块石按粒径界限10cm、20cm、50cm分为4类，长轴比界限1、1.5、2分为4类，则最终构建了4×4＝16个块石库；
②块石形态标准化，各个块石库内的块石轮廓对进行平移、旋转、缩放等标准化处理，输出块石中心与原点重合，长短轴与x、y轴重合，且块石面积、体积等于单位面积、体积的标准块石轮廓，便于后续数字化处理。
③三维球谐函数形态矩阵构建。在球坐标系下，基于球谐函数将单个三维标准块石轮廓转换为形态列向量，首先，如式(1)所示，球坐标系下，三维标准块石轮廓表面节点V(x,y,z)到块石中心点O(0,0,0)距离的球谐函数级数表征为
其中θ∈[0,π]，为表面节点相对中心点的倾向和倾角，可以通过节点V(x,y,z)进行求解，为形态特征指数，表征整个三维块石的形态轮廓，m∈[-n,n]，因此总共有(n+1)2项，含有(n+1)2个形态指数同时为n阶球谐函数，其如式(2)所示
其中为勒让德函数，可以通过MATLAB内置的Legendre函数直接求解。研究表明n＝15时，球谐函数可以精确表征块石形态，此时含有256个形态指数而通常块石3D扫描下，标准块石表面轮廓节点数>2000，对所有节点建立式(1)所示方程，并通过最小二乘法求解形态指数求解的为复数最终形态列向量如式(3)所示
Sαi含有512个变量，对形态库的所有块石进行球谐函数变换，最终构建单个块石库的三维球谐函数形态矩阵SMIJ＝[Sα1 Sα2 … SαN]，I,J分别为形态库粒径、长轴比的编号。
④二维傅立叶形态矩阵构建。在极坐标下，基于傅立叶转换为级数将单个二维块石轮廓转换为形态列向量，首先，极坐标系下，对二维土石混合体标准石轮廓节点V(x,y)到块石中心点O(0,0)距离r(θ)如式(1)所示，
其中An、Bn代表特定形态的傅里叶形态指数，具有唯一性。根据研究An、Bn分别选取前10项就可以精确描述颗粒形态，而二维块石轮廓节点数量>200，故可以对所有节点建立式(4)所示方程，并通过最小二乘法求解形态指数An、Bn，形成块石的形态列向量：
Pα1＝[A1…An,B1…Bn]' (5)
其次，对整个块石形态库的所有块石形态进行傅里叶级数展开，生成形态库矩阵，如对N个颗粒进行形态描述，得到形态向量矩阵：
PM＝[Pα1 Pα2 … PαN] (6)
(3)三维块石随机生成
首先，可以从构建的块石形态库中，直接随机选择真实形态的标准块石轮廓作为输出。其次，可以对真实块石形态库的球谐函数形态矩阵进行PCA主成分分析，再基于形态库的真实块石特征，通过主成分向量和随机场进行随机块石生成，最后，如没有三维块石形态矩阵，仅含有二维块石矩阵，则可以通过二维块石形态库求解形态特征向量，并基于傅立叶逆变换和多变量高斯随机场生成随机三维块石轮廓。
包括以下内容：
①块石形态库块石随机选取。如三维块石形态库中块石含有量足够大(例如远远大于20)，即可根据块石编号，直接从块石形态库中随机选择三维块石块石。
②球谐函数块石随机生成。如三维块石形态库中块石含有量不足，则可以通过已有三维块石库的块石形态特征，通过球谐函数随机生成三维块石。步骤如下：
首先，对三维球谐函数形态矩阵SM(Mij)进行PCA主成分分析。其中SM(Mij)含有512×n个元素，n为形态库的块石个数。过程如下：
1)计算SMIJ(Mij)行的均质u(ui)，u尺寸为512×1，其中
2)计算SM(Mij)中心矩阵
SMcen＝SM-huT (8)
h为尺寸512×1，所有元素为1的列向量。
2)计算SM(Mij)自相关矩阵SMcon，并计算相关矩阵矩阵SMcon的特征向量P和特征值λ
其中特征向量Pi对应的特征值λi越大，特征向量Pi越能体现整个协方差矩阵SMcon的特征，从而选取特征值大的Pi作为分析主成分。根据研究表明，将特征向量Pi根据特征值λi的大小从小到大排列，其前50项的累计贡献为97.2％。故仅需要特征向量P前50项的特征向量Pi和对应特征值λi来表征整个形态矩阵，进而忽略掉其他非关键的特征向量，将形态矩阵的维度由512降为50。得到主成分向量PC50＝[P1 P2 P3 … P50]。
其次，根据主成分向量PC50，特征值λi随机生成新的随机块石。根据式(9)可以生成随机块石的球谐形态指数尺寸512×1。
其中，uT为块石形态矩阵的行均质列向量，yi为服从标准正太分布的随机数，为第i个特征值的平方根，Pi为对于特征值λi的特征向量。
③高斯随机场块石生成。当没有三维块石形态矩阵，仅含有二维傅立叶块石形态矩阵时，可以根据二维块石形态库矩阵的标准形态向量。球坐标系下块石轮廓节点V(x,y,z)到块石中心O(0,0,0)的距离为以为变量的二维随机场，通常含有均质、方差，以及概率密度函数。基于中心极限定理，大样本统计分布收敛于高斯分布，因此可以认为距离在平面内服从二维高斯正太分布，θ∈[0,π]，从而可以采用高斯随机场构建三维随机块石。具体流程为
首先，拟定三维随机块石轮廓含有2562个节点，傅立叶形态指数取前128项。将二维形态矩阵PM的均质列向量作为该形态库的标准块石形态参数，对标准二维块石形态指数进行傅立叶逆变换，得到关于角度的一维自相关函数
同理假定角度θ∈[0,π]的一维自相关函数与相同，则θ∈[0,π]的自相关函数为
且可以认为θ∈[0,π]，相对独立，联合自相关函数可假定为
其次，随机生成块石含有2562个节点，各节点存在唯一的平面坐标，则可以通过节点与节点之间的相位差Δθ,带代入式(13)中，生成节点间的自相关系数。最终建立尺寸2562×2562的自相关矩阵C。
最后，对自相关矩阵C计算特征值λ，特征向量P，并由公式(15)就可随机生成块石。
其中{R}为所有节点的半径，λi为自相关矩阵[ρ]的特征值，{Φi}为对应的特征向量P，ξi为均质为0，方差为1的正太分布随机数，{μ}为{C}均值。
(4)三维多因素定向分割
首先，块石的堆叠椭球生成，并统计堆叠椭球的长轴比、倾角、级配统计规律。其次，基于开发的椭球快速波前堆叠算法，按照长轴比、倾角、级配规律将生成的堆叠椭球进行高质量、低孔隙的逐一相切堆叠。也可采用颗粒流软件中球组成椭球clump进行堆叠，生成高质量椭球堆叠。其次，椭球多因素多边形生成，将堆叠椭球进行球分解(颗粒流clump分解)，由分解圆球生成Voronoi加权多面体，最后根据球与椭球的从属关系合并多面体，定多面体分割，生成的多边体将继承大部分的椭球特性，如长轴比、倾角，以及块石粒径，进而生成定多因素的多面体分割。
包括以下内容：
①块石的堆叠椭球生成，对获取的单个块石轮廓通过边框六面体——等效椭球——堆叠椭球的步骤进行求取，其中各部分的意义如下：
1.边框六面体，完全包裹块石的最小六面体，用于求解块石的倾角、长短轴、长轴比。
2.等效椭球，内切于边框六面体，且长轴比与六面体相等的椭球，其倾角、长短轴、长轴比与块石形同，用于等效块石。
3.堆叠椭球，基于含石量缩放的等效椭球，其倾角、长短轴、长轴比与块石形同，用于地层的椭球堆叠生成特定含石量的土石混合体。
生成各个块石的堆叠椭球后，统计其倾角、长轴比、级配分布规律，用于后续堆叠。
②堆叠椭球的高密度相切堆叠，首先根据统计规律，随机生成单个椭球，其次，将单个椭球放入堆叠空间，进行逐一相切堆叠。最后，堆叠方法可以采用椭球波前层层堆叠法进行堆叠，或采用颗粒流软件中椭球clump挤压堆叠；
③多因素块体分割，首先，将堆叠后的椭球进行球体分割，其次，根据加权Voronoi多面体，采用分割后的球体生成初始多面体，最后，根据椭球与分解球的从属关系，将来源于于同一椭球的多面体进行合并，生成分割块体。由于块体其将继承大部分的椭球粒径、倾角。长轴比特性，故生成了定多因素块体分割。
(5)多因素三维土石混合体生成
首先，多边形分类，将多因素多边形根据块石类别(例如长轴比、粒径)进行分类，其次，多边形内块石生成，根据不同的类别，采用步骤S2，在多边形内部生成满足含石量虚拟随机块石。最后，检查块石是否与多边形边界重叠，如不重叠则接受，否则重新生成。最终生成多因素土石混合体模型。
包括以下内容：
①多因素块体分类；分割后的块体具有特定属性，因此需要对应特定的块石形态库，故首先对分割块体进行分类，为更好的选择形态库，分类的级别与块石形态库相同。
②多边形内部虚拟块石生成，根据块体分类，对特定块体由分类选择合适块石形态库，采用步骤S2生成虚拟块石；
③虚拟块石重叠性，检查生成的虚拟块石是否与块体边界相交，如不相交则接受，相交则重新生成块石，最终生成多因素三维土石混合体模型。
可达到相同目的的其他替换方案：
第1部分中，对断面数码图片的处理过程中，也可通过人工绘制边界获取，而不用数字图像处理进行二值化、滤波、图像分割等操作。
第3部分中，椭球堆叠也可采用PFC、Yade等颗粒流软件中球组成椭球clump进行堆叠，生成椭球堆叠。与快速椭球堆叠相比，随着颗粒数量的增多，clump堆叠生成效率降低。
第3部分中，堆叠椭球的粒径、长轴比、倾角，也可采用人工自定义的方式给定，用于生成多种因素下的土石混合体，用于敏感性分析。
如图2的堆叠椭球生成示意图所示，对单个块石轮廓分别求取块石：
1.边框六面体，完全包裹块石的最小六面体，可以求解块石的长轴a、短轴b、中轴c，进而求解长轴比λ1＝a/b，中轴比λ2＝c/b并将六面体长轴对应的倾角θ、倾向等同于块石倾角θ、倾向
2.等效椭球，是内切于边框六面体，其长轴比与边框六面体相等的椭球，其倾角、长短轴、长轴比与块石形同，用于等效块石。
3.堆叠椭球，是基于含石量将等效椭球等比缩放的椭球，其倾角、长短轴、长轴比与块石形同，用于地层的椭球堆叠生成特定含石量的土石混合体。
由上获取块石中心点坐标(x,y,z)，块石体积V，块石长轴与x轴的夹角(作为块石倾角)，与z轴的倾角θ，以及块石长轴a、短轴b(作为块石的粒径)、中轴c，以及长轴比λ1＝a/b，中轴比λ2＝c/b。并最后对所有土石混合体中的块石进行处理，得到块石的长轴比、倾角、粒径的统计分布规律。结合含石量和可以得到堆叠椭圆的粒径统计分布规律如图3、图4、图5所示。
由于不同长轴比、不同尺寸的块石形态存在差异，为增加生成块石的真实性，将块石根据长轴比分为3类，分类阈值为1.5、2.0，粒径分为3类，分类阈值为4.5、9.0，总共3×3＝9个类别。
对土石混合体断面采用数字图像处理，可自动进行块石分类，如图6所示。可以看出真实块石未完全包括含有9个类别，仅含有7个块石类别。实际中的三维块石分类可以通过土石混合体筛分进行。
如图7三维块石的标准操作示意图所示，对每一类的所有块石形态轮廓，进行平移、三维旋转、缩放标准化操作，生成单位体积，中心点与原点重合，主惯性轴与x、y、z轴平行的标准块石形态。
如图8最终建立第4类块石形态库。
如图9三维块石形态转换示意图，将三维块石轮廓3维节点V(x,y,z)从xyz坐标系，转换为球坐标系并对块石的所有块石进行操作，并通过式1、2、3建立块石形态库的三维数字形态矩阵
如图10所示的三维块石生成流程所示，通过块石三维形态矩阵生成三维块石的方法有3种：
方法1，当块石形态库丰富，可以直接从块石形态库中，根据块石编号随机选取。
方法2，有三维形态库，但不丰富，可以基于球偕函数对形态库进行进行PCA主成分分析，再进行块石生成，具体实现如下所示。
首先，对三维数字形态矩阵进行PCA主成分分析，其中SM(Mij)含有512×n个元素，n为形态库的块石个数。过程如下：
1)计算SMIJ(Mij)行的均质u(ui)，u尺寸为512×1，其中
2)计算SM(Mij)中心矩阵
SMcen＝SM-huT (8)
h为尺寸512×1，所有元素为1的列向量。
2)计算SM(Mij)自相关矩阵SMcon，并计算相关矩阵矩阵SMcon的特征向量P和特征值λ
其中特征向量Pi对应的特征值λi越大，特征向量Pi越能体现整个协方差矩阵SMcon的特征，从而选取特征值大的Pi作为分析主成分。根据研究表明，将特征向量Pi根据特征值λi的大小从小到大排列，其前50项的累计贡献为97.2％。故仅需要特征向量P前50项的特征向量Pi和对应特征值λi来表征整个形态矩阵，进而忽略掉其他非关键的特征向量，将形态矩阵的维度由512降为50。得到主成分向量PC50＝[P1 P2 P3 … P50]。
其次，根据主成分向量PC50，特征值λi随机生成新的随机块石。根据式(9)可以生成随机块石的球谐形态指数尺寸512×1。
方法3，高斯随机场块石生成，对没有三维形态库，仅有二维形态库时采用，方法如下：
首先，拟定三维随机块石轮廓含有2562个节点，傅立叶形态指数取前128项。将二维形态矩阵PM的均质列向量作为该形态库的标准块石形态参数，对标准二维块石形态指数进行傅立叶逆变换，得到关于角度的一维自相关函数
同理假定角度θ∈[0,π]的一维自相关函数与相同，则θ∈[0,π]的自相关函数为
且可以认为θ∈[0,π]，相对独立，联合自相关函数可假定为
其次，随机生成块石含有2562个节点，各节点存在唯一的平面坐标，则可以通过节点与节点之间的相位差Δθ,带代入式(13)中，生成节点间的自相关系数。最终建立尺寸2562×2562的自相关矩阵C。
最后，对自相关矩阵C计算特征值λ，特征向量P，并由公式(15)就可随机生成块石。
其中{R}为所有节点的半径，λi为自相关矩阵[ρ]的特征值，{Φi}为对应的特征向量P，ξi为均质为0，方差为1的正太分布随机数，{μ}为{C}均值。
如图10的椭圆层层堆叠示意图所示，根据堆叠椭圆的长轴比、倾角、粒径统计规律，进行了单个椭圆的层层堆叠，最终生成了定长轴比、倾角、粒径的高密度椭圆堆叠。
如图11椭球堆叠图，采用图10所示的层层堆叠法，生成了高密度的椭球相切堆叠。
如图12块体分割效果图，根据椭圆堆叠，采用椭圆切割方法，生成了定长轴比、倾角、粒径的块体切割。可以看出最终生成的块体基本继承了块体内部椭圆的长轴比、倾角、粒径等信息。
如图13的土石混合体生成示意图所示，首先，在最终切割的块体中根据块体的长轴比、倾角、粒径分类；其次，根据分类选择对应的块石形态库；然后，有图10所示的随机块石生成方法，生成满足块石形态库特征的随机块石；随后，将生成的随机块石，根据含石量采用缩放、旋转等还原手段，放入块石内部；最后，对所有块体进行块石生成和放置，生成最终的多因素土石混合体模型。
如图14的最终二维土石混合体模型所示，可以看出根据图3、图4、图5的统计信息生成了基本满足实际需要的定长轴比、倾角、粒径土石混合体模型。
如图15～17，对整个土石混合体模型进行了高质量的渐进网格划分，以及对土石混合体模型的FLAC块石网格划分模型。
最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (6)

1.多因素三维土石混合体生成方法，其特征在于：该方法包括以下步骤：
S1：块石信息获取；
S2：块石形态库构建；
S3：三维块石随机生成；
S4：多因素块体定向分割；
S5：多因素土石混合体生成；
S6：将步骤S1-S5整合到一个软件，使得基于真实块体形态信息和土石混合体信息生成细观的土石混合体结构模型，且进一步开发，能够生成直接用于计算的土石混合体细观计算模型。
2.根据权利要求1所述的多因素三维土石混合体生成方法，其特征在于：所述步骤S1为：首先，对土石混合体开挖断面进行数码成像，并对断面图像进行二值化、滤波和图像分割方法，编程快速提取土石混合体的整体信息，包括含石量、块石粒径和长轴比，以及块石的二维形态轮廓和平面分布；随后，对土石混合体试样通过断层CT扫描，获取各个断面土石混合体图像，并通过三维重构，批量得到土石混合体试样内块石的三维形态轮廓；最后，对中、小型块石进行逐一3D扫描获取块石的三维形态轮廓；
具体包括以下内容：
①断面的成像及数字图像处理；首先对土石混合体断面进行高清成像，其次对断面图像采用图像二值化、滤波和图像分割方法，编程快速提取所有块石的二维形态轮廓；最后，根据块石轮廓编程计算块石面积、惯性矩阵、中心、长轴比和倾角几何信息，最终输出土石混合体含石量，块石的粒径、倾角、长轴比分布规律；
②试样的断层CT扫描；首先对中、小型土石混合体试样进行CT断层扫描；其次对获得的断层图片，进行三维重构，建立土石混合体试样的三维重构模型；最后，从三维重构模型批量获取块石的三维轮廓信息；
③块石的3D扫描；首先确定土石的分界阈值，若设阈值为1cm，则表示1cm以上为块石，以下为基质土；其次，对土石混合体进行筛分，将不同粒径下的块石区分，并得到块石粒径级配曲线，最后，对不同粒径下的块石进行除尘去泥，并逐一进行3D扫描获取不同粒径下的块石3维轮廓信息。
3.根据权利要求1所述的多因素三维土石混合体生成方法，其特征在于：所述步骤S2为：首先，对获取的块石二维形态轮廓、三维形态轮廓根据粒径、长轴比和块石成分进行分类，初步建立真实块石形态库；其次，对所有分类的二维、三维块石进行平移、旋转和缩放块石标准化处理，将块石转化为面积或体积等于1，中心位于原点、主轴与x、y、z轴平行的标准化块石；随后，在球坐标系向下，对标准化三维块石轮廓进行球谐函数转换，得到三维块石轮廓的球谐形态列向量，进而得到三维形态库的球谐函数形态矩阵；最后，在极坐标系下，对标准化二维块石轮廓进行傅立叶变换，得到二维块石轮廓的形态列向量，进而得到二维块石形态库的傅立叶形态矩阵；
具体包括以下内容：
①真实块石形态库构建；对所有获得的二维、三维块石形态轮廓，根据粒径、长轴比和成分因素进行分类，初步建立真实块石库；即将块石按粒径界限10cm、20cm、50cm分为4类，长轴比界限1、1.5、2分为4类，则最终构建了4×4＝16个块石库；
②块石形态标准化，各个块石库内的块石轮廓对进行平移、旋转和缩放标准化处理，输出块石中心与原点重合，长短轴与x、y轴重合，且块石面积、体积等于单位面积、体积的标准块石轮廓，便于后续数字化处理；
③三维球谐函数形态矩阵构建；在球坐标系下，基于球谐函数将单个三维标准块石轮廓转换为形态列向量，首先，如式(1)所示，球坐标系下，三维标准块石轮廓表面节点V(x,y,z)到块石中心点O(0,0,0)距离的球谐函数级数表征为
其中θ∈[0,π]，为表面节点相对中心点的倾向和倾角，通过节点V(x,y,z)进行求解，为形态特征指数，表征整个三维块石的形态轮廓，m∈[-n,n]，因此总共有(n+1)2项，含有(n+1)2个形态指数同时为n阶球谐函数，其如式(2)所示
其中为勒让德函数，通过MATLAB内置的Legendre函数直接求解；研究表明n＝15时，球谐函数精确表征块石形态，此时含有256个形态指数而通常块石3D扫描下，标准块石表面轮廓节点数>2000，对所有节点建立式(1)所示方程，并通过最小二乘法求解形态指数求解的为复数最终形态列向量如式(3)所示
Sαi含有512个变量，对形态库的所有块石进行球谐函数变换，最终构建单个块石库的三维球谐函数形态矩阵SMIJ＝[Sα1 Sα2 … SαN]，I,J分别为形态库粒径、长轴比的编号；
④二维傅立叶形态矩阵构建；在极坐标下，基于傅立叶转换为级数将单个二维块石轮廓转换为形态列向量，首先，极坐标系下，对二维土石混合体标准石轮廓节点V(x,y)到块石中心点O(0,0)距离r(θ)如式(1)所示，
其中An、Bn代表特定形态的傅里叶形态指数，具有唯一性；根据研究An、Bn分别选取前10项精确描述颗粒形态，而二维块石轮廓节点数量>200，则对所有节点建立式(4)所示方程，并通过最小二乘法求解形态指数An、Bn，形成块石的形态列向量：
Pα1＝[A1…An,B1…Bn]' (5)
其次，对整个块石形态库的所有块石形态进行傅里叶级数展开，生成形态库矩阵，即对N个颗粒进行形态描述，得到形态向量矩阵：
PM＝[Pα1 Pα2 … PαN] (6)。
4.根据权利要求1所述的多因素三维土石混合体生成方法，其特征在于：所述步骤S3为：首先，从构建的块石形态库中，直接随机选择真实形态的标准块石轮廓作为输出；其次，对真实块石形态库的球谐函数形态矩阵进行PCA主成分分析，再基于形态库的真实块石特征，通过主成分向量和随机场进行随机块石生成，最后，若没有三维块石形态矩阵，仅含有二维块石矩阵，则通过二维块石形态库求解形态特征向量，并基于傅立叶逆变换和多变量高斯随机场生成随机三维块石轮廓；
具体包括以下内容：
①块石形态库块石随机选取；若三维块石形态库中块石含有量足够大，即远远大于20，则根据块石编号，直接从块石形态库中随机选择三维块石块石；
②球谐函数块石随机生成；若三维块石形态库中块石含有量不足，则通过已有三维块石库的块石形态特征，通过球谐函数随机生成三维块石；步骤如下：
首先，对三维球谐函数形态矩阵SM(Mij)进行PCA主成分分析；其中SM(Mij)含有512×n个元素，n为形态库的块石个数；过程如下：
1)计算SMIJ(Mij)行的均质u(ui)，u尺寸为512×1，其中
2)计算SM(Mij)中心矩阵
SMcen＝SM-huT (8)
h为尺寸512×1，所有元素为1的列向量；
2)计算SM(Mij)自相关矩阵SMcon，并计算相关矩阵矩阵SMcon的特征向量P和特征值λ
其中特征向量Pi对应的特征值λi越大，特征向量Pi越能体现整个协方差矩阵SMcon的特征，从而选取特征值大的Pi作为分析主成分；根据研究表明，将特征向量Pi根据特征值λi的大小从小到大排列，其前50项的累计贡献为97.2％；故仅需要特征向量P前50项的特征向量Pi和对应特征值λi来表征整个形态矩阵，进而忽略掉其他非关键的特征向量，将形态矩阵的维度由512降为50；得到主成分向量PC50＝[P1 P2 P3 … P50]；
其次，根据主成分向量PC50，特征值λi随机生成新的随机块石；根据式(9)生成随机块石的球谐形态指数尺寸512×1；
其中，uT为块石形态矩阵的行均质列向量，yi为服从标准正太分布的随机数，为第i个特征值的平方根，Pi为对于特征值λi的特征向量；
③高斯随机场块石生成；当没有三维块石形态矩阵，仅含有二维傅立叶块石形态矩阵时，根据二维块石形态库矩阵的标准形态向量；球坐标系下块石轮廓节点V(x,y,z)到块石中心O(0,0,0)的距离为以为变量的二维随机场，通常含有均质、方差，以及概率密度函数；基于中心极限定理，大样本统计分布收敛于高斯分布，认为距离在平面内服从二维高斯正太分布，θ∈[0,π]，从而采用高斯随机场构建三维随机块石；具体流程为
首先，拟定三维随机块石轮廓含有2562个节点，傅立叶形态指数取前128项；将二维形态矩阵PM的均质列向量作为该形态库的标准块石形态参数，对标准二维块石形态指数进行傅立叶逆变换，得到关于角度的一维自相关函数
同理假定角度θ∈[0,π]的一维自相关函数与相同，则θ∈[0,π]的自相关函数为
且认为θ∈[0,π]，相对独立，联合自相关函数假定为
其次，随机生成块石含有2562个节点，各节点存在唯一的平面坐标，则通过节点与节点之间的相位差Δθ,带代入式(13)中，生成节点间的自相关系数；最终建立尺寸2562×2562的自相关矩阵C；
最后，对自相关矩阵C计算特征值λ，特征向量P，并由公式(15)随机生成块石；
其中{R}为所有节点的半径，λi为自相关矩阵[ρ]的特征值，{Φi}为对应的特征向量P，ξi为均质为0，方差为1的正太分布随机数，{μ}为{C}均值。
5.根据权利要求1所述的多因素三维土石混合体生成方法，其特征在于：所述步骤S4为：首先，块石的堆叠椭球生成，并统计堆叠椭球的长轴比、倾角、级配统计规律；其次，基于开发的椭球快速波前堆叠算法，按照长轴比、倾角、级配规律将生成的堆叠椭球进行高质量、低孔隙的逐一相切堆叠；也能够采用颗粒流软件中球组成椭球clump进行堆叠，生成高质量椭球堆叠；其次，椭球多因素多边形生成，将堆叠椭球进行球分解，包括颗粒流clump分解，由分解圆球生成Voronoi加权多面体，最后根据球与椭球的从属关系合并多面体，定多面体分割，生成的多边体将继承大部分的椭球特性，包括长轴比、倾角，以及块石粒径，进而生成定多因素的多面体分割；
具体包括以下内容：
①块石的堆叠椭球生成，对获取的单个块石轮廓通过边框六面体——等效椭球——堆叠椭球的步骤进行求取，其中各部分的意义如下：
1.边框六面体，完全包裹块石的最小六面体，用于求解块石的倾角、长短轴、长轴比；
2.等效椭球，内切于边框六面体，且长轴比与六面体相等的椭球，其倾角、长短轴、长轴比与块石形同，用于等效块石；
3.堆叠椭球，基于含石量缩放的等效椭球，其倾角、长短轴、长轴比与块石形同，用于地层的椭球堆叠生成特定含石量的土石混合体；
生成各个块石的堆叠椭球后，统计其倾角、长轴比、级配分布规律，用于后续堆叠；
②堆叠椭球的高密度相切堆叠，首先根据统计规律，随机生成单个椭球，其次，将单个椭球放入堆叠空间，进行逐一相切堆叠；最后，堆叠方法采用椭球波前层层堆叠法进行堆叠，或采用颗粒流软件中椭球clump挤压堆叠；
③多因素块体分割，首先，将堆叠后的椭球进行球体分割，其次，根据加权Voronoi多面体，采用分割后的球体生成初始多面体，最后，根据椭球与分解球的从属关系，将来源于于同一椭球的多面体进行合并，生成分割块体；由于块体其将继承大部分的椭球粒径、倾角；长轴比特性，故生成了定多因素块体分割。
6.根据权利要求1所述的多因素三维土石混合体生成方法，其特征在于：所述步骤S5为：首先，多边形分类，将多因素多边形根据块石类别进行分类，包括长轴比和粒径；其次，多边形内块石生成，根据不同的类别，采用步骤S2，在多边形内部生成满足含石量虚拟随机块石；最后，检查块石是否与多边形边界重叠，若不重叠则接受，否则重新生成；最终生成多因素土石混合体模型；
具体包括以下内容：
①多因素块体分类；分割后的块体具有特定属性，因此需要对应特定的块石形态库，故首先对分割块体进行分类，为更好的选择形态库，分类的级别与块石形态库相同；
②多边形内部虚拟块石生成，根据块体分类，对特定块体由分类选择合适块石形态库，采用步骤S2生成虚拟块石；
③虚拟块石重叠性，检查生成的虚拟块石是否与块体边界相交，若不相交则接受，相交则重新生成块石，最终生成多因素三维土石混合体模型。
CN201811325980.7A 2018-11-08 2018-11-08 多因素三维土石混合体生成方法 Active CN109509251B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811325980.7A CN109509251B (zh) 2018-11-08 2018-11-08 多因素三维土石混合体生成方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811325980.7A CN109509251B (zh) 2018-11-08 2018-11-08 多因素三维土石混合体生成方法

Publications (2)

Publication Number Publication Date
CN109509251A true CN109509251A (zh) 2019-03-22
CN109509251B CN109509251B (zh) 2023-02-28

Family

ID=65747822

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811325980.7A Active CN109509251B (zh) 2018-11-08 2018-11-08 多因素三维土石混合体生成方法

Country Status (1)

Country Link
CN (1) CN109509251B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111737858A (zh) * 2020-05-29 2020-10-02 西安理工大学 一种堆石混凝土中随机堆石的生成方法
CN112200754A (zh) * 2020-10-30 2021-01-08 中国矿业大学 一种随机矸石块体三维形状参数自动获取方法
CN112347647A (zh) * 2020-11-06 2021-02-09 西北核技术研究所 一种基于Voronoi图的高含石率土石混合体模型构建方法
CN116644646A (zh) * 2023-05-17 2023-08-25 中山大学 岩石的细观重构方法、装置、设备及存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20110114010A (ko) * 2010-04-12 2011-10-19 전주대학교 산학협력단 경계가 없는 한국 전통 돌담 모델링을 위한 3차원 표면 생성 방법 및 이를 위한 기록매체
CN104036538A (zh) * 2014-06-25 2014-09-10 清华大学 土石混合体三维细观结构重建与分析方法及系统
CN104050717A (zh) * 2014-06-27 2014-09-17 清华大学 土石混合体三维细观结构生成方法及系统
CN104567723A (zh) * 2015-01-28 2015-04-29 清华大学 基于Kinect的在野外环境下物体快速三维扫描方法及装置
AU2015258226A1 (en) * 2009-09-29 2015-12-03 Keystone Retaining Wall Systems, Inc. Wall blocks, veneer panels for wall blocks and method of constructing walls
CN107315009A (zh) * 2017-07-07 2017-11-03 重庆大学 室内土石混合体试样的断层数码成像装置及三维数值逆向重构方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2015258226A1 (en) * 2009-09-29 2015-12-03 Keystone Retaining Wall Systems, Inc. Wall blocks, veneer panels for wall blocks and method of constructing walls
KR20110114010A (ko) * 2010-04-12 2011-10-19 전주대학교 산학협력단 경계가 없는 한국 전통 돌담 모델링을 위한 3차원 표면 생성 방법 및 이를 위한 기록매체
CN104036538A (zh) * 2014-06-25 2014-09-10 清华大学 土石混合体三维细观结构重建与分析方法及系统
CN104050717A (zh) * 2014-06-27 2014-09-17 清华大学 土石混合体三维细观结构生成方法及系统
CN104567723A (zh) * 2015-01-28 2015-04-29 清华大学 基于Kinect的在野外环境下物体快速三维扫描方法及装置
CN107315009A (zh) * 2017-07-07 2017-11-03 重庆大学 室内土石混合体试样的断层数码成像装置及三维数值逆向重构方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
宋岳等: "三峡库区滑坡土石混合体剪切强度试验研究", 《三峡大学学报(自然科学版)》 *
徐文杰等: "土石混合体细观结构及力学特性数值模拟研究", 《岩石力学与工程学报》 *
徐文杰等: "土石混合体随机细观结构生成系统的研发及其细观结构力学数值试验研究", 《岩石力学与工程学报》 *
徐文杰等: "基于真实块石形态的土石混合体细观力学三维数值直剪试验研究", 《岩石力学与工程学报》 *
李世海等: "三维离散元土石混合体随机计算模型及单向加载试验数值模拟", 《岩土工程学报》 *
李亮等: "土石混合边坡块石随机生成方法与稳定性分析", 《湖南大学学报(自然科学版)》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111737858A (zh) * 2020-05-29 2020-10-02 西安理工大学 一种堆石混凝土中随机堆石的生成方法
CN111737858B (zh) * 2020-05-29 2024-03-08 西安理工大学 一种堆石混凝土中随机堆石的生成方法
CN112200754A (zh) * 2020-10-30 2021-01-08 中国矿业大学 一种随机矸石块体三维形状参数自动获取方法
CN112347647A (zh) * 2020-11-06 2021-02-09 西北核技术研究所 一种基于Voronoi图的高含石率土石混合体模型构建方法
CN112347647B (zh) * 2020-11-06 2023-10-20 西北核技术研究所 一种基于Voronoi图的高含石率土石混合体模型构建方法
CN116644646A (zh) * 2023-05-17 2023-08-25 中山大学 岩石的细观重构方法、装置、设备及存储介质
CN116644646B (zh) * 2023-05-17 2024-03-22 中山大学 岩石的细观重构方法、装置、设备及存储介质

Also Published As

Publication number Publication date
CN109509251B (zh) 2023-02-28

Similar Documents

Publication Publication Date Title
CN109509251A (zh) 多因素三维土石混合体生成方法
Latham et al. Three-dimensional particle shape acquisition and use of shape library for DEM and FEM/DEM simulation
Maréchal Advances in octree-based all-hexahedral mesh generation: handling sharp features
Peytavie et al. Arches: a framework for modeling complex terrains
Hanaor et al. 3D printable geomaterials
Cavagna et al. The STARFLAG handbook on collective animal behaviour: Part II, three-dimensional analysis
CN109241646B (zh) 基于椭圆堆叠和随机场的多因素二维土石混合体生成方法
Jagnow et al. Stereological techniques for solid textures
Wang et al. Lidar data segmentation and classification based on octree structure
CN108280290A (zh) 一种混凝土骨料数值模型重建方法
CN106599053B (zh) 一种三维模型检索方法
WO2006058292A9 (en) Methods for retrieving shapes and drawings
CN111524229A (zh) 岩石颗粒的三维几何形貌信息提取系统及方法
CN101201845A (zh) 基于中轴点集分层螺旋信息的三维模型检索方法
Birdal et al. A point sampling algorithm for 3d matching of irregular geometries
Van De Weygaert et al. Alpha shape topology of the cosmic web
Zhou 3D urban modeling from city-scale aerial LiDAR data
Xie et al. Geometric modeling of Rosa roxburghii fruit based on three-dimensional point cloud reconstruction
CN106560814A (zh) 一种基于离散元的三维形态特征集料生成方法
Riley et al. Void space and scaffold analysis of packed particles: applications in granular biomaterials
Tanaka et al. Accuracy-based sampling and reconstruction with adaptive grid for parallel hierarchical tetrahedrization
Johnen Indirect quadrangular mesh generation and validation of curved finite elements
Aragón-Calvo et al. Geometry of the cosmic web: Minkowski functionals from the delaunay tessellation
Bhunre et al. Topological analysis of voxelized objects by discrete geodesic reeb graph
Pranav Persistent holes in the Universe: a hierarchical topology of the cosmic mass distribution

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant