弥散加权成像(DWI)^[1-4]是通过施加多方向的弥散磁敏感梯度将每一弥散敏感梯度下组织中水分子的扩散运动转化为图像的灰度或其他参数信息的技术，该技术可以用于活体心脏心肌纤维建模^[5]、大脑纤维、中枢神经系统病变、肝脏纤维化^[6]等疾病的辅助诊断。随着远程医疗诊断技术的普及，使得越来越多的DWI数据应用于远程诊断和科学研究中，也让原本在医院单机上存储和使用的DWI图像需要通过网络进行传输和使用。为了防止未授权的DWI图像被使用或遭受篡改，提出对DWI数据进行保护具有重要的实用价值和研究意义。数字水印技术^[7-8]作为保护数字媒体版权信息的一种重要技术，可以通过在DWI图像中加入特定版权标识信息从而增强DWI数据的信息安全性。虽然当前医学图像的版权保护通常采用可逆水印技术^[9]，但大多数可逆水印面对各种有意或无意的攻击时，几乎不具备鲁棒性。因此部分学者尝试研究具有鲁棒可逆性质的水印算法。

Zhang等人^[10]提出将图像通过对冗余小波变换选取中频子带系数进行随机投影压缩，使用压缩矩阵相邻系数进行水印嵌入。该算法对部分几何攻击与信号攻击具有较好的鲁棒性，但需要携带观测矩阵使得需要传输的图像数据增大。Bekkouch等人^[11]提出将图像的离散小波系数高频分量进行离散余弦变换，所获取的离散余弦系数采用奇异值分解，通过在奇异值矩阵中进行修改，从而完成水印的嵌入。该算法在旋转、裁剪、JPEG压缩、椒盐噪声均能较好地提取原始水印信息。An等人^[12]提出使用属性启发式像素调整(PIPA)算法对子带系数进行调整，用于防止水印嵌入导致像素值溢出的情况。该算法在不同类型的图像中嵌入时，均具有较好的可视质量且提出使用聚类算法应用于水印算法的提取。

由于鲁棒可逆水印算法在提取水印信息时，通常需要携带观测矩阵，且图像在受到攻击时提取的图像无法保证无损可逆。因此使得一些医学图像信息保护领域的学者们开始探索和尝试将鲁棒零水印算法的思想应用于医学图像版权信息保护中。隋淼等人^[13]提出一种基于离散余弦变换(DCT)的医学图像鲁棒零水印算法。该算法有效利用DCT低中频系数构建鲁棒零水印信息且具有较强的鲁棒性。但零水印^[13-15]无法对图像的版权进行有效保护。

在与DWI成像领域相关学者讨论研究后得到如下的研究前提条件：在DWI计算转换得到的弥散张量成像(DTI)模型中，通过将计算参数改变范围有效地进行控制，从而保证计算得到的DTI模型对医生的诊断信息不会产生较为明显的差别。在此前提下本文提出将鲁棒水印算法应用于DWI医学图像版权信息保护中。本文提出一种基于DWI图像鲁棒水印算法。该算法首先使用最大类间方差分割算法和面积控制阈值将满足条件的切片作为待嵌入区域，并将该区域拟合为矩形，其次对待嵌入区域使用整数小波变换获取其低频子带系数，使用固定步长对低频子带系数进行直方图统计，并对系数直方图簇的比值关系进行调整从而嵌入鲁棒水印信息，最后对嵌入水印后图像的DWI表观系数与DTI弥散张量值比值关系作为可逆密钥。再一次对已嵌入图像进行加密，进一步增强本算法对DWI图像的版权信息保护。

继Hahn^[16]于1956年提出水弥散对磁共振信号具有一定影响之后，Stejskal和Tanner^[17]提出Stejskal-Tanner方程用于定量测量水分子扩散系数, Bihan等人^[18]首次提出将DWI应用于生物组织中，DTI是一种基于DWI技术的发展与深化，是当前唯一通过非侵入性检查可以有效进行脑部白质纤维和心肌纤维等纤维束的跟踪的成像方法，对于揭示脑神经病变和心脏病变的原理具有重要研究意义。Kingsley^[19]系统地介绍了关于张量成像相关参数的数学计算方法，并阐述噪声对于张量成像中相关参数的影响。

常规MRI(magnetic resonance imaging)是一种以T1W1和T2W1为主的磁共振成像技术，主要显示人体器官或组织的形态结构及其信号强度变换。由于传统的MRI需要较长的成像时间，因此出现一些以快速成像序列为主的MRI成像技术，DWI成像技术就是其中一种。

DWI成像技术是对同一断层图像施加单个非零强度值的多方向弥散梯度，在不同方向弥散梯度下影响体内水分子流动并完成待检测部位影像采集。当强度值为零时，所获取的图像等同于传统MRI图像中的T2W1图像。

DTI图像通常采用多个梯度方向获取的DWI图像使用非线性最小二乘法建模计算得到图像中每个体素扩散张量，通过扩散张量进一步计算得到平均弥散程度${D_{{\rm{an}}}}$、扩散各项同性${D_{{\rm{av}}}}$、主轴方向角${\alpha _A}$。因此DTI不仅可以提供器官的像素值信息，还可以提供相关结构信息，如心肌纤维、脑部神经纤维走向等结构信息。

DWI和DTI的定义和计算过程在文献[15-18]中有详细介绍，本文不做过多赘述，现将本文要用到的关于DIW和DTI的各个参数和计算公式在表 1中进行说明。根据文献[15-18]相关计算，表 1列出DTI建模中将涉及的相关符号。

本文所述算法通过将${{D}_{\text{av}}}$、${{D}_{\text{an}}}$、${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}$的参数调整幅值分别控制在一定范围内，从而保证DTI模型的该变量在可接受的范围内。图 1是通过多幅DWI图像构建出来的DTI图像结构图，通常使用椭圆球体或圆柱体表征其结构，长轴方向为${{\alpha }_{\text{A}}}$角度所指示的方向。

最大类间方差分割法是一种基于全局求解最佳分割阈值的算法。假设图像为$\mathit{\pmb{I}}$，图像平均灰度${{I}_{\text{avg}}}$，前景平均灰度${{I}_{\text{f}}}$，背景平均灰度${{I}_{\text{b}}}$，属于前景像素点的概率${{P}_{\text{f}}}$和属于背景像素点概率${{P}_{\text{b}}}$，通过求取类间方差为$V$，最大类间方差为${{V}_{\text{max}}}$。则通过式(1)和式(2)，可以求出使得前景与背景方差最大值作为最佳分割阈值${{T}_{1}}$，即类间方差最大值

$ V = {P_{\rm{f}}} \times {\left( {{I_{\rm{f}}} \leftarrow {I_{{\rm{arg}}}}} \right)^2} + {P_{\rm{b}}} \times {\left( {{I_{\rm{b}}} - {I_{{\rm{arg}}}}} \right)^2} $

(1)

$ {V_{\max }} = \arg \max (V) $

(2)

将每一切片在不同弥散方向下的图像通过最大类间方差分割法进行分割，从而提取出图像的前景区域，但提取的前景区域并不一定是规则形状，因此本文将前景区域内的所有像素点的$x$坐标值和$y$坐标值进行排序，提取$\left( {{x}_{\min }}, {{y}_{\min }} \right)$，$\left( {{x}_{\max }}, {{y}_{\min }} \right)$，$\left({{x}_{\min }}, {{y}_{\max }} \right)$，$\left( {{x}_{\max }}, {{y}_{\max }} \right)$4个点将前景区域拟合为矩形并作为待嵌入区域。假设DWI图像切片数为$nu{{m}_{\text{c}}}$，弥散梯度方向数为$nu{{m}_{\text{d}}}$，任意切片中前景像素点个数为$n$。通过面积控制阈值$T$获取最大可嵌入区域计算步骤如下：

1) 本文读取一组患者DWI图像$\mathit{\pmb{I}}=\left\{ {{\mathit{\pmb{I}}}^{c, d}}\left| c\in \left[ 1, nu{{m}_{\text{c}}} \right], \right.d\in \left[ 1, nu{{m}_{\text{d}}} \right] \right\}$，DWI图像的长和宽为$m$、$n$，$m\times n$为图像的面积。

2) 对步骤1)中所有的切片使用最大类间方差分割法提取扩散方向为${{B}_{1}}$的前景区域，并计算其面积${{S}_{\text{ROI}}}$，设定前景区域${{S}_{\text{ROI}}}$与DWI图像面积之比${{P}_{\text{ROI}}}$大于面积控制阈值${{T}_{\text{s}}}$，则选取该范围切片作为待嵌入切片。${{P}_{\text{ROI}}}$的计算为

$ {P_{{\rm{ROI}}}} = \frac{{{S_{{\rm{ROI}}}}}}{{m \times n}} $

(3)

3) 获取所有切片多个扩散方向的前景区域像素点坐标位置集合${{\mathit{\pmb{L}}}_{\text{c}}}$，将每一个切片的前景区域的$\left[ \left( {{x}_{\min }}, {{y}_{\min }} \right), \left( {{x}_{\max }}, {{y}_{\min }} \right), \left( {{x}_{\min }}, {{y}_{\max }} \right), \left( {{x}_{\max }}, {{y}_{\max }} \right) \right]$4个像素点坐标组成一个集合${{\mathit{\pmb{L}}}_{\text{c}}}$，${{\mathit{\pmb{L}}}_{\text{c}}}=\left\{ \left( x_{d, c}^{\min }, y_{d, c}^{\min } \right), \left( x_{d, c}^{\max }, y_{d, c}^{\min } \right), \left( x_{d, c}^{\min }, y_{d, c}^{\max } \right), \left( x_{d, c}^{\max }, y_{d, c}^{\max } \right) \right\}$，作为所有切片下的待嵌入区域, 将该待嵌入区域记为${\mathit{\pmb{I}}}_{f}^{c, d}$。通过该方法提取到的待嵌入区域与原始切片图像如图 2所示。

灰度统计直方图是指通过统计一幅图像中所出现的灰度级个数，反映每个灰度级出现的频率。由于在DWI图像中${{B}_{1}}$方向为不施加弥散梯度时获取的MRI图像，与其他施加弥散梯度方向${{B}_{2}}$~${{B}_{13}}$图像相比具有较大灰度值，在遭受攻击时, 其灰度直方图变形较大，因此本文选取已施加弥散梯度方向梯度图像$\left( {{B}_{2}}, {{B}_{3}}, \ldots , {{B}_{13}} \right)$用于小波系数统计。

文献[20]中已证明，当图像遭受几何攻击时，虽然像素点的空间位置发生改变，但其统计直方图的主体形状不会发生较大改变。由离散小波变换基本原理可知，图像的低频子带系数是原图在最大尺度与最小分辨率下的最优逼近。因此低频子带系数的统计直方图对于几何攻击也应具有较好的几何不变性。

本文提出基于整数小波变换的嵌入算法，其具体步骤如下：

输入待嵌入DWI图像${\mathit{\pmb{I}}}_{f}^{c, d}$, 未调制的一维二值水印向量$\mathit{\pmb{W}}=\left\{ {{w}_{1}}, {{w}_{2}}, \ldots {{w}_{{{n}_{w}}}} \right\}$，水印的长度为${{n}_{w}}$，区域面积控制参数${{T}_{\text{s}}}$，子带系数范围阈值$\alpha $。

输出嵌入后DWI图像$\tilde{\mathit{\pmb{I}}}_{f}^{c, d}$。

1) 使用随机函数生成一个二值混沌向量$\mathit{\pmb{O}}=\left\{ {{o}_{1}}, {{o}_{2}}, \ldots {{o}_{{{n}_{w}}}} \right\}$与1维二值水印向量按位异或生成水印信号$W_{\left( k \right)}^{'}=\left\{ o_{1}^{'}, o_{2}^{'}, \ldots o_{{{n}_{w}}}^{'} \right\}, k\in \left[ 1, {{n}_{w}} \right]$。

2) 将DWI图像中的每一切片的不同扩散方向${\mathit{\pmb{I}}}_{f}^{c, d}$作为一个嵌入单元进行整数小波变换获取低频子带系数$\mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right)$, 由于$\mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right)$的系数间隔较小，为控制水印嵌入量，本文以步长$\delta <1$进行低频系数统计。为将水印嵌入到直方图主体区域，因此求出该切片下的低频子带系数最大值，记为$\arg \max \left( \mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right) \right)$，以及低频子带系数平均值，记为$\arg \;\text{avg} \left( \mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right) \right)$。本文利用统计子带系数平均值与子带系数最大值比值大小, 用于设定待嵌入的子带系数范围，其系数范围如式(4)所示

$ \begin{array}{*{20}{c}} {{R_{\rm{c}}} = \left[ {(1 - \alpha )\arg avg \left( {\mathit{\boldsymbol{R}}_{LL}^{c,d}(u,v)} \right),} \right.}\\ {\left. {(1 + \alpha )\arg avg \left( {\mathit{\boldsymbol{R}}_{LL}^{c,d}(u,v)} \right)} \right]} \end{array} $

(4)

式中，$\alpha \text{=}\frac{\arg \text{avg}\left( \mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right) \right)}{\arg \max \left( \mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right) \right)}$，为保证水印鲁棒性, 设置水印的待嵌入区域${{R}_{\text{c}}}\ge 2{{n}_{w}}$。

3) 将$\mathit{\pmb{R}}_{LL}^{c, d}$中的系数进行两两分组作为一个簇，假设某相邻两灰度统计直方图灰度范围值为$\mathit{\pmb{R}}_{LL}^{c, d1}$、$\mathit{\pmb{R}}_{LL}^{c, d2}$时，其对应的灰度统计个数为${{H}_{1}}$、${{H}_{2}}$，本文算法通过调整统计灰度个数嵌入水印信号，调节强度$T$值可以平衡水印鲁棒性和透明性。

当$w{{'}_{\left( k \right)}}==1$，且${{H}_{1}}$的个数比的个数小于$T$时，按照2.3节所述规则选取系数值为$\mathit{\pmb{R}}_{LL}^{c, d}\left( {{u}_{2}}, {{v}_{2}} \right)$的${{n}_{\text{rand}}}$个子带系数修改至$\mathit{\pmb{R}}_{LL}^{c, d}\left( {{u}_{1}}, {{v}_{1}} \right)$系数范围内，使得修改后灰度统计个数$H_{1}^{'}$、$H_{2}^{'}$关系为$\frac{H_{1}^{'}}{H_{2}^{'}}>T$。

当$w{{'}_{\left( k \right)}}==0$，且${{H}_{2}}$的个数比${{H}_{1}}$的个数小于$T$时，选取系数值为$\mathit{\pmb{R}}_{LL}^{c, d}$的${{n}_{\text{rand}}}$个子带系数修改至$\mathit{\pmb{R}}_{LL}^{c, d2}$系数范围内，使得修改后灰度统计个数$H_{1}^{'}$、$H_{2}^{'}$关系为$\frac{H_{1}^{'}}{H_{2}^{'}}<T$。

4) 将修改后子带系数使用整数小波变换逆变换，生成已嵌入水印图$\bar{I}_{f}^{c, d}$。

由表 1中公式可知，DWI图像在修改灰度值后对DTI求解到的结构参数有较大影响。为保证将嵌入后DWI图像所求得的DTI结构参数尽可能控制在一定范围内，本文提出一种基于标记的待调整系数选取规则。其具体步骤如下：

1) 将输入小波子带系数$\mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right)$按照待修改小波系数范围进行统计，并记录该范围低频子带系数位置$\mathit{\pmb{l}}\left( {{u}_{k}}, {{v}_{k}} \right)$, 这样的系数个数有${{n}_{c}}$个，且${{n}_{c}}>{{n}_{w}}$。

2) 获取$\mathit{\pmb{l}}\left( {{u}_{k}}, {{v}_{k}} \right)$处的低频子带系数$\mathit{\pmb{R}}_{LL}^{c, d}\left( u, v \right)$，按照2.2节所述嵌入算法进行调整，使得调整后的低频系数值$\mathit{\pmb{R}}_{LL}^{'c, d}\left( u, v \right)$满足水印嵌入算法要求。

3) 对修改后的子带系数进行重构，并计算是否存在有DTI角度参数${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}>5{}^\circ $。若有则恢复其低频子带系数值；若${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}\le 5{}^\circ $，将修改后值存入修改子带系数矩阵，其具体的流程图如图 3所示。

利用表 1中的公式，可对灰度值和张量矩阵之间的计算结果在特定情况下进行推导，因此本文算法对施加梯度矢量$\mathit{\pmb{b}}_{ij}^{\mathit{\pmb{v}}}$的灰度值${{S}_{\text{b}}}$进行修改，称为针对灰度值的梯度平移，假设对灰度值进行放大2倍处理后，修改不同弥散梯度方向矢量$\mathit{\pmb{b}}_{ij}^{\mathit{\pmb{v}}}$值所对应的灰度值$\mathit{\pmb{S}}'$，如式(5)所示。

$ S' = \left( {\begin{array}{*{20}{c}} {{S_1}}\\ \vdots \\ {2 \times {S_{12}}}\\ {2 \times {S_{13}}} \end{array}} \right) $

(5)

弥散梯度相关参数在采集时已经确定，即弥散梯度强度${{b}_{ij}}\text{=}0$时${{B}_{1}}$对应${{S}_{1}}$，由于${{S}_{1}}$灰度值较大，将式(5)代入表 1的公式中可得

$ \mathit{\boldsymbol{S'}} = \left( {\begin{array}{*{20}{c}} {{S_1}}\\ \vdots \\ {2 \times {S_{12}}\exp \left( {\sum\limits_{i = {x_{12}},{y_{12}},{z_{12}}} {\sum\limits_{j = {x_{12}},{y_{12}},{z_{12}}} {{b_{ij}}{{D'}_{ij}}} } } \right)}\\ {2 \times {S_{13}}\exp \left( {\sum\limits_{i = {x_{13}},{y_{13}},{z_{13}}} {\sum\limits_{j = {x_{13}},{y_{13}},{z_{13}}} {{b_{ij}}{{D'}_{ij}}} } } \right)} \end{array}} \right) $

(6)

由式(6)可知，不同弥散梯度方向下的灰度值发生改变，但${{B}_{0}}$方向下弥散梯度方向矢量相关参数无法修改，因此其灰度值不发生改变。将${{S}_{b}}={{S}_{1}}\text{exp}\left(-\sum\limits_{i\in DL}{\sum\limits_{j\in DL}{{{b}_{ij}}{{D}_{ij}}}} \right), DL=\left\{ x, y, z \right\}$进行变形可获得关系如式(7)所示。式(7)中${{b}_{ij}}$为施加的弥散方向矢量的值，${{D}_{ij}}$为张量矩阵中的一个张量值。

$ \ln \frac{{{S_1}}}{{{{S'}_b}}} = \sum\limits_{i \in {\mathit{\boldsymbol{D}}_L}} {\sum\limits_{j \in {\mathit{\boldsymbol{D}}_L}} {{b_{ij}}{{D'}_{ij}}} } ,{\mathit{\boldsymbol{D}}_L} = \left\{ {x,y,z} \right\} $

(7)

联立式(7)和表 1中公式，并代入式(5)中仅存在一个扩散方向时的值可得，当且仅当有一个扩散方向时其改变值与原值的关系如式(8)所示

$ \begin{array}{*{20}{c}} {\ln \left( {\frac{{{S_1}}}{{S_b^\prime }} - \ln \frac{{{S_1}}}{{{S_b}}}} \right) = \ln \left( {\frac{{{S_1}}}{{S_b^\prime }} \times \frac{{{S_b}}}{{{S_1}}}} \right) = \ln \left( {\frac{{{S_1}}}{{2 \times {S_b}}} \times \frac{{{S_b}}}{{{S_1}}}} \right) = }\\ { - \ln 2 + \ln 1 = \sum\limits_{i \in DL} {\sum\limits_{j \in DL} {\left( {{b_{ij}}D_{ij}^\prime - {b_{ij}}{D_{ij}}} \right)} } ,}\\ {{\mathit{\boldsymbol{D}}_L} = \{ x,y,z\} } \end{array} $

(8)

将其扩展至多方向时，由式(8)可推导出

$ \ln \frac{{{S_1}}}{{S_b^\prime }} = \left( {\begin{array}{*{20}{c}} {\ln \frac{{{S_1}}}{{{S_2}}} - \ln 2}\\ \vdots \\ {\ln \frac{{{S_1}}}{{{S_{12}}}} - \ln 2}\\ {\ln \frac{{{S_1}}}{{{S_{13}}}} - \ln 2} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\ln \frac{{{S_1}}}{{{S_2}}}}\\ \vdots \\ {\ln \frac{{{S_1}}}{{{S_{12}}}}}\\ {\ln \frac{{{S_1}}}{{{S_{13}}}}} \end{array}} \right) - \left( {\begin{array}{*{20}{c}} {\ln 2}\\ \vdots \\ {\ln 2}\\ {\ln 2} \end{array}} \right) $

(9)

由于${{b}_{ij}}$的值并不随着灰度值的改变而改变，因此当弥散梯度方向矢量下灰度值增大$n$倍时，同理可推论得出单点对应的张量值应减小为原值$-\ln n$。因此当已知当前的张量关系与表观系数增加倍数为$n$，可以通过以上推导求得原始扩散图像灰度值，也可证明在此情况下，张量矩阵的变换与原张量矩阵之间具有一定关联，可利用三重积公式进行证明，其具体证明如附录所示。因此，本文利用此推论提出一种基于DWI图像可逆加密算法，通过该加密算法中密钥${{k}_{n}}$对已嵌有水印的DWI图像$\mathit{\pmb{I}}_{f}^{c, d}$再次加密，获得可逆密钥加密后图像$\widehat{{\mathit{\pmb{I}}}_{f}^{c, d}}$。在对嵌入水印信息且生成的DTI图像解密时，利用上述推导获取的密钥${{k}_{n}}$可直接计算获得无损的DWI灰度值矩阵$\mathit{\pmb{S}}$。

1) 利用面积控制阈值${{T}_{\text{s}}}$提取已嵌入水印的切片的拟合区域$\widehat{{\mathit{\pmb{I}}}_{f}^{c, d}}$。

2) 由步骤1获得嵌有可逆密钥与水印图像$\widehat{{\mathit{\pmb{I}}}_{f}^{c, d}}$，从接收到的图像文件中提取到可逆密钥${{k}_{n}}$对带弥散矢量加权图像灰度使用下式进行解密，从而获得无可逆密钥图像$\widehat{{\mathit{\pmb{I}}}_{f}^{c, d}}$。

$ \mathit{\boldsymbol{\widetilde I}}_f^{c,d} = \frac{{\widehat {\mathit{\boldsymbol{I}}_f^{c,d}}}}{{{k_n}}} $

(10)

1) 通过2.5.1节计算获得可逆密钥，使用该密钥可对加密后待提取水印信号的图像解密，获得图像$\tilde{\mathit{\pmb{I}}}_{f}^{c, d}$。

2)对步骤1)中获取到的$\tilde{\mathit{\pmb{I}}}_{f}^{c, d}$进行整数小波变换，获取低频子带系数，使用子带系数范围$\tilde{\mathit{\pmb{R}}}_{LL}^{c, d}\left( u, v \right)$提取嵌入水印子带系数范围区间，将低频子带以步进为$\delta $进行低频系数统计。假设有两个相邻统计子带系数为$\tilde{\mathit{\pmb{R}}}_{LL}^{c, d}\left( {{u}_{1}}, {{v}_{1}} \right)$、$\tilde{\mathit{\pmb{R}}}_{LL}^{c, d}\left( {{u}_{2}}, {{v}_{2}} \right)$，其灰度统计个数为$H_{1}^{'}$、$H_{2}^{'}$。

当相邻子带系数统计直方图之比$\frac{H_{1}^{'}}{H_{2}^{'}}\ge 1$时，该相邻统计直方图组中嵌入水印值$\mathit{\pmb{w}}_{\left( c, k \right)}^{''}=1$；否则当对应的相邻统计直方图簇个数比值$\frac{H_{1}^{'}}{H_{2}^{'}}<1$时，对应的统计直方图组中所嵌入的水印值$w_{\left( c, k \right)}^{''}=0$。

3) 由步骤2)获得$num_{c}^{'}$个版本的$w_{\left( c, k \right)}^{''}=0$，其中$num_{c}^{'}\le nu{{m}_{\text{c}}}$，当且仅当${{T}_{\text{s}}}\le \arg \ \min \ {{P}_{\text{ROI}}}$时，$num_{c}^{'}=nu{{m}_{\text{c}}}$。然后利用投票公式(如式(11))获得加密水印序列$\mathit{\pmb{w}}_{\left( k \right)}^{''}$，将加密水印序列与用水印嵌入调制时使用的同一个混沌序列进行按位异或获得提取的水印信息序列$\mathit{\pmb{w}}_{\left( k \right)}^{\text{end}}$。

$ \mathit{\boldsymbol{w}}_{\left( k \right)}^{{\rm{end}}} = \left\{ \begin{array}{l} 1\;\;\;\;\;\sum\limits_{c = 1}^{nu{{m'}_c}} {{{w''}_{\left( {c,k} \right)}}} \ge \frac{{nu{{m'}_c}}}{2}\\ 0\;\;\;\;\;aaa \end{array} \right. $

(11)

在本文中使用在不同切片中反复重嵌入相同水印值，从而有效降低提取水印时误判率，提高提取后水印的鲁棒性。

本文使用MATLAB 2017b作为实验平台进行实验，本实验室于不同设备中采集到约42位病人，每位病人采集25~54个切层，采集方向数量为13~65个方向，共计约20 000余幅图像。在DWI图像数据库中任选一组病人的DWI图像用于实验结果的说明。通过多次实验后, 本文设置面积阈值${{T}_{\text{s}}}$为0.7可以较有效地去除心脏图像中的心尖与心底图像。为保证子带系数具有较好的集合不变性，选取子带系数范围阈值$\alpha $为0.6用于提取子带系数统计图主体部分。其中嵌入强度是平衡算法鲁棒性与不可见性的关键指标，为平衡指标通过反复实验选取嵌入强度$T$为4，在DWI图像中低频子带系数变化约为每隔0.25会出现一个子带系数值，为平衡鲁棒性与角度参数，选取步长$\delta $为0.5。

在DTI成像中，主轴方向反映体素的纤维走向是医生用于诊断纤维病变的关键依据，因此保证纤维主方向变化夹角在一定范围内是本文算法中对于医生诊断不影响的最大保证。使用表 1中公式计算原始主轴方向矢量$\mathit{\pmb{\varepsilon}} =\left( {{\varepsilon }_{x}}, {{\varepsilon }_{y}}, {{\varepsilon }_{z}} \right)$与改变后的单位主方向矢量$\mathit{\pmb{\varepsilon}} '=\left( \varepsilon _{x}^{'}, \varepsilon _{y}^{'}, \varepsilon _{y}^{'} \right)$，计算夹角改变量${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}$、${{D}_{\text{av}}}$、${{D}_{\text{an}}}$三个参数改变前后超出额定指标在多个切片中统计的个数，统计结果见表 2所示。

由表 2可知，利用本文所述算法，对33个切层进行水印嵌入后，较少出现像素点嵌入前后夹角${{\alpha }_{\text{A}}}$变化超过5°像素点，并且平均弥散程度${{D}_{\text{av}}}$和各项同性${{D}_{\text{an}}}$变换量很小，满足DWI向DTI模型转换后对DTI中各参数改变量的要求。

任意选取一组图像中某一切层弥散加权图像对各参数进行计算, 如图 4所示。图 4(a)(d)(g)为未嵌入图像前的各参数热图，图 4(b)(e)(h)为嵌入水印后各参数热图，图 4(c)(f)(i)为嵌入前后各参数变化量绘制的对应热图。由图 4可知，各参数虽存在少量变化，但没有超出临床医生所要求的标准，且出现位置皆在非心肌纤维处，因此本算法可以适用于临床。

为评价本文所提算法的性能，本文将从视觉质量和水印信息鲁棒性两个指标对算法性能进行客观分析和测量。

峰值信噪比(PSNR)是一种通过客观统计进行水印嵌入前后图像质量标准的方法。当嵌入前后的图像的差异越小时，PSNR的值越大。通常认为当PSNR大于35 dB时图像的差异通过肉眼已经无法进行判别。

$ PSNR = 10{\rm{lg}}\left( {\frac{{m \times n \times {2^{12}} \times {2^{12}}}}{{\sum\limits_{i = 0}^{m - 1} {\sum\limits_{j = 0}^{n - 1} {{{\left( {\mathit{\boldsymbol{\tilde I}}_f^{c,d}(i,j) - \mathit{\boldsymbol{I}}_f^{c,d}(i,j)} \right)}^2}} } }}} \right) $

(12)

归一化相关系数(NC)用于判断在嵌入图像遭受攻击前后水印信号的差别，以反映水印嵌入算法的鲁棒性，具体为

$ NC = \frac{{\sum\limits_{i = 1}^{nw} {{{\mathit{\boldsymbol{W'}}}_{\left( k \right)}}\left( i \right)} \times {\mathit{\boldsymbol{w}}^{{\rm{end}}}}\left( i \right)}}{{\sum\limits_{i - 1}^{nw} {{{\left( {{{\mathit{\boldsymbol{W'}}}_{\left( k \right)}}\left( i \right)} \right)}^2}} }} $

(13)

旋转攻击是在图像攻击中一种常见的几何攻击方式，它利用旋转变换，改变原有灰度所处的空间位置，其中旋转变换分为旋转变换和带裁剪的旋转变换。其中带裁剪的变换是指由于对原始图像进行旋转攻击后图像尺寸发生变化，因此为保证旋转后图像尺寸与原始图像尺寸保持一致，对旋转后图像进行边缘裁剪。在本文所提算法的实验中，当旋转角度在1°~ 5°时，$NC$= 1。当旋转角度为15°时，$NC$ = 0.98，虽然随着旋转角度的增大，本文算法的鲁棒性有所下降，但始终保持一个较好的质量, 因此本算法对于一定强度的几何攻击具有较好的鲁棒性。具体实验结果如表 3所示。

缩放攻击是对图像进行大小缩放变换后生成新的不同或与原图像大小相同的图像。该攻击通常分为等比例缩放和不等比例缩放。本文中对嵌入水印的DWI图像进行0.7, 0.9, 1.2和1.5倍等比例缩放攻击。当缩放攻击为0.9倍时$NC$值等于0.67，当缩放攻击为1.5倍时$NC$值等于0.81。具体实验结果如表 4所示。

噪声攻击是指在原始图像信号中加入一定的噪声信号，使得原始图像信号遭到篡改。本文使用椒盐噪声和高斯噪声作为主要的噪声攻击方式，本文选择对所有切层下所有方向图片分别加入一定量的高斯噪声和椒盐噪声。实验结果如表 5所示。通过实验可知，本文所提算法在不同噪声攻击下的$NC$值皆大于0.79，在方差分别为0.01和0.02的椒盐噪声攻击时，本文所提算法的$NC$值可以达到1。由此可见，本算法对常见的噪声攻击具有较好的鲁棒性。

通常在远程医疗传输接收端为减少信道中存在高斯噪声对医生诊断的影响会使用低通滤波器对图像进行预处理，因此在DWI的水印中能抵御低通滤波有助于在图像传输后在接收端依旧可以提取版权信息。本算法在进行高斯低通滤波攻击时，NC值为1，因此本文所提算法针对高斯低通滤波具有很强的鲁棒性。

信号增强攻击是指通过亮度、对比增强或直方图均衡化对原始图像进行篡改，从而影响原始水印信息于图像之间的同步性，影响水印信息的有效性。通过表 6可知本文所提算法在直方图均衡化和对比度增强10%时NC值均能保持0.8左右，因此本文所提算法对信号增强攻击具有一定鲁棒性。

本文选取切片数量为54层，施加弥散梯度方向数为12，大小为128×128像素的人体弥散加权心脏图像作为测试图像，但当前尚无针对弥散加权成像的不可见鲁棒水印算法。本文选取文献[20]的鲁棒水印算法作为对比实验组，将该算法应用于本文所使用的弥散加权图像中，对比两种算法在弥散加权成像中的结果。

本文首次将鲁棒水印应用于DWI图像的版权保护中。(为此本文提出相关算法可适用于弥散加权成像的前程)：需要使得嵌入水印后的图像不仅满足在视觉上使得人眼无法感知，还需要使得嵌入后DWI建模形成的DTI参数中相应纤维主方向变化度数${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}$不得超过5°，平均弥散程度与各向异性改变量不得超过0.01，因此在不同算法的比较中首先针对该算法的适用性进行比较。

将本文算法与文献[20]中提出的低频水印算法进行对比。图 5(a)(c)(e)为原始切层图像在${{B}_{1}}$方向下的图像，图 5(b)(d)(f)为文献[20]所提算法在不同切层下图像的角度变化情况，使用红色“+”对主轴方向改变量${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}$超过5°的坐标点进行标记。使用文献[20]所提出算法，在心脏图像中多处出现角度偏差较大情况，且部分出现在心肌纤维处，该算法可能会严重影响心肌纤维处原始纤维走向诊断。因此本文算法与文献[20]相比，在角度控制方面性能更好。

将本算法与文献[20]算法对心脏非心尖、心底区域所在的33个切片，先计算每一个切片的PSNR值，然后再求平均得到的整个嵌入切片的平均PSNR值进行比较，其结果如表 7所示。

根据反复实验，在该情况下，水印信号无法被人眼直接察觉，但从表 8可知，在同样的切片数和方向数下，本文所提算法在视觉评价指标中明显好于文献[20]中的算法。

在测试鲁棒性的实验中，将本文算法与文献[20]低频算法进行鲁棒性对比试验，测试在高斯噪声、旋转等情况下的鲁棒性比较。实验结果如表 8所示。

由表 8可知本文算法在旋转攻击、缩放攻击、高斯滤波、椒盐噪声等情况下，相比文献[20]算法具有较好的鲁棒性。

在参数改变量对比实验中，本文将文献[20]中算法与本文算法应用于图像库中同一组图像数据后，统计各参数变量在所有切片该变量的均值。由表 9所示，文献[20]所提出的鲁棒水印算法，角度参数改变量${{\alpha }_{\text{A }\!\!\_\!\!\text{ c}}}$超过5°的个数为149个，本文算法角度改变平均个数为1个；在反映纤维各项同性参数上本文略好于文献[20]；在各项异性参数上，两种算法均不会产生较大影响。

本文提出一种基于直方图统计修改的弥散加权图像鲁棒水印。该算法首先利用最大类间方差法提取患者所有切片中的前景区域，并使用面积阈值选取待嵌入切片。其次将待嵌入切片中所有弥散方向前景区域图像进行整数小波变换，对指定范围内的低频小波区间使用固定步长统计子带系数直方图。然后调整相邻系数统计直方图簇的比值关系完成水印信息的嵌入。最后使用表观系数和弥散张量的关联，设计可逆密钥。该算法在图像遭受几何攻击和噪声攻击时，具有一定鲁棒性且该算法在生成的DTI角度参数中改变较小，因此本算法适用于DWI图像水印嵌入。且本算法在与参考文献[20]所提出的算法进行对比后，在鲁棒性、角度参数和各项同性参数上具有较明显的优势，能够有效控制DTI参数，且本文通过二次嵌入使得未获得密钥时也无法正确提取相应信息，从而增强图像安全性。在本文中对于参数的控制只能通过反馈网络进行，为自适应地进行修改，将提出使用神经网络等算法进行修改。

已知原始DWI表观系数为

$ \mathit{\boldsymbol{S}} = \left( {\begin{array}{*{20}{c}} {{S_1}}\\ \vdots \\ {{S_{N - 1}}}\\ {{S_N}} \end{array}} \right) $

(14)

修改系数变量$v$, 则当前表观系数为

$ \mathit{\boldsymbol{S'}} = \left( {\begin{array}{*{20}{c}} {{S_1}}\\ \vdots \\ {v \times {S_{N - 1}}}\\ {v \times {S_N}} \end{array}} \right) $

(15)

根据文献[10]中所述，使用表 1中公式运用加权线性最小二乘法对矩阵$\mathit{\pmb{B}}$进行拟合求解张量矩阵$\mathit{\pmb{\alpha}} $。

其中，修改表观系数后为

$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}} = diag \left( {\frac{{{{\left( {v \times {S_i}} \right)}^2}}}{{\sigma _i^2}}} \right) $

(16)

因此修改后矩阵${{\mathit{\pmb{B}}}^{\text{T}}}{{\mathit{\pmb{\varSigma}}}^{-1}}\mathit{\pmb{B}}$为

$ {\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}\mathit{\boldsymbol{B}} =\\ \left( {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^N {b_{xxi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {{b_{xxi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{xxi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{xxi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {\sum\limits_{i = 1}^N {{b_{yyi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {b_{yyi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{yyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{yyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{zzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{xyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{xzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{yzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{yzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {b_{yzi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{yzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ { - \sum\limits_{i = 1}^N {{b_{xxi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{yyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{ - 2\sum\limits_{i = 1}^N {{b_{yzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}} } \end{array}} \right) $

(17)

$ \mathit{\boldsymbol{x'}} = \left( {\begin{array}{*{20}{c}} {\ln \left( {{S_1}} \right)}\\ \vdots \\ {\ln \left( {v \times {S_{N - 1}}} \right)}\\ {\ln \left( {v \times {S_N}} \right)} \end{array}} \right) $

(18)

因此联立式(17)和式(18)可知

$ \begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}} \right)\mathit{\boldsymbol{x'}} = }\\ {\left( {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^N {b_{xxi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {{b_{xxi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{xxi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{xxi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {\sum\limits_{i = 1}^N {{b_{yyi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {b_{yyi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{yyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{yyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{2\sum\limits_{i = 1}^N {{b_{zzi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{zzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {{b_{xyi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{xyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {{b_{xzi}}} {b_{yzi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{xzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ {2\sum\limits_{i = 1}^N {{b_{yzi}}} {b_{xxi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{2\sum\limits_{i = 1}^N {{b_{yzi}}} {b_{yyi}}\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{4\sum\limits_{i = 1}^N {b_{yzi}^2} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - 2\sum\limits_{i = 1}^N {{b_{yzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}\\ { - \sum\limits_{i = 1}^N {{b_{xxi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {{b_{yyi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}& \cdots &{ - 2\sum\limits_{i = 1}^N {{b_{yzi}}} \frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}}&{ - \sum\limits_{i = 1}^N {\frac{{{{\left( {S_i^\prime } \right)}^2}}}{{\sigma _i^2}}} } \end{array}} \right)}\\ {\left( {\begin{array}{*{20}{c}} { - {b_{xx1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - {b_{xx2}}\frac{{{{\left( {v \times {S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - {b_{xxN}}\frac{{{{\left( {v \times {S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ { - {b_{yy1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - {b_{yy2}}\frac{{{{\left( {v \times {S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - {b_{{\rm{yyN}}}}\frac{{{{\left( {v \times {S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ \vdots & \vdots & \cdots &{}\\ { - 2{b_{{\rm{xy1}}}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - 2{b_{{\rm{xy2}}}}\frac{{{{\left( {v \times {S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - 2{b_{{\rm{xyN}}}}\frac{{{{\left( {v \times {S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ { - 2{b_{{\rm{xz1}}}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - 2{b_{{\rm{xz2}}}}\frac{{{{\left( {v \times {S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - 2{b_{{\rm{xzN}}}}\frac{{{{\left( {v \times {S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ {\frac{{S_1^2}}{{\sigma _1^2}}}&{\frac{{{{\left( {v \times {S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{\frac{{{{\left( {v \times {S_N}} \right)}^2}}}{{\sigma _N^2}}} \end{array}} \right)} \end{array} $

(19)

$ {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}} = \left( {\begin{array}{*{20}{c}} { - {b_{xxi}}}& \cdots &{ - {b_{xxN}}}\\ { - {b_{yyi}}}& \cdots &{ - {b_{yyN}}}\\ { - {b_{zzi}}}& \cdots &{ - {b_{zzN}}}\\ { - {b_{xyi}}}& \cdots &{ - {b_{xyN}}}\\ { - {b_{xzi}}}& \cdots &{ - {b_{xzN}}}\\ { - {b_{yzi}}}& \cdots &{ - {b_{yzN}}}\\ 1& \cdots &1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} { - {b_{xxi}}}&{ - {b_{yyi}}}&{ - {b_{zzi}}}&{ - 2{b_{xyi}}}&{ - 2{b_{xzi}}}&{ - 2{b_{yzi}}}&1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ { - {b_{xxN}}}&{ - {b_{yyN}}}&{ - {b_{zzN}}}&{ - 2{b_{xyN}}}&{ - 2{b_{xzN}}}&{ - 2{b_{yzN}}}&1 \end{array}} \right) $

(20)

$ \mathit{\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\frac{{S_1^2}}{{\sigma _1^2}}}&0& \cdots &0 \end{array}}\\ {\begin{array}{*{20}{c}} 0&{\frac{{S_2^2}}{{\sigma _2^2}}}& \cdots &0 \end{array}}\\ \vdots \\ {\begin{array}{*{20}{c}} 0&0&0&{\frac{{S_N^2}}{{\sigma _N^2}}} \end{array}} \end{array}} \right) $

(21)

$ {\mathit{\boldsymbol{M}}_1} = \mathit{\boldsymbol{M}} * \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1& \cdots &0&0 \end{array}}\\ {\begin{array}{*{20}{c}} 0&{{v^2}}& \cdots &0 \end{array}}\\ \vdots \\ {\begin{array}{*{20}{c}} 0&0&0&{{v^2}} \end{array}} \end{array}} \right] $

(22)

$ \left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right) = {\mathit{\boldsymbol{B}}^{\rm{T}}} \times \left( {\mathit{\boldsymbol{M}} * \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1& \cdots &0&0 \end{array}}\\ {\begin{array}{*{20}{c}} 0&{{v^2}}& \cdots &0 \end{array}}\\ \vdots \\ {\begin{array}{*{20}{c}} 0&0&0&{{v^2}} \end{array}} \end{array}} \right]} \right) \times \mathit{\boldsymbol{B}} $

(23)

根据三重积公式可得

$ \left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right) = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1& \cdots &0&0 \end{array}}\\ {\begin{array}{*{20}{c}} 0&{{v^2}}& \cdots &0 \end{array}}\\ \vdots \\ {\begin{array}{*{20}{c}} 0&0&0&{{v^2}} \end{array}} \end{array}} \right]*\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{MB}}} \right) $

(24)

对该矩阵求逆可得

$ {{{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}} = {{\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1& \cdots &0&0 \end{array}}\\ {\begin{array}{*{20}{c}} 0&{{v^2}}& \cdots &0 \end{array}}\\ \vdots \\ {\begin{array}{*{20}{c}} 0&0&0&{{v^2}} \end{array}} \end{array}} \right]}^{ - 1}} * {{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{MB}}} \right)}^{ - 1}}} $

(25)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}} = \left( {\begin{array}{*{20}{c}} 1&{{v^2}}&{{v^2}}&{{v^2}}\\ 1&{{v^2}}&{{v^2}}&{{v^2}}\\ 1&{{v^2}}&{{v^2}}&{{v^2}}\\ 1&{{v^2}}&{{v^2}}&{{v^2}} \end{array}} \right) * }\\ {\left( {\begin{array}{*{20}{c}} { - {b_{xx1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - {b_{xx2}}\frac{{{{\left( {{S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - {b_{xxN}}\frac{{{{\left( {{S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ { - {b_{yy1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - {b_{yy2}}\frac{{{{\left( {{S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - {b_{yyN}}\frac{{{{\left( {{S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ \vdots & \vdots & \vdots & \vdots \\ { - 2{b_{xy1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - 2{b_{xy2}}\frac{{{{\left( {{S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - 2{b_{xyN}}\frac{{{{\left( {{S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ { - 2{b_{xz1}}\frac{{S_1^2}}{{\sigma _1^2}}}&{ - 2{b_{xz2}}\frac{{{{\left( {{S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{ - 2{b_{xzN}}\frac{{{{\left( {{S_N}} \right)}^2}}}{{\sigma _N^2}}}\\ {\frac{{S_1^2}}{{\sigma _1^2}}}&{\frac{{{{\left( {{S_2}} \right)}^2}}}{{\sigma _2^2}}}& \cdots &{\frac{{{{\left( {{S_N}} \right)}^2}}}{{\sigma _N^2}}} \end{array}} \right)} \end{array} $

(26)

首先求解实对称矩阵的逆矩阵即求${{\left[ \begin{array}{*{35}{l}} 1 & \ldots & 0 & 0 \\ 0 & {{v}_{2}} & \ldots & 0 \\ {} & \ \ \ \ \ \ \vdots & {} & {} \\ 0 & 0 & 0 & {{v}_{2}} \\ \end{array} \right]}^{-1}}$, 由于该矩阵为对角矩阵因此该矩阵的逆矩阵为$\left[ \begin{array}{*{35}{l}} 1 & \ldots & 0 & 0 \\ 0 & \frac{1}{{{v}_{2}}} & \ldots & 0 \\ {} & \ \ \ \ \ \ \vdots & {} & {} \\ 0 & 0 & 0 & \frac{1}{{{v}_{2}}} \\ \end{array} \right]$

再次使用三重积性质可得

$ {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{{\mathit{\Sigma '}}^{ - 1}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{{\mathit{\Sigma '}}^{ - 1}}} \right) = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}} \right) $

(27)

因此

$ \begin{array}{*{20}{c}} {\alpha ' - \alpha = {{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}} \right)x' - }\\ {{{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\Sigma }^{ - 1}}} \right)\mathit{\boldsymbol{x}}} \end{array} $

(28)

又因为

$ \mathit{\boldsymbol{x'}} = \left( {\begin{array}{*{20}{c}} {\ln \left( {{S_1}} \right)}\\ {\ln \left( {v \times {S_2}} \right)}\\ \vdots \\ {\ln \left( {v \times {S_N}} \right)} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\ln \left( {{S_1}} \right)}\\ {\ln \left( {{S_2}} \right)}\\ \vdots \\ {\ln \left( {{S_N}} \right)} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {\ln 1}\\ {\ln v}\\ \vdots \\ {\ln v} \end{array}} \right) $

(29)

综上所述，当已知变换倍数时，可以求得原状态下的张量矩阵。