电阻抗成像（electrical impedance tomography，EIT）是功能性成像技术，其原理是通过对附在目标表面的电极施加安全的激励电流，测量由此电流产生的目标表面的电压信号，重构出被测目标内部的阻抗图像分布^[1]. 由于其具有成像迅速、无损伤、价格低廉等特点，在水文监测和医学诊断方面具有广泛的应用前景^[2-3]. 在生物EIT研究中，为了便于分析，大多认为生物组织的电导率是确定性的量. 但是Kuwahara等^[4-8]的研究表明，生物组织的电导率不仅会随着组织本身的变化（如组织癌变、肿瘤或脓肿等病变）而变化，还会受外部某些物理环境（如成像频率、电刺激或温度等）的影响. 因此，生物组织的电导率不是固定不变的量，而是变化的、不确定性的量，研究电导率分布的不确定性对EIT的影响具有重要意义.

不确定性量化（uncertainty quantification，UQ）在实际问题的优化和决策过程中，对降低不确定性参数的影响起着关键作用^[9]. UQ主要研究的是模型中数值和物理参数的变化如何影响模拟的结果，其中“量化”指将模型中所有不确定性信息最终以数学或计算的形式表达出来的能力. UQ方法大体分为2类：概率框架下和非概率框架下^[10]. 常用的概率UQ方法主要有蒙特卡洛（Monte Carlo, MC）法、混沌多项式展开（polynomial chaos expansion, PCE）法和稀疏网格配点法等. MC法是经典的UQ方法，具有原理简单、适用性强的特点，但在运算中存在计算量大、收敛速度慢的问题，从而发展出其他改进抽样方法的MC法，如拉丁超立方抽样^[11]. PCE法由Ghanem等^[12]提出，并应用于动力学不确定性问题中，其具有指数收敛速度，因此逐渐替代MC法. Xiu等^[13]结合MC法和随机Galerkin法的优点，提出一类高阶配点法，其中应用最为广泛的是稀疏网格配点法，该方法具有较高的精度和收敛速度，且更适用于处理高维不确定性问题. Nobile等^[14]对稀疏网格配点法进行具体分析，通过实际算例验证稀疏网格配点法处理高维问题的有效性.

在EIT的研究中，有关UQ方法的应用，Hyvönen等^[15]采用随机Galerkin有限元法对EIT正问题进行数值求解. Sun等^[16]采用混沌多项式展开法构造出EIT电导率场的替代模型，研究不同测量误差对EIT图像重构的影响. 有关EIT中不确定性参数的研究，Kolehmainen等^[17]研究EIT模型边界形状的不确定性，利用Teichmuller映射理论从EIT的测量数据中恢复不确定性信息. Boyle等^[18]研究电极运动的不确定性对EIT的影响，通过计算电极位置的雅可比矩阵重构电极运动. Nissinen等^[19]研究EIT电极接触阻抗的不确定性，通过将其建模为噪声过程，使用贝叶斯建模误差方法补偿未知的电极接触阻抗引起的误差. 本研究针对EIT正问题中电导率分布的不确定性对输出电位的影响，基于灵敏度分析提出新的UQ方法：改进的稀疏网格配点法. 以4层同心圆头模型为算例，将改进方法的计算结果与MC法、PCE法和稀疏网格配点法相比，验证所提改进方法的高效性.

MC法是典型的随机抽样法，根据随机输入变量的概率分布选取大量样本，通过反复迭代计算使估算结果收敛至原模型函数的真实值. MC法原理简单，当选取的样本点足够多时，其能够较准确地还原真实原模型，因此常用MC法作为基准检验其他不确定性量化方法的有效性和精度^[20].

MC法的关键是选取样本点，在随机输入变量 ${\boldsymbol{X}} = [{X_1},{X_2},\cdots,{X_d}]$的整个变化范围内随机选取N个样本点 ${{\boldsymbol{x}}_j} = [{x_{j1}},{x_{j2}},\cdots,{x_{jd}}]$， $j = 1,2,\cdots,N$；逐次将样本点 ${{\boldsymbol{x}}_j}$代入原模型函数 $f(X)$中，计算模型输出值 ${f_j}$；根据所有输出值求解模型函数的统计信息（如均值和标准差）.

PCE法是随机展开法，它用多组有限阶截断展开的正交多项式表示模型输出响应，根据多项式的正交性质，直接从展开系数中计算模型输出的统计特性^[21]. 针对PCE模型，其展开表达式为

(1) $ \begin{split} &f = {\alpha _0}{H_0} + \sum\limits_{{j_1} = 1}^p {{\alpha _{{j_1}}}{H_1}({\xi _{{j_1}}})} + \\ & \quad \sum\limits_{{j_1} = 1}^p {\sum\limits_{{j_2} = 1}^{{j_1}} {{\alpha _{{j_1}{j_2}}}{H_2}({\xi _{{j_1}}},{\xi _{{j_2}}})} } {\text{ + }} \cdots {\text{ + }} \\ & \quad \sum\limits_{{j_1} = 1}^p {\sum\limits_{{j_2} = 1}^{{j_1}} { \cdots \sum\limits_{{j_d} = 1}^{{j_{d - 1}}} {{\alpha _{{j_1}{j_2}\cdots{j_d}}}{H_p}({\xi _{{j_1}}},{\xi _{{j_2}}},\cdots,{\xi _{{j_d}}})} } }. \end{split} $

式中： $ {\alpha _{{j_1}{j_2}\cdots{j_d}}} $为多项式展开系数，p为多项式截断阶数， ${H_p}(\xi )$为p阶的埃尔米特正交多项式， $\xi = ({\xi _{{j_1}}},{\xi _{{j_2}}},\cdots,{\xi _{{j_d}}})$为标准随机变量.

PCE法将原模型输出响应表示为p阶的PCE模型，再在标准随机变量空间中选取样本点；通过线性变换将样本点转换至原模型空间并计算各样本点在原模型上的函数值；通过估算多项式展开系数，求解模型输出的统计信息（如均值和标准差）.

稀疏网格配点法源于smolyak算法^[22]，利用多项式函数逼近原随机模型，通过在样本点上进行线性回归求得多项式函数的系数，构造出原随机模型的替代模型，并在替代模型的基础上进行不确定性量化分析. 其基本原理是：从一维形式出发，利用特殊的张量积运算将一维形式的样本点和基函数扩展到多维情况^[23].

含有不确定性参数 ${\boldsymbol{\lambda}} \in {{\text{R}}^d}$的原随机模型 $f({\boldsymbol{\lambda}} )$，通过构造替代模型将原随机模型的隐式表达式转化为显式表达式，可简写为

(2) $ f({\boldsymbol{\lambda}} ) \approx \hat f({\boldsymbol{\lambda}} ) = \sum\limits_{n = 1}^N {{{\boldsymbol{b}}_n} {\varPhi _n}} ({\boldsymbol{\lambda}} ). $

式中： ${{\varPhi} _n}({\boldsymbol{\lambda}} ){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}$为d维基函数，b_n（n₁=1,···N）为系数向量.

稀疏网格配点法定义多指数组合i=[i₁,···,i_d]，其中i_j (j=1,···,d)决定了第j维上所配置的样本点和基函数的数目，其取值满足：

(3) $ d \leqslant \left| i \right| \leqslant k + d{\kern 1pt} . $

式中：|i|=i₁+···+i_d，其中i₁,···,i_d为对应于各维变量的指数，称为多指数；d为原随机模型的维数；k为各维变量的精度水平. 在稀疏网格配点法中，每个维度(1,2,···,d)被分配相同的精度水平k，即K=[k, k,···,k]. 通常，所选精度水平越高所得计算结果的精度越高.

通过稀疏网格配点法构造的替代模型的显式表达式为

(4) $ \begin{gathered} \hat f({\boldsymbol{\lambda}} ) = \sum\limits_{d \leqslant \left| i \right| \leqslant k + d} {{{( - 1)}^{k + d - \left| i \right|}}\left[ \begin{gathered} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} d - 1 \hfill \\ k + d - \left| i \right| \hfill \\ \end{gathered} \right]} \times \hfill \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\boldsymbol{U}}^{{i_1}}}({\lambda _1}) \otimes \cdots \otimes {{\boldsymbol{U}}^{{i_d}}}({\lambda _d}). \hfill \\ \end{gathered} $

式中： ${{\boldsymbol{U}}^{{i_j}}}({\lambda _j}) = \sum\limits_{\ell {\text{ = }}1}^{m({i_j})} {{{\boldsymbol{b}}_\ell }{\boldsymbol{\varPhi}} _\ell ^{{i_j}}({\lambda _j}){\kern 1pt} {\kern 1pt} {\kern 1pt} } {\kern 1pt} (j = 1,\cdots,d)$为一维多项式函数， ${\boldsymbol{\varPhi}} _\ell ^{{i_j}}({\lambda _j})$为一维基函数， ${{\boldsymbol{b}}_\ell }$为系数， $ \otimes $为张量积运算符.

灵敏度分析法是定量地研究各输入不确定性参数对模型输出响应的影响程度的方法^[24]. 全局灵敏度分析法的应用广泛，其考虑了整个模型内所有输入参数对模型输出响应的影响，以及不同输入参数间的耦合影响. 本研究采用基于方差的全局灵敏度分析法，该方法采用方差来衡量输入参数的不确定性对输出响应的影响，将输出响应的方差分解为各输入参数以及各输入参数间的耦合效应的共同影响^[25].

将含有d个独立输入参数的原随机模型 $f({\boldsymbol{\lambda}} )$展开为多个子项之和的形式，展开式为

(5) $ \begin{split} f({\boldsymbol{\lambda}} ) =& {f_0}{\text{ + }}\sum\limits_{p = 1}^d {{f_p}({\lambda _p})} + \sum\limits_{1 \leqslant p \leqslant q \leqslant d} {{f_{pq}}({\lambda _p},{\lambda _q})} + \cdots +\\ & {f_{12 \cdots d}}({\lambda _1},{\lambda _2},\cdots,{\lambda _d}). \hfill \\ \end{split} $

若函数 $f({\boldsymbol{\lambda}} )$平方可积且展开式中每项对其包含的任意参数的积分值均为零，即：

(6) $ \int_0^1 {{f_{{p_1},\cdots,{p_z}}}({\lambda _{{p_1}}},{\lambda _{{p_2}}},\cdots,{\lambda _{{p_z}}}){\rm{d}}{\lambda _{{p_t}}}} = 0;{\kern 1pt} {\kern 1pt} {\kern 1pt} 1 \leqslant t \leqslant z. $

则表明式（5）中各项均相互正交，并可唯一求得，即 ${f_0} = E(f),$ ${f_p} = E(f|{\lambda _p}) - {f_0} $, ${f_{pq}} = E(f|{\lambda _p},{\lambda _q}) - {f_0} - {f_p} - {f_q},$以此类推.

对式（5）两侧同时取方差，可得：

(7) $ {\text{Var}}(f) = \sum\limits_{p = 1}^d {{V_p}} + \sum\limits_{1 \leqslant p \leqslant q \leqslant d} {{V_{pq}}} + \cdots + {V_{12 \cdots d}}. $

式中： ${V_p} = {\text{Var}}({f_p}) = {\text{Var}}(E(f|{\lambda _p}))$， ${V_{pq}} = {\text{Var}}({f_{pq}}) = $ $ {\text{Var}} (E(f|{\lambda _p},{\lambda _q})) - {V_p} - {V_q}$，以此类推. $ E(·) $、 $ \text{Var}(·) $分别为期望、方差算子.

输入参数的灵敏度指标可表示为

(8) $ {S}_{p} = \frac{{{V_p}}}{{{\rm{Var}}(f)}},{S_{pq}} = \frac{{{V_{pq}}}}{{{\rm{Var}}(f)}},\cdots,{\kern 1pt} {\kern 1pt} {\kern 1pt} {S_{12 \cdots d}}{\text{ = }}\frac{{{V_{12 \cdots d}}}}{{{\rm{Var}}(f)}}. $

式中： ${S}_{p}$为参数 ${\lambda _p}$对输出响应的影响程度， ${S_{pq}}$为 ${\lambda _p}$、 ${\lambda _q}$的耦合效应对输出响应的影响程度. 各灵敏度指标的和满足：

(9) $ \sum\limits_{p = 1}^d {{S_p}} + \sum\limits_{1 \leqslant p \leqslant q \leqslant d} {{S_{pq}}} + \cdots + {S_{12 \cdots d}} = 1. $

灵敏度指标的值越大，说明该输入参数的变化对输出响应的影响程度越大.

稀疏网格配点法在处理不确定性问题时，采用特殊的选点规则，其样本点数目相比于传统张量积法有所减少，但存在一定的局限. 稀疏网格配点法没有考虑模型中各维变量的重要性程度，认为每个维度同等重要，并给各维变量配置相同的精度水平k（即相同数量的样本点），这会导致计算成本的浪费. 通常，各维变量对模型输出响应的影响程度是不同的，因此有重要维度和非重要维度之分，若能够减少某些非重要维度上的样本点，给重要的维度上分配较多样本点，使样本点可以更合理地配置，将有效节省计算成本.

本研究提出改进的稀疏网格配点法，该方法结合灵敏度分析允许对模型中各维变量进行不同的处理，可以根据灵敏度分析结果中每个维度不同的灵敏度指标配置不同的精度水平，即灵敏度指标高的维度上配置高的精度水平，将更多的计算资源分配到更重要维度上，从而提高计算精度和效率. 改进稀疏网格配点法的基本原理与稀疏网格配点法类似，不同的是改进方法考虑各维变量的重要性程度，并将不同的精度水平k₁,···,k_d分配给每个维度，即K=[k₁,···,k_d]. 改进方法在式（3）的基础上提出新的规则，以确定满足条件的多指数组合i=[i₁,···,i_d]：

(10) $ d \leqslant \left| i \right| \leqslant d + {k^{\max }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {i_j} \leqslant {k_j} + 1{\kern 1pt} {\kern 1pt} .{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} $

式中：k^max=(k₁,···,k_d)，其中k₁,···,k_d为对应于各维(1,···,d)的精度水平. 通常，精度水平选取得越高，维度上的样本点数就越多，所得计算结果也越准确. 改进方法需要求得各维变量对模型输出响应的影响程度. 通过基于方差的全局灵敏度分析法获得的结果，按各灵敏度指标的大小排序选择每个维度的精度水平k₁,···,k_d，并基于此求出满足式（10）的所有可能的多指数组合i=[i₁,···,i_d]，计算出相应的一维样本点和基函数. 本研究采用切比雪夫多项式作为一维基函数，切比雪夫多项式的极值作为一维样本点. 其中，各维上的一维样本点和基函数的数目m(i_j)和对应的多指数i_j (j=1,···,d)的关系为

(11) $ m({i_j}) = \left\{ \begin{gathered} {\text{1, }}\qquad\quad\,{i_j} = 1; \hfill \\ {2^{{i_j} - 1}} + 1, \;\;{i_j} > 1. \hfill \\ \end{gathered} \right. $

基于式（11）求得的一维样本点均分布在[−1,1]，因此需将其投影到原模型参数的分布区间内，则所得一维样本点为 ${\theta ^{{i_1}}},{\theta ^{{i_2}}}, \cdots ,{\theta ^{{i_d}}}$，对这些一维样本点进行特殊的张量积运算，得到d维样本点集 ${{\boldsymbol{\theta}} ^d}{\text{ = }}\left\{ {{\theta _1},{\theta _2},\cdots,{{{\theta}} _N}} \right\}$：

(12) $ {{\boldsymbol{\theta}} ^d} = \bigcup\limits_{d \leqslant \left| i \right| \leqslant d + {k^{\max }},\left\{ {{i_j} \leqslant {k_j}{\text{ + 1}}} \right\}_{j = 1}^d} {{{\boldsymbol{\theta}} ^{{i_1}}} \otimes {{\boldsymbol{\theta}} ^{{i_2}}} \otimes \cdots \otimes {{\boldsymbol{\theta}} ^{{i_d}}}} {\kern 1pt} . $

基于求得的一维基函数 ${{\boldsymbol{\varPhi}} ^{{i_1}}}({\lambda _1}),{{\boldsymbol{\varPhi}} ^{{i_2}}}({\lambda _2}), \cdots , {{\boldsymbol{\varPhi}} ^{{i_d}}} ({\lambda _d})$构造一维多项式函数 ${{\boldsymbol{U}}^i}({\boldsymbol{\lambda}} )$：

(13) $ {{\boldsymbol{U}}^i}({\boldsymbol{\lambda}} ) = \sum\limits_{\ell = m(i - 1){\text{ + }}1}^{m(i)} {{{\boldsymbol{b}}_\ell }} {\boldsymbol{\varPhi}} _\ell ^i({\boldsymbol{\lambda}} ). $

在各一维多项式函数的基础上进行张量积运算，得到其中一组多指数对应的d维多项式函数形式：

(14) $ \begin{split} & {{\boldsymbol{U}}^{{i_1}}}({\lambda _1}) \otimes \cdots \otimes {{\boldsymbol{U}}^{{i_d}}}({\lambda _d}) = \sum\limits_{{\ell _1} = m({i_1} - 1) + 1}^{m({i_1})} \cdots \\ & \sum\limits_{{\ell _d} = m({i_d} - 1) + 1}^{m({i_d})} {{{\boldsymbol{b}}_{{\ell _1}\cdots{\ell _d}}} \cdot {\boldsymbol{\varPhi}} _{{\ell _1}}^{{i_1}}({\lambda _1}) \otimes \cdots \otimes {\boldsymbol{\varPhi}} _{{\ell _d}}^{{i_d}}({\lambda _d})} . \end{split} $

式中： ${\boldsymbol{\varPhi}} _{{\ell _1}}^{{i_1}}({\lambda _1}) \otimes \cdots \otimes {\boldsymbol{\varPhi}} _{{\ell _d}}^{{i_d}}({\lambda _d})$为一组d维基函数， ${{\boldsymbol{b}}_{{\ell _1}\cdots{\ell _d}}}$为系数向量.

将所有满足要求的多指数组合对应的d维多项式函数形式，基于式（10）进行特殊的加权求和，由改进稀疏网格配点法构造的替代模型的显式表达式为

(15) $ \hat f({\boldsymbol{\lambda}} ) = \sum\limits_{d \leqslant \left| i \right| \leqslant d + {k^{\max }},\left\{ {{i_j} \leqslant {k_j}{\text{ + 1}}} \right\}_{j = 1}^d} {{{\boldsymbol{U}}^{{i_1}}}({\lambda _1}) \otimes \cdots \otimes {{\boldsymbol{U}}^{{i_d}}}} ({\lambda _d}). $

将d维样本点 ${\theta _1},{\theta _2},\cdots,{\theta _N}$输入原模型 $f(\lambda )$中，得到各样本点处的原模型输出值 ${\boldsymbol{F}} = [ f({\theta _1}),f({\theta _2}),\cdots, $ $ {f({\theta _N})} ]^{\rm{T}}$，将N个d维样本点和相应的原模型输出值代入式（2），得

(16) $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varPhi}}_1}\left( {{\theta _1}} \right)}&{{{\boldsymbol{\varPhi}}_2}\left( {{\theta _1}} \right)}& \cdots &{{{\boldsymbol{\varPhi}}_N}\left( {{\theta _1}} \right)}\\ {{{\boldsymbol{\varPhi}}_1}\left( {{\theta _2}} \right)}&{{{\boldsymbol{\varPhi}}_2}\left( {{\theta _2}} \right)}& \cdots &{{{\boldsymbol{\varPhi}}_N}\left( {{\theta _2}} \right)}\\ \vdots & \vdots &{}& \vdots \\ {{{\boldsymbol{\varPhi}}_1}\left( {{\theta _N}} \right)}&{{{\boldsymbol{\varPhi}}_2}\left( {{\theta _N}} \right)}& \cdots &{{{\boldsymbol{\varPhi}}_N}\left( {{\theta _N}} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{b_1}}\\ {{b_2}}\\ \vdots \\ {{b_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {f\left( {{\theta _1}} \right)}\\ {f\left( {{\theta _2}} \right)}\\ \vdots \\ {f\left( {{\theta _N}} \right)} \end{array}} \right]. $

可简写为

(17) $ {\boldsymbol{Ab}} = {\boldsymbol{F}}. $

式中：A为基函数矩阵. 利用最小二乘回归法求得系数向量b为

(18) $ {\boldsymbol{b}} = {({{\boldsymbol{A}}^{\rm{T}}}{\boldsymbol{A}})^{ - 1}}{{\boldsymbol{A}}^{\rm{T}}}{\boldsymbol{F}}. $

EIT正问题是利用给定的边界激励信号和已知的电阻抗分布求解EIT模型，得出求解区域内部及边界的电位分布. 考虑边界条件Γ，EIT求解场域Ω的数学模型可以描述为如下的拉普拉斯方程：

(19) $ \left. \begin{array}{l} \varOmega :\nabla \cdot \sigma \left( {\nabla {\varphi _0}} \right) = 0,\\ {\varGamma _1}:{\varphi _0} = \varphi ,\\ {\varGamma _2}:\sigma \frac{{\partial {\varphi _0}}}{{\partial n}} = - {J_n}. \end{array} \right\} $

式中： $\sigma $为模型内的电导率分布， ${\varphi _0}$为场域内的电位， $\varphi $为边界的电位， ${J_n}$为注入电流密度. 因为求解的目标大多形状不规则，EIT正问题很难有解析解，所以本研究采用有限元法对EIT正问题进行求解.

本研究针对EIT正问题中电导率分布的不确定性对输出电位的影响，因此EIT模型不确定性量化的输出变量为边界电位. 输入不确定性参数为电导率分布，并假定电导率分布服从随机均匀分布.

基于稀疏网格配点法的EIT不确定性量化，即构造显式表达的替代模型使之逼近原隐式表达的EIT数学模型，主要有如下步骤。

1）根据指定的各维精度水平K=[k,···,k]计算满足式（3）的所有多指数组合i=[i₁,···,i_d].

2）根据多指数组合i和各维输入电导率的概率分布 $\rho ({\boldsymbol{\lambda}} )$，计算一维样本点和基函数.

3）采用张量积运算，将一维形式的样本点和基函数扩展到多维形式.

4）利用最小二乘回归法计算多项式函数的系数向量，基于式（4），得到稀疏网格配点法构造的替代模型 $\hat f({\boldsymbol{\lambda}} )$.

5）根据需要求解输出电极电位的均值 $\mu $及标准差s：

(20) $ \mu (\hat f({\boldsymbol{\lambda}} )) = \int\limits_\Omega {\hat f({\boldsymbol{\lambda}} )\rho ({\boldsymbol{\lambda}} ){\rm{d}}{\boldsymbol{\lambda}} } ， $

(21) $ s(\hat f({\boldsymbol{\lambda}} )) = \sqrt {\int\limits_\Omega {{{(\hat f({\boldsymbol{\lambda}} ) - \mu (\hat f({\boldsymbol{\lambda}} )))}^2}\rho ({\boldsymbol{\lambda}} ){\rm{d}}{\boldsymbol{\lambda}} } } . $

改进稀疏网格配点法结合灵敏度分析法，引入新的选点规则. 其EIT不确定性量化过程的主要步骤如下。

1）采用基于方差的全局灵敏度分析法计算EIT模型中各维电导率分布的灵敏度指标，根据所得数据中各灵敏度指标的大小排序来选择模型中各维的精度水平K=[k₁,···,k_d].

2）根据各维的精度水平和各维输入电导率的概率分布 $\rho ({\boldsymbol{\lambda}} )$，计算一维样本点和基函数.

3）基于新的选点规则式（10），通过式（12）、（14）计算多维样本点和基函数.

4）利用最小二乘回归法计算多项式函数的系数向量，得到式（15）原EIT模型的替代模型，求解输出电极电位的统计信息（如均值和标准差）.

基于上述步骤，可以得到EIT替代模型的具体数学表达式、输出变量的均值及标准差，能够定量地研究电导率分布的不确定性对输出电极电位影响的规律特性.

在研究人体头部的电阻抗成像时，采用4层同心圆模型对头部进行仿真，从里到外不同层分别模拟头部的大脑、脑脊液、颅骨和头皮，利用有限元法对整个场域进行单元剖分，所得如图1所示. 4层同心圆头模型中各层的半径r和电导率 $\sigma $^[26]如表1所示. 模型表面均匀放置着16个电极，电流注入为相对激励模式，激励电流加在1、9号电极上，电流从1号电极流入，从9号电极流出，激励电流的频率为50 kHz，幅值为1 mA. 将9号电极的电位作为零参考电位，由有限元法求得的边界上各电极的电位分布如图2所示. 图中， $\varphi $为电极电位，L为电极序号.

在4层同心圆头模型中，输入不确定性参数为头模型中4层组织的电导率值，以表1中所示电导率值为均值，假定各层电导率分布服从变化率为 $ \pm $20%的随机均匀分布. 同时选取1、2、3、4号电极的电位作为输出变量. 分别采用本研究的改进稀疏网格配点法和稀疏网格配点法对EIT正问题中电导率分布的不确定性进行量化研究. 为了便于对比分析，引入MC法和PCE法. PCE法的多项式截断阶数选为3阶（p=3）. MC法作为基准，其计算精度与样本点数N有关，1号电极电位均值 $\mu $的收敛趋势如图3所示. 可知当MC法样本点数达到1×10⁵时，电位均值以趋于收敛，因此，MC法的样本点数选取1×10⁵个.

根据MC法获得的量化结果，采用基于方差的全局灵敏度分析法对4层同心圆头模型中各层电导率进行灵敏度分析，所得灵敏度指标如表2所示. 表中，S₁~S₄分别为1号~4号电极的灵敏度.根据灵敏度分析结果，采用改进稀疏网格配点法选取各层的精度水平K=[k₁,k₂,k₃,k₄]，其中k₁,k₂,k₃,k₄分别代表大脑层、脑脊液层、颅骨层和头皮层的精度水平. 可见，头皮层和颅骨层的灵敏度指标较大，脑脊液层和大脑层的灵敏度指标则很小，表明头皮层和颅骨层的电导率变化对输出电极电位的影响最大，脑脊液层和大脑层的影响较小，采用改进方法选取的一组精度水平可以为K=[3,3,4,4]. 模型中各层电导率的灵敏度指标大小排序为脑脊液层<大脑层<颅骨层<头皮层，采用改进方法选取的另一组精度水平可以为K=[3,2,4,5]. 稀疏网格配点法只能对各层变量赋予相同的精度水平，选取K=[4,4,4,4]、[5,5,5,5]进行对照. 基于改进稀疏网格配点法、稀疏网格配点法、MC法和PCE法的计算结果，1、2、3、4号电极电位 $\varphi $的概率密度函数（probability density function，PDF）分布分别如图4~7所示，各电极电位的统计信息如表3所示. 表中， $\mu $为均值， $s$为标准差，N为样本点数. 可知，当头部各层电导率均匀变化时，输出电极电位的分布规律基本一致趋于正态分布. 由4种UQ方法计算得到的1~4号电极电位的均值与图2中相应电极上的电位分布的有限元法结果基本接近. 电极电位的标准差与电极位置有关，随着各电极位置逐渐远离1号激励电极，各电极电位的标准差在逐渐减小.

将MC法的仿真结果作为基准，定义3种相对于MC法的误差， ${\varepsilon _1}$、 ${\varepsilon _2}$分别为均值误差和标准差误差， ${\varepsilon _3}$为PDF误差，各误差的表达式为

(22) $ {\varepsilon _1} = \left| {\hat \mu - \mu '} \right|, $

(23) $ {\varepsilon _2} = \left| {\hat s - s'} \right|, $

(24) $ {\varepsilon _3}{\rm{ = }}\sqrt {{{\sum\limits_{n = 1}^N {\left( {\hat f({\theta _n}) - f'({\theta _n})} \right)} } \mathord{\left/ {\vphantom {{\sum\limits_{n = 1}^N {\left( {\hat f({\theta _n}) - f'({\theta _n})} \right)} } {\sum\limits_{n = 1}^N {{{f'}^2}({\theta _n})} }}} \right. \kern-\nulldelimiterspace} {\sum\limits_{n = 1}^N {{{f'}^2}({\theta _n})} }}} . $

式中： $\;\mu '$、 $s'$、 $f'({{\theta} _n})$分别为MC法的均值、标准差和PDF； $\;\hat \mu $、 $\hat s$、 $\hat f({{\theta} _n})$分别为其余3种UQ方法的均值、标准差和PDF. 各误差值越小表明该UQ方法越精确，特别是误差 ${\varepsilon _3}$反映PDF曲线的拟合程度，其值越小，说明所获得的计算结果越接近MC法. 1号电极点处稀疏网格配点法、改进稀疏网格配点法和PCE法的各误差值如表4所示.

由表3、4可知，对于稀疏网格配点法，随着各维精度水平的提高（从 $ {\boldsymbol{K}} = \left[ {4,4,4,4} \right] $变化到K=[5,5,5,5]），各误差值虽均有所降低，但计算所需的样本点数增长显著（从401变化到1105）. 对于本研究所提的改进方法，由于是根据输入变量各自不同的重要性程度来分配各维精度水平（从K=[3,3,4,4]变化到K=[3,2,4,5]），其样本点数增长得更加平稳（从385变化到845），并且误差值明显减少，这种优势随着模型维数的增加变得更加突出. 进一步来看，采用K=[3,2,4,5]的改进方法与采用K=[5,5,5,5]的稀疏网格配点法（或采用K=[3,3,4,4]的改进方法与采用K=[4,4,4,4]的稀疏网格配点法），不同精度水平下的这2种方法可以获得几乎相同的精确结果，即2种方法的各误差值均相差很小，但改进方法的计算成本更低，其样本点数少于稀疏网格配点法，因此与稀疏网格配点法相比，本研究改进方法可以用更少的样本点数获得与稀疏网格配点法相当精度的计算结果. 采用最少样本点数的K=[3,3,4,4]改进方法与PCE法相比，发现PCE法的各误差值和样本点数均高于改进方法，表明改进方法能够以更少的样本点数获得比PCE法更加精确的结果. MC法作为基准，计算结果最为准确，但其样本点数相比上述方法都高出几个数量级.

（1）将不确定性量化的概念引入EIT中，研究电导率分布的不确定性对EIT正问题输出电位的影响. 以4层同心圆头模型为算例，发现当各层电导率均匀变化时，输出电位的分布呈现正态分布.

（2）为了进一步提升EIT不确定性量化方法的计算效率，在传统稀疏网格配点法的基础上，结合灵敏度分析提出改进的稀疏网格配点法. 改进方法引入新的选点规则，不同于传统稀疏网格配点法只能给各维变量指定相同的精度水平，改进方法可根据各维变量不同的重要性程度合理配置各维精度水平,将非重要维度上多余的样本点分配到更重要的维度上，使样本点可以更合理地配置.

（3）通过与传统稀疏网格配点法、MC法和PCE法相比，发现改进方法具有较高的计算精度，且计算量更小,计算效率更高，更适应于EIT不确定性问题的量化研究. 在EIT中，正问题的计算结果为逆问题的求解提供依据. 电导率分布的不确定性会使EIT正问题产生不确定的计算结果，影响逆问题的图像重构. 采用不确定性量化方法对影响EIT正问题的不确定性变量进行量化，对进一步优化EIT逆问题的成像算法、降低不确定性变量的影响具有重要意义.

（4）需要指出的是，改进的稀疏网格配点法中各维精度水平的选取是根据灵敏度分析的计算结果按照经验选择的，为了更合理地确定精度水平，可以在灵敏度指标与精度水平间指定一些映射函数，这将在今后研究中做进一步探索.