CN107330569A - 基于模拟退火的静力触探土层自动识别方法 - Google Patents
基于模拟退火的静力触探土层自动识别方法 Download PDFInfo
- Publication number
- CN107330569A CN107330569A CN201710707758.2A CN201710707758A CN107330569A CN 107330569 A CN107330569 A CN 107330569A CN 201710707758 A CN201710707758 A CN 201710707758A CN 107330569 A CN107330569 A CN 107330569A
- Authority
- CN
- China
- Prior art keywords
- soil layer
- soil
- layer number
- boundary depth
- max
- 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.)
- Pending
Links
- 239000002689 soil Substances 0.000 title claims abstract description 287
- 238000000034 method Methods 0.000 title claims abstract description 27
- 238000002922 simulated annealing Methods 0.000 title claims abstract description 18
- 230000003068 static effect Effects 0.000 title claims abstract description 15
- 238000011835 investigation Methods 0.000 claims abstract description 20
- 238000013459 approach Methods 0.000 claims abstract description 7
- 239000013598 vector Substances 0.000 claims description 25
- 239000011159 matrix material Substances 0.000 claims description 13
- 238000012797 qualification Methods 0.000 claims description 8
- 238000005457 optimization Methods 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 239000000463 material Substances 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 2
- 238000005259 measurement Methods 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 abstract description 5
- 238000013461 design Methods 0.000 abstract description 4
- 230000007547 defect Effects 0.000 abstract 1
- 238000012360 testing method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000010606 normalization Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000035515 penetration Effects 0.000 description 3
- 238000005553 drilling Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 241001269238 Data Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 238000013524 data verification Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 239000005445 natural material Substances 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000000638 stimulation Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
- G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
- G06F18/24155—Bayesian classification
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Economics (AREA)
- Human Resources & Organizations (AREA)
- Strategic Management (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Business, Economics & Management (AREA)
- Data Mining & Analysis (AREA)
- Tourism & Hospitality (AREA)
- Marketing (AREA)
- Health & Medical Sciences (AREA)
- Probability & Statistics with Applications (AREA)
- Game Theory and Decision Science (AREA)
- Development Economics (AREA)
- Agronomy & Crop Science (AREA)
- Animal Husbandry (AREA)
- Marine Sciences & Fisheries (AREA)
- Mining & Mineral Resources (AREA)
- Operations Research (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Entrepreneurship & Innovation (AREA)
- Quality & Reliability (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于模拟退火的静力触探土层自动识别方法，属于岩土工程勘察领域。本发明在考虑土体空间变异性条件下，根据贝叶斯理论基于一套静力触探数据建立给定土层数目条件下土层边界深度后验分布和土层数目后验分布的贝叶斯方程；采用模拟退火算法自动识别出给定土层数目条件下的最可能土层；根据贝叶斯模型选择方法选出最可能土层数目；采用拉普拉斯逼近方法得到土层边界深度的标准差，进而定量表征土层边界深度的不确定性。本发明解决了现有技术存在的无法在考虑土体空间变异性条件下根据CPT数据划分土层，土层划分结果可靠性未知等缺陷；本发明更加科学、严谨、合理可行，为岩土工程分析和设计提供了科学依据。
Description
技术领域
本发明涉及一种基于模拟退火的静力触探土层自动识别方法，属于岩土工程勘察领域。
背景技术
根据勘探资料划分土层是岩土工程勘察中一项重要任务，为岩土工程分析与设计提供必要的土层剖面信息，包括土层数目和相邻土层的界面位置。静力触探试验(conepenetration test,CPT)因具有良好测试连续性、精度和重现性等优势，在我国岩土工程勘察领域被广泛应用，其主要应用之一就是根据静力触探数据划分土层。
然而，静力触探试验不提供钻探土样，若要更加精确地划分土层需要进行额外的钻探试验，这就增加了勘察成本。目前在工程实践中，仅根据探测数据划分土层，依赖于工程师的工程经验和主观判断。土体材料是一种天然材料，具有空间变异性，测得的有限数据具有离散性，同时，相邻土层界面的物理组成复杂，CPT测试数据受到上覆土层和下卧土层的影响，交叉影响深度因相邻土层土体类型和力学性质差异性而异，导致影响区域内测试结果不准确。以上等因素导致所划分的土层不可避免的具有不确定性。为合理地考虑不确定性对土层划分的影响，目前已有模糊集法、聚类分析法、贝叶斯系统识别方法和小波分析法等。以上方法可根据CPT数据确定土层边界深度，但均无法表征划分土层的不确定性。综上所述，对于给定的CPT数据，合理地考虑空间变异性对土层划分的影响，量化土层边界深度的不确定性，能够有效地反映所得土层剖面的可靠性，为岩土工程分析和设计提供参考。如何在考虑土体材料空间变异性条件下根据CPT数据划分土层，并定量表征土层边界深度的不确定性仍不为所知。
发明内容
本发明的目的在于克服现有技术存在的缺点，寻求一种在考虑土体材料空间变异性条件下基于静力触探数据的土层自动识别方法，并定量表征土层边界深度的不确定性。
为了达到上述目的，本发明提供一种基于模拟退火的静力触探土层自动识别方法，包括如下步骤：
步骤一：由现场的岩土勘察报告获取一套CPT数据，用向量ξ＝[ξ 1,ξ 2,...,ξ N]T表示，ξ为整体CPT数据向量，ξ 1,ξ 2,...,ξ N表示第1,2,…,N层土内的CPT数据，N为土层数目；
步骤二：确定描述CPT数据的空间变异性的概率模型MP(·)及模型参数θ：采用随机场或随机变量概率模型描述土体的空间变异性，假设不同土层之间相互独立；
步骤三：基于贝叶斯理论确定在给定土层数目和CPT数据条件下土层边界深度的后验分布的贝叶斯公式：
P(D N|ξ,N)＝KNP(ξ|D N,N)P(D N|N) 公式一
式中：P(D N|ξ,N)表示给定土层数目和CPT数据条件下土层边界深度后验分布；
KN为归一化常数，与土层边界深度无关；
P(ξ|D N,N)表示给定土层数目和边界深度条件下CPT数据的概率密度函数，称作似然函数；
P(D N|N)为土层边界深度的先验分布，假设整体边界深度向量DN在空间Ω＝{0<D1<D2,…,<DN-1<H}上均匀分布，H为总的探测深度；
D N＝[D1,D2,…,DN-1]为第1至N层土的下边界深度组成的整体边界向量，D1,D2,…,DN-1表示第1,2,…,N-1层土的下边界深度，第N层土的下边界深度DN是确定的，其值等于总的探测深度H；ξ＝[ξ 1,ξ 2,...,ξ N]T表示第1至N层土内的CPT数据组成的整体CPT数据向量；
根据全概率理论和不同土层之间相互独立的假设，似然函数P(ξ|D N,N)可以写成：
式中：P(ξ n|D N,N)表示给定土层数目和边界深度条件下第n层土的CPT数据的概率密度函数；
P(ξ n|θ n,D N,N)为给定一套模型参数θ n条件下根据MP(θ n)建立的ξ n的联合概率密度函数；
P(θ n|D N,N)为模型参数的先验分布，根据岩土勘察报告和文献资料确定；
MP(·)表示步骤二中确定的概率模型；ξ n表示第n层土的CPT数据组成的向量；θ n表示第n层土的概率模型参数；n＝1,2…,N；
步骤四：基于贝叶斯理论确定在已知CPT数据条件下土层数目的后验分布的贝叶斯公式：
P(N|ξ)＝P(ξ|N)P(N)/P(ξ) 公式三
式中：N表示土层数目，N＝1,2,…,Nmax；Nmax表示可能的最大土层数目，根据岩土勘察报告和CPT数据进行预判；
P(N|ξ)表示土层数目后验分布；
P(ξ|N)为给定土层数目条件下CPT数据的概率密度函数，称作模型证据；
P(N)是土层数目的先验概率，在先验信息不充足的情况下，假设每个可能的N值具有相同的先验概率，即P(N)＝1/Nmax；
P(ξ)为归一化常数，与土层数目无关；
步骤五：确定给定土层数目条件下最可能土层边界深度D N,MPV，具体实现如下：
当N＝1时，最可能土层边界深度D N,MPV即为总的探测深度H；
当N≠1时，则需要采用模拟退火优化算法，模拟退火优化算法中优化的目标函数为公式一中土层边界深度后验分布的负对数即fobj＝-ln[P(D N|ξ,N)]，使目标函数fobj的值最小的土层边界深度为最可能的土层边界深度D N,MPV；
步骤六：获取每种可能土层数目的模型证据值，具体实现如下：
当N＝1时，模型证据值等于似然函数值；
当N≠1时，根据拉普拉斯逼近方法估计每种可能土层数目的模型证据值的计算公式：
P(ξ|N)≈P(ξ|D N,MPV,N)P(D N,MPV|N)(2π)(N-1)/2|detJ N(D N,MPV)|-1/2 公式四
式中：N表示土层数目，N＝2,…,Nmax，Nmax表示可能的最大土层数目；
D N,MPV表示土层数目为N时最可能的土层边界深度；
P(ξ|D N,MPV,N)表示土层数目为N、土层边界深度为D N,MPV条件下ξ的概率密度函数；
P(D N,MPV|N)表示土层数目为N条件下土层边界深度的先验分布；
detJ N(D N,MPV)表示土层数目为N条件下，目标函数fobj对最可能土层边界深度D N,MPV二阶偏导的Hessian矩阵的行列式的值；
Hessian矩阵J N的第i行第j列的元素是i,j＝1,2,…,(N-1)；fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
若N≠Nmax，将土层数目N改为N+1重复步骤五、六；若N＝Nmax，进入下面步骤；
步骤七：确定最可能土层数目N*，具体实现如下：
根据贝叶斯模型选择方法，比较已计算出的Nmax个模型证据值，选出最大模型证据值对应的土层数目即为最可能土层数目N*；
步骤八：最可能土层数目N*条件下土层边界深度的不确定性表征，具体实现如下：
根据拉普拉斯逼近方法，最可能土层数目N*条件下土层边界深度后验分布的协方差矩阵计算公式：
上式中：Hessian矩阵的第i行第j列的元素是
为最可能土层数目N*条件下土层边界深度后验分布的协方差矩阵；
表示土层数目为N*时最可能的土层边界深度；
fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
协方差矩阵的对角线即为土层边界深度的方差，进而计算出土层边界深度的标准差表征土层边界深度不确定性。
上述步骤一中获取的一套CPT数据为任意一套根据CPT测量的数据或由CPT测量数据计算得到的数据。
上述步骤四中可能的最大土层数目Nmax，根据模型证据值随N的变化趋势判断所取可能的最大土层数目Nmax是否足够大，若模型证据值随N的变化无下降趋势，则Nmax取值不够大，此时Nmax改为Nmax+1，使N＝Nmax+1，重复步骤五、六、七；若模型证据值随N的变化有下降趋势，则说明Nmax取值足够大。
本发明的有益效果是：
与现有技术相比，基于模拟退火的静力触探土层自动识别方法，在合理地考虑土体的空间变异性条件下确定最可能的土层数目，识别最可能的土层边界深度，并定量表征土层边界深度的不确定性。通过一套模拟的土体分类指数Ic数据验证表明，本发明所提出的方法更加科学、严谨、合理可行，不仅能够正确地识别出土层，还能够有效地反映所得土层剖面的可靠性，为岩土工程分析和设计提供参考。
附图说明
图1是本发明的流程框图。
图2是模拟的土体分类指数数据剖面示意图。
图3是描述土体分类指数空间变异性的随机场概率模型示意图。
图4是归一化土质分类图示意图。
图5是模拟数据土层划分结果示意图。
图6是不同土层数目的模型证据对数值结果示意图。
图7是土层边界深度标准差示意图。
具体实施方式
下面结合附图对本发明作进一步说明。
如图1所示，本发明所述的基于模拟退火的静力触探土层自动识别方法，本发明实施的具体过程为：
步骤一：由现场的岩土勘察报告获取一套CPT数据，用向量ξ＝[ξ 1,ξ 2,...,ξ N]T表示，ξ为整体CPT数据向量，ξ 1,ξ 2,...,ξ N表示第1,2,…,N层土内的CPT数据，N为土层数目；上述获取的一套CPT数据为任意一套根据CPT测量的数据或由CPT测量数据计算得到的数据，如锥尖阻力qc、侧壁摩擦阻力fs、孔隙水压力u2和土体分类指数Ic数据等。
步骤二：确定描述CPT数据的空间变异性的概率模型及模型参数：采用随机场或随机变量概率模型描述土体的空间变异性，确定的概率模型用MP(·)表示，确定概率模型MP(·)的模型参数θ，并假设不同土层之间相互独立；
步骤三：给定土层数目和CPT数据条件下土层边界深度后验分布的贝叶斯公式如公式(1)：
P(D N|ξ,N)＝KNP(ξ|D N,N)P(D N|N) (1)
上式中：P(D N|ξ,N)表示给定土层数目和CPT数据条件下土层边界深度后验分布；
D N＝[D1,D2,…,DN-1]为第1至N层土的下边界深度组成的整体边界向量，D1,D2,…,DN-1表示第1,2,…,N-1层土的下边界深度，第N层土的下边界深度D N是确定的；ξ＝[ξ 1,ξ 2,...,ξ N]T表示第1至N层土内的CPT数据组成的整体CPT数据向量；
P(D N|N)为土层边界深度的先验分布，(获得方法：在没有任何有关土层的先验信息条件下，可以认为当给定的土层数目条件下，所有可能的边界组合是等可能的。假设为均匀分布)反映了获取CPT数据之前关于土层边界深度的先验信息，假设整体边界向量D N在空间Ω＝{0<D1<D2,…,<DN-1<H}上均匀分布，H为总的探测深度；
P(ξ|D N,N)为似然函数，表示给定土层数目和边界条件下CPT数据的概率密度函数，反映了CPT数据提供的关于土层边界深度的信息；
KN为归一化常数(KN＝P(ξ|N)即为步骤五中提到的模型证据事先是一个未知的常数)，与土层边界深度无关。根据全概率理论和不同土层之间相互独立的假设，似然函数P(ξ|D N,N)可以写成如公式(2)：
上式中：
P(ξ n|θ n,D N,N)为给定一套模型参数θ n条件下根据MP(θ n)建立的ξ n的联合概率密度函数；
P(θ n|D N,N)为模型参数的先验分布，反映了关于θ n的先验信息，根据岩土勘察报告和文献资料确定；MP(·)表示步骤二中确定的概率模型；ξ n表示第n层土内的CPT数据组成的向量；θ n表示第n层土的概率模型参数；n＝1,2…,N；
步骤四：土层数目后验分布的贝叶斯公式如公式(3)：
P(N|ξ)＝P(ξ|N)P(N)/P(ξ) (3)
上式中：P(N|ξ)表示土层数目后验分布；N＝1,2,…,Nmax；
P(ξ|N)为给定N条件下ξ的概率密度函数，称作模型证据值，反映了ξ提供的关于含有N层土的分层模型的信息；P(N)是N的先验概率(是一个未知常数，求解过程中不需要知道其取值)，在先验信息不充足的情况下，可假设每个可能的N值具有相同的先验概率1/Nmax，即N在1到Nmax之间服从离散的均匀分布，根据岩土勘察报告和CPT数据进行预判，确定可能的最大土层数目Nmax；P(ξ)为归一化常数，与N无关。
步骤五：确定给定土层数目条件下最可能土层边界深度D N,MPV：
当假设只有一个土层N＝1，即CPT钻探深度内没有其他土层时，确定的边界即为探测深度位置；当N≠1，即假设探测深度范围内存在其他土层时，则需要采用模拟退火优化算法，模拟退火优化算法中优化的目标函数为公式(1)中边界后验分布的负对数即fobj＝-ln[P(D N|ξ,N)]，使后验概率即P(D N|ξ,N)的值最大即使目标函数fobj的值最小的土层边界深度为最可能的土层边界深度D N,MPV；
D N＝[D1,D2,…,DN-1]，D N为第1至N-1土层下边界深度组成的整体边界向量；MPV表示最可能值(most probable value)，D N,MPV表示最可能的土层边界深度。
步骤六：获取每种可能土层数目的模型证据值：
当土层数目N＝1时，模型证据值等于似然函数值；
当土层数目N≠1时，根据拉普拉斯逼近方法估计每种可能土层数目的模型证据值的计算公式如公式(4)：
P(ξ|N)≈P(ξ|D N,MPV,N)P(D N,MPV|N)(2π)(N-1)/2|detJ N(D N,MPV)|-1/2 (4)
上式中：N表示土层数目，N＝2,…,Nmax，Nmax表示可能的最大土层数目；D N,MPV表示土层数目为N时最可能的土层边界深度；
P(ξ|D N,MPV,N)表示土层数目为N、土层边界深度为D N,MPV条件下ξ的概率密度函数；
P(D N,MPV|N)表示土层数目为N条件下土层边界深度的先验分布；
detJ N(D N,MPV)表示土层数目为N条件下，目标函数fobj在最可能土层边界深度D N,MPV的二阶偏导的Hessian矩阵的行列式的值；
Hessian矩阵J N的第i行第j列的元素是i,j＝1,2,…,(N-1)；fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
若N≠Nmax，将土层数目N改为(N+1)重复步骤五、六；若N＝Nmax，进入下面步骤；
步骤七：确定最可能土层数目N*：
根据贝叶斯模型选择方法，根据比较已计算出的Nmax个模型证据值，选出最大值对应的土层数目即为最可能土层数目。
步骤八：最可能土层数目N*条件下土层边界深度的不确定性表征：
根据拉普拉斯逼近方法，则最可能土层数目条件下土层边界深度后验分布的协方差计算公式如公式(5)：
上式中：Hessian矩阵的第i行第j列的元素是
为最可能土层数目N*条件下土层边界深度后验分布的协方差矩阵；
表示土层数目为N*时最可能的土层边界深度；
fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
协方差矩阵的对角线即为土层边界深度的方差，进而计算出土层边界深度的标准差表征土层边界深度不确定性。
下面结合模拟土体分类指数数据进行实例说明。
图3为一套模拟的土体分类指数数据，该套数据模拟参数如表1所示。该虚拟场地由5层土构成，土层厚度分别为2m、3m、10m、20m和15m，对应的4个土层边界深度分别为2m、5m、15m和35m。每个土层采用对数正态随机场模拟土体分类指数Ic的空间变异性，随机场参数见表1。
表1.模拟土层厚度以及每层Ic的统计特征
本发明的具体实施过程如下：
步骤一、土体分类指数Ic的确定：
根据模拟参数如表(1)模拟一套具有五个土层的CPT数据如图2所示土体分类指数Ic。
步骤二、描述土体空间变异性的概率模型及模型参数的确定：采用对数正态随机场模型描述Ic的空间变异性，假设不同土层之间相互独立，对应的随机场模型参数θ为Ic的均值μ、Ic的标准差σ、lnIc波动范围λ，如图3所示。自相关系数计算采用指数型相关函数公式如公式(6)：
ρn＝exp(-2|ΔZ|/λn) (6)
上式中|ΔZ|为不同测点之间的距离；λn为第n层波动范围。
步骤三、给定土层数目条件下土层边界深度后验分布的贝叶斯公式的确定：如公式(1)所示，似然函数可以写成如公式(2)所示，其中ξ n可视作ξ的正态随机场的一次实现，则P(ξ n|θ n,D N,N)表达为如公式(7)：
上式中ξ＝[ξ 1,ξ 2,...,ξ N]T表示一套Ic数据，其中ξ n＝[ξn,1＝lnIc,1,ξn,2＝lnIc,2,...,ξn,kn＝lnIc,kn]T为第n层kn个测点Ic的对数值；θ n＝[μn，σn，λn]；和分别为第n层土中ξ的平均值和标准差，它们均为μn和σn的函数；为ξ n的协方差矩阵；R n为ξ n的相关系数矩阵，其元素根据式(6)计算。
模型参数先验分布的确定：
在现有信息不足时模型参数的先验分布假设为联合均匀分布如公式(8)：
上式中μn,max,σn,max和λn,max分别为μn,σn和λn的最大值；μn,min,σn,min和λn,min分别为μn,σn和λn的最小值。μn,min和μn,max分别取Robertson和Wride提出的由归一化锥尖阻力(Qtn)和归一化摩阻比(Fr)构成的归一化土质分类图如图4所示，上左上角角点(Qtn的最大值和Fr的最小值)和右下角角点(Qtn的最小值和Fr的最大值)对应的Ic值，由Ic计算公式可知，即Qtn-Fr归一化土质分类图上Ic的最小值0.52和最大值4.12；根据定义，标准差为非负值，因此σn,min取为0，σn,max取Ic服从0.52到4.12之间均匀分布时对应的标准差，即λn,min和λn,max可根据文献确定为0.1和1.2。
步骤四：土层数目后验分布的贝叶斯公式的确定：如公式(3)所示。
可能的最大土层数目的确定：假设可能的最大土层数目Nmax＝9。
步骤五、给定土层数目条件下最可能土层边界深度的确定：
当N＝1时，即CPT钻探深度内没有分层，唯一的边界深度即为探测深度位置。当N>1时，即假设探测深度范围内存在其他土层，则需要采用模拟退火优化算法，在给定土层数目条件下搜索使目标函数(fobj＝-ln[P(D N|ξ,N)])最小时的土层边界深度，即为最可能土层边界深度。优化结果如图5和表2所示。
步骤六、最可能土层数目的确定：
当土层数目N＝1时，模型证据值等于似然函数值；当土层数目N>1时，根据模型证据近似求解公式(4)得到不同土层数目模型的模型证据对数值，如表2所示。
步骤七：确定最可能土层数目N*：
如图6可知，N＝5时，lnP(ξ|N)最大，P(ξ|N)也最大，因此本算例中最可能的土层数目N*＝5。当N>5时，lnP(ξ|N)呈减小趋势，说明无需进一步增大N值，说明Nmax足够大。
步骤八、土层边界深度的不确定性表征：
根据公式(5)计算出土层边界深度后验分布的协方差矩阵，进而得到土层边界深度的标准差量化土层边界深度的不确定性如图7。
表2.贝叶斯模型选择与土层划分结果
Claims (3)
1.一种基于模拟退火的静力触探土层自动识别方法，其特征在于：包括如下步骤：
步骤一：由现场的岩土勘察报告获取一套CPT数据，用向量ξ＝[ξ 1,ξ 2,...,ξ N]T表示，ξ为整体CPT数据向量，ξ 1,ξ 2,...,ξ N表示第1,2,…,N层土内的CPT数据，N为土层数目；
步骤二：确定描述CPT数据的空间变异性的概率模型MP(·)及模型参数θ：采用随机场或随机变量概率模型描述土体的空间变异性，假设不同土层之间相互独立；
步骤三：基于贝叶斯理论确定在给定土层数目和CPT数据条件下土层边界深度的后验分布的贝叶斯公式：
P(D N|ξ,N)＝KNP(ξ|D N,N)P(D N|N) 公式一
式中：P(D N|ξ,N)表示给定土层数目和CPT数据条件下土层边界深度后验分布；
KN为归一化常数，与土层边界深度无关；
P(ξ|D N,N)表示给定土层数目和边界深度条件下CPT数据的概率密度函数，称作似然函数；
P(D N|N)为土层边界深度的先验分布，假设整体边界深度向量D N在空间Ω＝{0<D1<D2,…,<DN-1<H}上均匀分布，H为总的探测深度；
D N＝[D1,D2,…,DN-1]为第1至N层土的下边界深度组成的整体边界向量，D1,D2,…,DN-1表示第1,2,…,N-1层土的下边界深度，第N层土的下边界深度DN是确定的，其值等于总的探测深度H；ξ＝[ξ 1,ξ 2,...,ξ N]T表示第1至N层土内的CPT数据组成的整体CPT数据向量；
根据全概率理论和不同土层之间相互独立的假设，似然函数P(ξ|D N,N)可以写成：
式中：P(ξ n|D N,N)表示给定土层数目和边界深度条件下第n层土的CPT数据的概率密度函数；
P(ξ n|θ n,D N,N)为给定一套模型参数θ n条件下根据MP(θ n)建立的ξ n的联合概率密度函数；
P(θ n|D N,N)为模型参数的先验分布，根据岩土勘察报告和文献资料确定；
MP(·)表示步骤二中确定的概率模型；ξ n表示第n层土的CPT数据组成的向量；θ n表示第n层土的概率模型参数；n＝1,2…,N；
步骤四：基于贝叶斯理论确定在已知CPT数据条件下土层数目的后验分布的贝叶斯公式：
P(N|ξ)＝P(ξ|N)P(N)/P(ξ) 公式三
式中：N表示土层数目，N＝1,2,…,Nmax；Nmax表示可能的最大土层数目，根据岩土勘察报告和CPT数据进行预判；
P(N|ξ)表示土层数目后验分布；
P(ξ|N)为给定土层数目条件下CPT数据的概率密度函数，称作模型证据；
P(N)是土层数目的先验概率，在先验信息不充足的情况下，假设每个可能的N值具有相同的先验概率，即P(N)＝1/Nmax；
P(ξ)为归一化常数，与土层数目无关；
步骤五：确定给定土层数目条件下最可能土层边界深度D N,MPV，具体实现如下：
当N＝1时，最可能土层边界深度D N,MPV即为总的探测深度H；
当N≠1时，则需要采用模拟退火优化算法，模拟退火优化算法中优化的目标函数为公式一中土层边界深度后验分布的负对数即fobj＝-ln[P(D N|ξ,N)]，使目标函数fobj的值最小的土层边界深度为最可能的土层边界深度D N,MPV；
步骤六：获取每种可能土层数目的模型证据值，具体实现如下：
当N＝1时，模型证据值等于似然函数值；
当N≠1时，根据拉普拉斯逼近方法估计每种可能土层数目的模型证据值的计算公式：
P(ξ|N)≈P(ξ|D N,MPV,N)P(D N,MPV|N)(2π)(N-1)/2|detJ N(D N,MPV)|-1/2 公式四
式中：N表示土层数目，N＝2,…,Nmax，Nmax表示可能的最大土层数目；
D N,MPV表示土层数目为N时最可能的土层边界深度；
P(ξ|D N,MPV,N)表示土层数目为N、土层边界深度为D N,MPV条件下ξ的概率密度函数；
P(D N,MPV|N)表示土层数目为N条件下土层边界深度的先验分布；
detJ N(D N,MPV)表示土层数目为N条件下，目标函数fobj在最可能土层边界深度D N,MPV的二阶偏导的Hessian矩阵的行列式的值；
Hessian矩阵J N的第i行第j列的元素是i,j＝1,2,…,(N-1)；fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
若N≠Nmax，将土层数目N改为N+1重复步骤五、六；若N＝Nmax，进入下面步骤；
步骤七：确定最可能土层数目N*，具体实现如下：
根据贝叶斯模型选择方法，比较已计算出的Nmax个模型证据值，选出最大模型证据值对应的土层数目即为最可能土层数目N*；
步骤八：最可能土层数目N*条件下土层边界深度的不确定性表征，具体实现如下：
根据拉普拉斯逼近方法，最可能土层数目N*条件下土层边界深度后验分布的协方差矩阵计算公式：
上式中：Hessian矩阵的第i行第j列的元素是
为最可能土层数目N*条件下土层边界深度后验分布的协方差矩阵；
表示土层数目为N*时最可能的土层边界深度；
fobj＝-ln[P(D N|ξ,N)]；Di和Dj为向量D N中的第i，j个元素；
协方差矩阵的对角线即为土层边界深度的方差，进而计算出土层边界深度的标准差表征土层边界深度不确定性。
2.根据权利要求1所述的基于模拟退火的静力触探土层自动识别方法，其特征在于：
步骤一中获取的一套CPT数据为任意一套根据CPT测量的数据或由CPT测量数据计算得到的数据。
3.根据权利要求1所述的基于模拟退火的静力触探土层自动识别方法，其特征在于：步骤四中可能的最大土层数目Nmax，根据模型证据值随N的变化趋势判断所取可能的最大土层数目Nmax是否足够大，若模型证据值随N的变化无下降趋势，则Nmax取值不够大，此时Nmax改为Nmax+1，使N＝Nmax+1，重复步骤五、六、七；若模型证据值随N的变化有下降趋势，则说明Nmax取值足够大。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710707758.2A CN107330569A (zh) | 2017-08-17 | 2017-08-17 | 基于模拟退火的静力触探土层自动识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710707758.2A CN107330569A (zh) | 2017-08-17 | 2017-08-17 | 基于模拟退火的静力触探土层自动识别方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107330569A true CN107330569A (zh) | 2017-11-07 |
Family
ID=60224200
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710707758.2A Pending CN107330569A (zh) | 2017-08-17 | 2017-08-17 | 基于模拟退火的静力触探土层自动识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107330569A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109063285A (zh) * | 2018-07-18 | 2018-12-21 | 南昌大学 | 一种土坡钻孔布置方案设计方法 |
| CN111651879A (zh) * | 2020-05-28 | 2020-09-11 | 南京工业大学 | 一种基于压缩感知的地勘方案优化方法 |
| CN112396130A (zh) * | 2020-12-09 | 2021-02-23 | 中国能源建设集团江苏省电力设计院有限公司 | 静力触探试验岩层智能识别方法、系统、计算机设备及介质 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103741658A (zh) * | 2014-01-08 | 2014-04-23 | 江苏省水利科学研究院 | 采用探地雷达及静力触探仪联合勘测吹填砂量的方法 |
-
2017
- 2017-08-17 CN CN201710707758.2A patent/CN107330569A/zh active Pending
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103741658A (zh) * | 2014-01-08 | 2014-04-23 | 江苏省水利科学研究院 | 采用探地雷达及静力触探仪联合勘测吹填砂量的方法 |
Non-Patent Citations (2)
| Title |
|---|
| ZIJUN CAO ET AL.: "Probabilistic Site Characterization Using Cone Penetration Tests", 《SPRINGER》 * |
| 曹子君 等: "基于静力触探的土层自动划分方法与不确定性表征", 《岩土工程学报》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109063285A (zh) * | 2018-07-18 | 2018-12-21 | 南昌大学 | 一种土坡钻孔布置方案设计方法 |
| CN109063285B (zh) * | 2018-07-18 | 2022-12-02 | 南昌大学 | 一种土坡钻孔布置方案设计方法 |
| CN111651879A (zh) * | 2020-05-28 | 2020-09-11 | 南京工业大学 | 一种基于压缩感知的地勘方案优化方法 |
| CN111651879B (zh) * | 2020-05-28 | 2024-04-05 | 南京工业大学 | 一种基于压缩感知的地勘方案优化方法 |
| CN112396130A (zh) * | 2020-12-09 | 2021-02-23 | 中国能源建设集团江苏省电力设计院有限公司 | 静力触探试验岩层智能识别方法、系统、计算机设备及介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Manzocchi | The connectivity of two‐dimensional networks of spatially correlated fractures | |
| Scheidt et al. | Representing spatial uncertainty using distances and kernels | |
| CN110674841B (zh) | 一种基于聚类算法的测井曲线识别方法 | |
| CN105760673B (zh) | 一种河流相储层地震敏感参数模板分析方法 | |
| US10267934B2 (en) | System and method for generating a depositional sequence volume from seismic data | |
| CN107330569A (zh) | 基于模拟退火的静力触探土层自动识别方法 | |
| CN110909488A (zh) | 一种高效边坡可靠度分析方法 | |
| US20140297186A1 (en) | Rock Classification Based on Texture and Composition | |
| Olea | A practical primer on geostatistics | |
| CN104011566A (zh) | 用于分析地质构造的特性的基于小波变换的系统和方法 | |
| CN109800954A (zh) | 基于测井数据的储层评价新方法 | |
| Tang et al. | Multivariate statistical log log-facies classification on a shallow marine reservoir | |
| CN104992178A (zh) | 基于支持向量机模拟交会图的致密砂岩流体类型识别方法 | |
| Carvalho et al. | Soil classification system from cone penetration test data applying distance-based machine learning algorithms | |
| Naraghi et al. | 3D reconstruction of porous media from a 2D section and comparisons of transport and elastic properties | |
| Majdi et al. | Identification of well logs with significant impact on prediction of oil and gas reservoirs permeability using statistical analysis of RSE values | |
| Rodriguez et al. | New approach to identify analogue reservoirs | |
| Da Silva et al. | A new approach to soil classification mapping based on the spatial distribution of soil properties | |
| Yu et al. | Training image optimization method based on convolutional neural network and its application in discrete fracture network model selection | |
| CN111594156B (zh) | 一种天然气水合物饱和度计算方法及系统 | |
| CN111542819B (zh) | 用于改进的地下数据处理系统的装置和方法 | |
| CN114021644A (zh) | 一种基于K-means和去丛聚法的区域代表性地下水位计算方法 | |
| D'Windt et al. | Bayesian-Based approach for hydraulic flow unit identification and permeability prediction: a field case application in a tight carbonate reservoir | |
| Skjæveland et al. | Seismic Tiles, a data format to facilitate analytics on seismic reflectors | |
| Ismagilov et al. | Machine learning approach to open hole interpretation and static modelling applied to a giant field |
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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20171107 |
|
| RJ01 | Rejection of invention patent application after publication |