\chapter{知识模型设计} 数字孪生已经从理论研究发展到实际应用，其中模型是数字孪生的重要组成部分，也是成功数字孪生应用的先决条件。在虚拟空间中，数字孪生模型可以通过几何、物理、行为和规则四个模型维度来表达物理实体的属性\cite{ref17}。几何模型描述了物理实体的几何形状和装配关系。物理模型反映了物理实体的物理特性、特征和约束。行为模型表示物理实体对内部和外部机制的动态行为响应。规则模型结合历史数据，可以利用隐性知识，使数字孪生模型更加智能。通过融合多学科知识，多维数字孪生模型可以在数字世界中执行描述、诊断、预测、决策等功能。在本平台中，为了更方便统一的去构建模型，将行为与规则模型统一称为知识模型， 本章将讲述三种构建知识模型的方法。 \section{知识模型} 在建立数字孪生虚拟世界过程中，会有大大小小的规则、算法、逻辑等知识纳入其中，任何规则算法归纳来讲都是描述任意虚拟或实体对象相关作用关系，这些作用关系对于该系统内部世界来说就是一个个具体的知识。本文归纳了这些知识的共同特征，建立了一个一般性的知识模型，该模型是一个动态模型，其规则求解过程可根据特点纳入平台内运算或使用独立进程进行运算求解。 如图\ref{fig-f12}所示，基本模型需要先建立三个外部特征，触发域、输入域、输出域。触发域主要用于设定触发求解过程的条件。输入域是指定该模型需要哪些虚拟/实体对象的什么参数。输出域是指定该模型输入影响范围和参数。 对于不同的知识或者说规则，根据其作用域和功能特点，将其分为现实类知识、虚拟类知识、仿真类知识 \begin{figure}[h!] \centering \includegraphics[width=1\textwidth]{figure/f12.png} \caption{知识模型分类} \label{fig-f12} \end{figure} \begin{description} \item[现实类知识]该类知识是具体场景下最底层知识，其触发域和输入域仅能为实体对象数据，虚拟场景运行与否不影响该类知识触发和生效作用。 \item[仿真类知识]该类知识为虚拟环境中与实体对象对应的虚拟对象的数据的关联规则，其输出域可以设定于虚拟世界，也可作用于现实世界。 \item[虚拟类知识]该类知识其输出域仅能作用于虚拟世界。 \end{description} 如图\ref{fig-f9}所示，基本模型内部求解过程可以根据类型和功能特点选择提交公式或算法代码，或直接运行独立的求解程序。 本章后面以构建碰撞检测和刚体运动力学规则为例讲解代码和公式确定知识模型。 \begin{figure}[h!] \centering \includegraphics[width=0.8\textwidth]{figure/f9.png} \caption{知识模型上传过程} \label{fig-f9} \end{figure} \section{建模方法} 上节讲述了对于知识模型的一般性描述，本节将讲述本平台构建知识模型的三种建模方法，分别是公式建模、代码建模、程序建模。 \subsection{公式建模} 公式建模指直接在平台内定义好输入输出间公式关系和其触发条件，平台将在需要时自动根据其触发条件填入格式化输入数据并计算其结果，其流程如图\ref{fig-f9}所示，用户编写公式，并在测试后上传平台解析存储，平台在运行时检测触发条件触发后，将输入信息调入沙盒内，执行公式并计算其结果后编码回传平台。效果如图\ref{fig-f40}所示，其语法规则参考latex设计，对于常用函数做了部分精简，需要计算的部分用$?$表示，输入参数将以其名为变量出现。 \begin{figure}[h!] \centering \includegraphics[width=1\textwidth]{figure/f40.png} \caption{公式建模编辑测试界面} \label{fig-f40} \end{figure} \subsection{代码建模} 代码建模指用户自己编写一段函数代码，定义好输入输出，并在测试后提交到平台内，其执行过程同公式模型一致，如图\ref{fig-f9}流程所示，代码在用户上传完毕后会根据需要在平台上生成运行环境，并根据需要提供一定的运算资源和数据背景，提供更高阶的数据查询和接入方法。目前代码解析部分采取现成的语法，根据语言特性和场景需要，目前仅支持js和python，尽管采用了一定的沙盒机制去运行，但是仍然面临很大的安全问题，所以现在的代码机制仅供内部使用，后续根据需要开放出自定义代码或者代码检查。 \begin{figure}[h!] \centering \includegraphics[width=1\textwidth]{figure/f42.png} \caption{代码建模编辑界面} \label{fig-f42} \end{figure} \begin{figure}[h!] \centering \includegraphics[width=1\textwidth]{figure/f43.png} \caption{代码测试运行界面} \label{fig-f43} \end{figure} 如图\ref{fig-f42}所示，为代码编辑界面，用户需要创建main函数作为入口程序，获取参数通过$P$函数获取，其参数为运行场景内参数的id号，平台提供代码和参数提示功能，在此可直接输入相关中文参数路径，根据提示选择后会自动生成id号，右上角提供运行功能，如图\ref{fig-f43}其运行测试功能目前仅处于测试阶段，仅能在控制台查看。 \subsection{程序建模} 程序建模指用户自己编写代码并在自己指定设备运行，但是需要调用平台提供的api依赖库，这种方式自由度和保密性最高，但是不易分享，上述两种模型编写好后可以选择公开。如图\ref{fig-f9}所示，用户仅需要申请访问秘钥，调用本平台开发的依赖库即刻纳入到通信网络中作为一个计算节点执行计算任务，与内部的计算节点无任何差别，唯一要考虑的是跨地区网络后消息实时性无法满足要求，不适合去做实时性分析人物。目前支持的依赖库python/go/c/rust。 \subsection{解析算法} 在构建知识模型其中很重要的一步是支持公式和代码输入，对于一条规则，我们可以提供节点数据并转换成合适的格式，比如对于某一理想形状体，我给予100N的力10s，如何通过规则模型解析公式获得该物体10s内的运动状态，首先我们需要去解析公式，确定输入变量数量和格式，在确定输出变量数量和格式，根据公式以一定的频率去输出结果。 在公式撰写规则上直接采用改造的latex语法去语义化编辑公式，采用抽象语法树（AST）算法去解析公式结构\cite{ref37}。 一般来说，语法分为具体语法和抽象语法。第一种是语法的文本形式，它是人类可读和可输入的，这种语法表示是使用字符串(数组或字符列表)实现的，其优点是它可以很容易地被人类阅读，这种语法涉及一个基于字符串的简单计算模型。然而，这种表现方式的缺点是众多而严重的。具体语法包含太多对许多操作不重要的信息，如空白、中缀/前缀表示法和关键字，重要的计算信息，如递归结构、函数-参数关系、项-子项关系等，没有被显式表示\cite{ref38}。具体语法的计算成本可以通过将具体语法解析为抽象语法树(AST)来克服，其实现使用的是标记树或链表，并使用构造函数和析构函数(例如Lisp中的car/cdr/cons)来处理。AST可以用来更好的表示递归结构，并且可以省略很多的用于方便人类编写具体语法而产生的结构，如缩进、间距、关键词等，其解析算法也更易于计算机处理和内存存储。 在对代码和公式进行解析前，我们需要做一些预处理，比如去除空格，去除注释，去除换行符，去除多余的括号等，其中最重要的是对于公式表达式的转换。公式表达式分为前缀、中缀以及后缀表示法，其中中缀表示法即为人类阅读的格式，前缀以及后缀表达式也成为波兰表示法和逆波兰表示法，通过Kasprzyk\cite{ref39}等人的工作已经这些年语法处理技术的发展，逆波兰表示法成为更适于计算机处理的公式表示法。 因此对于公式的处理过程是先进行预处理，将公式表达式转换为逆波兰表示法，然后通过逆波兰表示法构建抽象语法树，最后通过抽象语法树去解析公式，得到输入变量和输出变量的数量和格式，然后通过公式去计算结果。 将中缀表达式转换为逆波兰表达式（RPN）的算法\cite{ref40}如下: \begin{enumerate} \item 检查中缀表达式的当前元素 \item 如果当前元素是操作数，则将其发送到输出并转到步骤6 \item 如果当前元素是左括号，请按该元素并转到步骤6 \item 如果当前元素是运算符，则执行以下操作：如果它的优先级高于堆栈顶部推送该运算符，否则从堆栈中弹出运算符，将其发送到输出并重复步骤4 \item 如果当前元素是右括号，则从堆栈中弹出运算符并将它们发送到输出，直到弹出左括号 \item 如果输入表达式有更多元素，请转到步骤2，否则，弹出堆栈的其余部分，将其发送到输出并停止。 \end{enumerate} 以该公式为例: \begin{equation} y = x1 / 3 + (5 + sin(x1)) * x2 \label{eq3} \end{equation} 将右侧表达式转换为RPN表达式 \begin{equation} x1\ 3\ /\ 5\ x1\ sin\ +\ x2\ *\ + \label{eq4} \end{equation} 其RPN表达式存储结构为 \begin{lstlisting}[ language={python}, caption={RPN表达式存储结构}, label={code-c1}, ] { "type": "Program", "body": [ { "type": "ExpressionStatement", "expression": { "type": "BinaryExpression", "operator": "+", "left": { "type": "BinaryExpression", "operator": "/", "left": { "type": "Identifier", "name": "x1" }, "right": { "type": "Literal", "value": 3, "raw": "3" } }, "right": { "type": "BinaryExpression", "operator": "*", "left": { "type": "BinaryExpression", "operator": "+", "left": { "type": "Literal", "value": 5, "raw": "5" }, "right": { "type": "CallExpression", "callee": { "type": "Identifier", "name": "sin" }, "arguments": [ { "type": "Identifier", "name": "x1" } ] } }, "right": { "type": "Identifier", "name": "x2" } } } } ], "sourceType": "script" } \end{lstlisting} 在得到公式的RPN表达式后，就是构建抽象语法树，上步构建的RPN表达式\ref{eq4}时存储在栈中，将元素出栈，并根据操作符右序构建即可，其操作步骤如下： \begin{table} \centering \caption{RPN表达式\ref{eq4}构建AST\ref{tik1}步骤} \label{tab-c6} \begin{tabular}{ccl} \toprule Token & 类型 & 步骤 \\ \midrule + & 二元操作符（BinaryExpression） & 出栈，栈顶为root，二元操作符，产生左右两子节点 \\ * & 二元操作符 & 出栈，root的右子节点，二元操作符，产生两个子节点 \\ x2 & 标识符（Identifier） & 出栈，*的右节点，标识符，表示为变量 \\ + & 二元操作符 & 出栈，*的左节点，二元操作符，产生两个子节点\\ sin & 调用操作符（CallExpression） & 出栈，+ 的右节点，调用操作符，产生一个子节点\\ x1 & 标识符 & 出栈，sin的唯一子节点，标识符\\ 5 & 数字 & 出栈，+的左节点，数值，无节点\\ / & 二元操作符 & 出栈， 栈顶root的左节点， 二元操作符， 产生两个子节点\\ 3 & 数字 & 出栈，/ 的右节点， 数值，五子节点 \\ x1 & 标识符& 出栈， / 的左节点， 标识符， 无子节点 \\ \bottomrule \end{tabular} \end{table} 根据RPN表达式\ref{eq4}构建出的AST语法树为 \begin{center} \begin{tikzpicture} [thick,scale=1, every node/.style={scale=1}] \node {+} child {node {/} child {node {x1}} child [missing] {} child {node {3}} } child [missing] {} child [missing] {} child [missing] {} child { node {*} child {node {+} child {node {5}} child {node {sin} child [missing] {} child {node {x1}} child [missing] {} } } child [missing] {} child {node {x2}} }; \end{tikzpicture} \label{tik1} \end{center} \begin{figure}[h!] \centering \includegraphics[width=1\textwidth]{figure/f41.png} \caption{公式及代码解析构建AST树} \label{fig-f41} \end{figure} 如图\ref{fig-f40}图\ref{fig-f41}所示，以AST树构建公式算法，并根据节点内容渲染不同类型公式以及计算，输入目前支持浮点数及矩阵格式，支持常用函数表达，需要注意的是设计节点间响应的具体公式需要去做输入输出数据格式转换。 \section{代码模型案例-碰撞检测} 碰撞检测目前定义为虚拟类规则，暂时其输出域数据无法作用于现实世界，主要用于虚拟世界用户交互检测。 为了简化物体之间的碰撞检测运算，通常会对物体创建一个规则的几何外形将其包围。在本系统中，碰撞检测中将物体分为三种检测模型，点、AABB、球体。其中，AABB（axis-aligned bounding box）包围盒被称为轴对齐包围盒。 轴对齐包围盒是判断两个物体是否重叠的最快算法，物体被包裹在一个非旋转的（因此轴对齐的）盒中，并检查这些盒在三维坐标空间中的位置，以确定它们是否重叠。 由于性能原因，轴对齐是有一些约束的。两个非旋转的盒子之间是否重叠可以通过逻辑比较进行检查，而旋转的盒子则需要三角运算，这会导致性能下降。如果你有旋转的物体，可以通过修改边框的尺寸，这样盒子仍可以包裹物体，或者选择使用另一种边界几何类型，比如球体 (球体旋转，形状不会变)。 在该例中，根据不同需要适合使用代码片段去定义该规则。 \subsection{点与AABB} 如果检测到一个点是否在 AABB 内部就非常简单了 — 我们只需要检查这个点的坐标是否在 AABB 内; 分别考虑到每种坐标轴。如果假设 $P_x$, $P_y$ 和 $P_z$ 是点的坐标, $B_{minX} - B_{maxX}$, $B_{minY} - B_{maxY}$, 和$B_{minZ}–B_{maxZ}$ 是 AABB 的每一个坐标轴的范围，我们可以使用以下代码计算两者之间的碰撞是否发生： \begin{lstlisting}[ language={}, label={code-js-sample-1}, ] function isPointInsideAABB(point, box) { return (point.x >= box.minX && point.x <= box.maxX) && (point.y >= box.minY && point.y <= box.maxY) && (point.z >= box.minY && point.z <= box.maxZ); } \end{lstlisting} \subsection{AABB与AABB} 检查一个 AABB 是否和另一个 AABB 相交类似于检测两个点一样。我们只需要基于每一条坐标轴并利用盒子的边缘去检测。 \begin{lstlisting}[ language={}, label={code-js-sample}, ] function intersect(a, b) { return (a.minX <= b.maxX && a.maxX >= b.minX) && (a.minY <= b.maxY && a.maxY >= b.minY) && (a.minZ <= b.maxZ && a.maxZ >= b.minZ); } \end{lstlisting} \section{公式模型案例-刚体运动力学模型} 在虚拟世界中为关联现实设备运动状态，需要根据已有的传感器数据，如加速度、里程信息等估计实体对象运动姿态，或者由虚拟对象指导影响实体对象运动，两者相互作用皆需要实现基本的刚体运动。刚体运动的典型计算特征可用于验证公式规则模型，其评估效率和验证性能直观有效。 刚体的运动主要基于牛顿三大定律来模拟: \begin{itemize} \item 1.惯性 物体在不受力时，总是保持速度不变. \item 2.力，质量，加速度力在物体上产生加速度，满足$F = ma$. \item 3.力的作用是相互的 \end{itemize} 基于牛顿三大定律，对于每个物体视作刚体运动，如图\ref{fig-f44}使用迭代逼近的方式来求解稳态位姿: \begin{figure}[h!] \centering \includegraphics[height=4cm]{figure/f44.png} \caption{刚体稳态求解流程图} \label{fig-f44} \end{figure} 最后将所有运动分解为粒子运动和旋转运动，每秒以一定频率分别计算其空间位置和空间姿态，最终合成为完整运动分析过程。 \subsection{空间位置} 在做初步解算时，先不考虑物体的形状和旋转，把物体当成理想粒子来对待，根据牛顿定律来循环计算物体的速度和位置: \[ dt = t_{i+1} - t_{i} \] \[ v(t_{i+1}) = v(t_{i}) + (\frac{f(t_{i})}{m})dt \] \[ p(t_{i+1}) = v(t_{i}) + v(t_{i+1})dt \] 将其泰勒展开 \[ p(t_{i+1}) = p(t_{i}) + p^{'}(t_{i})dt + p^{''}(t_{i})\frac{dt^{2}}{2!} + p^{'''}(t_{i})\frac{dt^{3}}{3!} + ... \] 对于浏览器按像素点实时位置渲染精度而言，3阶泰勒展开已符合精度要求。 因此对于理想点的运动规则模拟可以提交公式 \begin{lstlisting}[ language={}, label={code}, ] p(t_{i+1}) = p(t_{i}) + p^{'}(t_{i})dt + p^{''}(t_{i})\frac{dt^{2}}{2!} + p^{'''}(t_{i})\frac{dt^{3}}{3!} \end{lstlisting} \subsection{空间姿态} 要计算物体的空间姿态，即物体绕某一定点在三自由度上的旋转量，可以用四元数q来表示刚体在极短时间内的旋转量。 在三维空间中任取一点$Q$ 及一个直角坐标系$Q_{xyz}$，可以得到物体的惯性张量: \begin{equation} I ={ \left[ \begin{array}{ccc} I_{xx} & I_{xy} & I_{xz}\\ I_{yz} & I_{yy} & I_{yz}\\ I_{zx} & I_{zy} & I_{zz} \end{array} \right ]} \end{equation} 对角元素$I_{xx}$,$I_{yy}$,$I_{yy}$是物体分别相对于$x$，$y$，$z$轴的转动惯量。 \[ I_{xx} = \int{(r_{y}^{2}+r_z^2)dm} \] \[ I_{yy} = \int{(r_{x}^{2}+r_z^2)dm} \] \[ I_{zz} = \int{(r_{x}^{2}+r_y^2)dm} \] 计算惯量积 \[ I_{xy} = I_{yx} = - \int{(r_x r_y)dm} \] \[ I_{xz} = I_{zx} = - \int{(r_x r_z)dm} \] \[ I_{yz} = I_{zy} = - \int{(r_z r_y)dm} \] 根据惯性张量计算力矩$\tau$，角将速度$d\omega$，角速度$\omega$。 \[ \tau = \frac{dL}{dt} = I \frac{d\omega}{dt} \] \[ d\omega = I^{-1} \tau dt \] \[ \omega(t_{i+1}) = \omega(t_i) + (I^{-1} \tau(t_i)) dt \] 用四元数$q$来表示刚体的旋转状态，刚体在$dt$时间内沿着$\omega$旋转了$|\omega dt|$的角度，得到 \begin{equation} q(t_{i+1}) = q(t_i) * [\cos{\frac{|\omega dt|}{2}},\sin{\frac{|\omega dt|}{2}}\frac{\omega}{|\omega|}] \end{equation} 对于获取刚体旋转状态的规则特征可以提交为 \begin{lstlisting}[ language={}, label={code2}, ] q(t_{i+1}) = q(t_i) * [\cos{\frac{|\omega dt|}{2}},\sin{\frac{|\omega dt|}{2}}\frac{\omega}{|\omega|}] \end{lstlisting}