专利名称：：一种局部快速行进模式的弥散张量核磁共振图像配准方法
技术领域：
：本发明涉及一种局部快速行进模式的弥散张量核磁共振图像fffi准方法，属于脑科学，计算机视觉、图像理解以及模式识别等领域。适合于基于DTI图像的大脑结构与功能的研究，大脑神经性相关疾病的诊断与治疗，及临床神经外科手术的术前规划等领域。
背景技术：
：弥散张量核磁共振振图像(DTI)使得三维大脑白质神经纤维的可视化及量化分析成为可能。其在多发性硬化，中风，老年痴呆症，精神分裂K，癫痫，脑肿瘤等神经性相关疾病的研究中的成功应用，使得DTI技术成为近15年来脑研究方面最流行的核磁共振成像技术之一。在理论研究及实践应用中，为了比较从不同的对象或者同一个对象不同时间点获取的DTI图像，需要首先将这些图像配准到同一个参考空间中，所以体内和体外的DTI图像的空间配准是一项重要的任务。因为DTI图像的基本像素单位是张量，所以DTI图像的配准和普通的尺度图像的配准有很大的不同。以往，人们多利用从张量图像提出的一些尺度信息，如各向异性分数(FractionalAnisotropy，FA)、表观弥散系数(ApparentDiffusionCoefficient，ADC)及平均扩散系数(Meandiffiisivity，MD)等，或者只利用当前像素的张量信息，来比较像素之间的相似性。此类方法的缺点在于，前者只利用了部分张量信息，没有充分利用张量图像所提供的大脑白质的结构、方向信息；后者虽然利用了张量信息，却又只是基于当前像素来计算的，没有考虑邻域的相关性，也在很大程度上缺失了DTI图像所提供的最为重要的关于白质神经纤维走向的信息。针对以上方法的不足，在此发明中，我们采用一种基于局部快速行进(FastMarching,FM)的算法来提取像素周边的形态特征，即在整个大脑中，对每个感兴趣像素计算其局部FM时间图，作为张量图像配准的特征。相对于前面提到的特征，这个基于局部FM算法的弥散时间图不仅考虑到张量方向，还利用每一个像素的相邻张量信息作为独特的属性，从而反映了更加完整和全面的图像特征，可以更好的作为相似性度量特征来指导图像的配准。
发明内容为了避免现有技术的不足之处，本发明提出一种^于局部快速行进模式的弥散张量核磁共振图像配准方法，在参考快速行进(FastMarc1ing，FM)算法的基础上提出的DTI图像的配准方法。一种基于局部快速行进模式的弥散张量图像配准方法，其特征在于步骤如下步骤l:以分辨率/对DTI图像进行下采样，所述的分辨率为W-5;步骤2:当采样像素的部分各向异性值FA大于阈值时，利用快速行进模式P'astMarching-FM的方法模拟水分子在其邻域内的传播模式，计算该采样像素的局部FM时间图，具体步骤为步骤a:以行进速度^(1")=《，10+7《rt'a(r)自点y'推进到周l邻域其他候选点y，其中《，(0为主要行进速度、《^00为平滑行进速度、/7为主要行进速度和平滑行进速度之间的折中系数；所述的vT，卜w,ZVr，w是消除在脑脊液内部快速行进的系数，「y'—y，,为r的转置，Dy'为点y'的弥散张量；所述的w=1/(l+exp(-".(W-仰)，其中"和p是调节FA对快速行进进行控制的常量；所述的v广a(r)-Vy.(i;,.r)，其中!y和vy分别是当前点y'的冇进方向和速度，"."表示两个向量点积；步骤b:计算y'外部的每一个相邻点y的到达时间,y^y.+""y，其中，,/是到达y'的时间，Vy是主要行进速度；步骤c:以y'外部中最小到达时间的那个相邻点y作为新的前向点，然后从步骤a开始重复进行得到从像素x向前演化的局部FM时间图；重复的范围为以像素x点为球心的，半径为16mm的一个圆球内；步骤3:利用得到的局部FM时间图和部分各向异性值FA图作特征，建立DTI模板图像z^和目标图像z^之间变形场f的能量方程其中^,y)是模板图像r里面从像素x到像素y的行进时间，^f(x)+x，f(y)+y)是考虑变形场f后目标图S里面像素f(x)+x到像素f(y)+y的行进时间，l丄力是沿着行进轨迹的从x到y的距离，WW代表像素x的邻域，/^W是模板图像r里面像素x的部分各项异性值，/^.(f(x)+x)是考虑变形场f后目标图S里面像素f(x)+x的部分各项异性值，//是折中系数；步骤4:变形场f采用三次B样条配准模型来建模，对每一个三次B样条的控制点c，估计能量方程对于控制点的梯度e^，用v(c)-v(c)^%更新控制点的变形值，并更新受这一控制点影响的变形场；对所有的控制点进行迭代，以更新变形场；步骤5:重复步骤4直到收敛；完成分辨率〃下配准之后，如果r:则转到第6步；否则进行三次B样条上采样，在一个更高分辨率上计算它的相关变形场，在获得变形场之后，采用Xu的重定向方法对张量图像进行重定向操作，并设置r〃-l，然后转到第2步进行下一个分辨率处理；+//(i^(f(x)+x)-叫(x))2！.步骤6:用产生的三次B样条模型计算最终的变形场；并釆.目Xu的重定向方法重定向目标图像到模板空间上。所述的部分各向异性值FA大于阈值为0.25。所述的7选择0.9。所述的"选择50。所述的"选择0.3。所述的A选择l.O。本发明提出的基于局部快速行进模式的弥散张量核磁井振图像配准方法，在配准过程中不仅仅考虑到张量方向，还利用每一个像素相邻的张J3"信息作为独特的属性；同时采用基于三次样条的配准方法来计算变形场。在具体实现过程中，针对当变形场变化之后，需要重新计算时间图，计算量太大的特点，本发明还提出了一种时间图更新的快速算法，可以在保证配准精度的情况下，大大的提高计算效率。大量的实验证明，相对于其它方法，本发明提出的基于FM模式的配准方法可以提取更丰富、更有效的张量特征，从而更好地界定图像之间点的对应关系，可以产生更稳定、更精确的配准结果。图l:基于张量特征和基于FM时间图的相似性度量的比较(a)模板图(FA)及一个参考点；(b)目标图(FA);(c)利用张量特征计算得到的距离图；(d)利用FM特征计算得到的距离图2:局部快速行进(FM)时间图(a)局部FM时间图；(b)(c)和(d)说明在不同时间的演变行进路线；图3:简化算法的示例说明图4:模拟图像实例(a)模板图；(b)(C)(d)模拟图像；图5:配准误差分布图图6:白质神经纤维束标记实例(a)模板图像中的白质神经纤维；(b)在配准目标图像中，提取的白质神经纤维；(C)在原始目标图像中手工提取的白质神经纤维。具体实施例方式现结合附图对本发明做进一步描述(1)问题的描述本文提出的DTI图像配准算法与普通的尺度图像配准算法很类似，目标是求出一个变形场f，以达到配准DTI模板图像^和目标图像^的目的。能量方程定义如下￡(f)=￡，(",.,Ds,f)"五咖(f)，(1)其中，fi、,JZV,A，f)是根据当前变形场f计算的两幅图像之间的相似性度量，￡f。(f)是变形场的光滑约束。在本发明中，变形场f采用三次B样条配准模型来建模，可以保证变形场的光滑，因此可以忽略第二项。从而，详细的能量函数￡=￡,,重新定义为￡=4Z)s(f(X)+X)]—/f[Q(f(X))￡V(x)27(f(X))]||2(2)其中x是模板图像区域Q内的一个像素，Z),W是在模板图像中像素x处的张量，A、-(fOO+x)是目标图像中相对应像素的张量，2(f(x))是在像素x处根据变形场f计算出来的旋转矩阵。与普通尺度图像配准方法不同的是，在张量图像的配准过程中，为了与目标图像相对应的张量进行比较，模板图像的张量需要重定向，e(fw)正是为张量重定向而计算得到的旋转矩阵。4]表达了以一个像素为中心的弥散张量的特征提取操作，且M[A]-4A]p给出了一组特征向量的距离。如果4]=/，一个像素的图像特征就是张量本身，并且因此M[A]-4A]II1A-A1卩变成一种简单的张量距离。此种情况下，公式2所表达的含义就和普通的基于张量相似性的配准算法完全一样了。在本发明中，邻域张量特征4]定义为通过执行局部FM算法而得到时间图。虽然单个像素领域内的张量是一种高维系统，但通过一种简单的基于FM的算法，我们可以有效的提取像素周围的张量模式。因此，本发明的主要想法，就是利用局部Bf算法所得到的更有效的张量特征，来提高DTI图像配准的精度。除了以局部FM为基础的张量特征，FA值也作为一个特征向量元素以更多强调当前张量点。为了阐明为什么以FM为基础的张量特征有更好的可区分性，图l分别展示了以简单张量距离和基于FM张量特征计算的距离图。此处，FM算法在以每一个像素为球心的，半径为16mm的一个圆球内进行。可以看出，如果仅使用张量距离，在整个模板图像中有大量的与选定的参考像素非常相似的点；而基于FM的张量特征则能更好的帮助区分参考像素与其他像素。因此，这种独特的特征可用于配准以取得更稳定的结果。(2)基于局部FM的弥散张量特征提取在本发明中，对每一感兴趣像素，我们运用局部FM算法模拟水分子在其邻域内传播模式，作为当前像素的张量特征用于指导图像配准。首先，演进前面从当前像素x开始向周围邻域进行演进，假设在某个时间点到达了点y'，然后从y'推进到其他周边候选点y，演进的速度可以定义为<formula>formulaseeoriginaldocumentpage9</formula>其中，、(r)是主要演进速度，f:y'—y，并且其中，w是在脑脊液内部消除快速演进的系数。FA值小则w小，反之亦然。因此w可以通过<formula>formulaseeoriginaldocumentpage9</formula>计算得至IJ，其中"和/是控制FA效果的常量。公式3的第二项用来保证演进的平滑<formula>formulaseeoriginaldocumentpage9</formula>其中r,和vy.分别是当前点y'的演进方向和速度，"."表两个向I:点积。FM算法中，演进前面的演进到达时间和演进速度的关系可以用Eikonal方程来描述<formula>formulaseeoriginaldocumentpage9</formula>其中，vr是从一个点到另一个点演进所需时间，「是它们之间的演进速度，"是这些点之间的距离，因此到达前向点y的时间,y是<formula>formulaseeoriginaldocumentpage10</formula>其中，V是到达y'的时间，vy是公式3所定义的演进速度。对于y'外部的每一个相邻点，我们利用公式7计算到达它时间，选择具有最小到达时间的点作为新的前向点。这样，从像素x向前演化，x邻域点的到达时间可以形成一个包含有周围张量信息特征的图，称之为局部FM时间图。在发明中，〃选择0.9，"选择50，"选择0.3。图2显示了一个给定像素的局部FM时间图，可以看出，时间图的分布反映了该像素周围张量分布信息，也反映了给定像素邻域内白质神经纤维f勺走向分布信息。在本文方程2中，局部FM时间图被用作是张量特征4]。详细的实施将在下一节描述。(3)基于局部FM时间图的配准算法根据以上的方法，对整个大脑的所有感兴趣像素都可以得到一个局部FM时K图。从图2可以看出弥散速度或者到达时间沿着不同的方向非常不同，这取决于DTI图像本身。事实上，沿着白质神经纤维的方向演进快而其他方向炉。因此，从基于张量FM得到的时间图反映了一个像素周围的神经纤维的存在模式。，'面介绍配准算法的实现回到公式2，用局部FM时间图来计算张量特征，因此公式2可以写成<formula>formulaseeoriginaldocumentpage10</formula>其中，《(x,y)是目标图S里面从像素x到像素y的演进时间，W)z^(f)，f(x,y)是考虑到变形场f及张量重定向，变形后模板图从像素x到像素y的时间，W(x)代表像素x的邻域。明显，一旦变形场f改变了，每一个模板像素周围的时间图相应的也必然会改变，所以新的时间图就必须在变形和重定向后模板图像的基础上重新计算，这是一个非常耗时的工作，因此我们对^(，^(f)，f(X，y)的实现做了简化，具体实现为、(f)D^r(f),f(x,y)纡(x,y)d+1")",)—:'二厂iy—x"),(9)其中，i."y是沿着行进轨迹的从x到y的距离，最终，简化的能量函数可以被定义为:qi:^(f(x)+x,f(y)+y)—,7.(x，y)(l+if(勿^:如『)xenL"ww、'(10)+/^s(f(x)+x)-F4(x))2}图3给出了这种简化的示例说明对光滑的弹性变形场而言，沿着FM轨迹的局部角变化相对很小。比如，从左下角红色张量演进至右上角张量的FM轨迹，在光滑弹性变形前后，几乎一样。在本发明的方法中，假定沿着这些演进轨迹的角变化很小，因此，时间图的改变可以近似的根据变形前后距离的变化来计算。(4)综合以上，基于局部快速演进模式的弥散张量图像配准方法可以总结如下步骤l，选择分辨率R，对输入的DTI图像进行下采样，用丁多分辨率配准；并设置当前分辨率^W;步骤2，在分辨率r下，对每一个像素，如果它的FA值大于给定阈值，计算它的局部FM时间步骤3，对每一个三次B样条的控制点c，估计目标函数对于控制点的梯度，如%，用*)=一-￡%更新控制点的变形值，并更新受这一控制点影响的变形场；对所有的控制点进行迭代，以更新变形场。步骤4，重复步骤3直到收敛。完成分辨率r下配准之后，如果n则转到步骤5，否则进行三次B样条上采样，在一个更高分辨率上计算它的相关变形场。在获得变形11场之后，采用Xu的重定向方法对张量图像进行重定向操作。设置"rl，然后转到步骤2进行下一个分辨率处理。步骤5，用产生的三次B样条模型计算最终的变形场；如果必要，则重定向目标图像到模板空间上。注意，在多分辨率的实现中，为了消除由于简化算法而产生的误差，我们在每一个分辨率下都对变形模板图像重新计算局部FM时间图。(5)配准结果及分析为了评估本发明提出的配准算法，进行了两组实验。在第一组实验中，用已知的变形场来模拟数据，然后通过将己知变形场和本发明提出的算法得到的变形场进行比较，来评估算法的精确性。第二个实验中，在真正的弥散张量图像上，我们应用本算法来进行配准，以检测它的性能。A实验一基于仿真数据的实验首先，应用Xue的基于统计模型的弹性仿真算法来模拟IO个DTI图像。具体方法是先通过Xue的方法模拟10个弹性变形场，然后用Xu的DTI图像配准算法将一己知的模板图像，通过前面模拟的IO个变形场配准到目标空间中。这样，我们就可以模拟得的到10个DTI图像，并且从模板图像到这些模拟图像的弹性变形场是己知的。同时计算各向异性分数图像(FA)，作为此算法的一维特征用于图像配准中。图四显示了模板图像及三个随机选取的模拟图像的FA图，可以看到这些通过模拟得到的FA图像显著的不同于模板图像。因为这些模拟图像和原始的模板图像之间的夷实变形场是已知的，通过比较真实变形场与利用此发明所提出算法得到的模板图像和真实图像之间的变形场，可以进行算法的评价。作为比较，应用传统的DTI图像配准算法，即利用高维自由形式的配准算法配准尺度FA图像，来得到作为算法比较的变形场。完成模板图像和仿真图像的配准之后，计算配准变形场和真实变形场之间的变形误差分布。假定真实变形场为f，计算的变形场为f，两者的差值图像M为f和f'之间每一像素点的欧式距离。图5显示了应用本发明提出的配准算法及传统配准算法于10个仿真图像及模板图像之后所得到平均变形误差分布图。可见，应用了邻域张量信息的本发明算法具有更小的变形误差和更精确的配准结果。B实验二真实的DTI图像实验第二组实验的目的是检查本文所提算法在真实DTI图像上的表现。为了能够可视化评估算法的表现，我们用所提算法把14个正常成年人的DTI图像自动配准到一个模板图像上，从模板图像中手动选定某个具体纤维束作为参考纤维束(记为纤维束A)，然后根据Hausdorff距离在配准的目标图像中选取与纤维束A距离相近的纤维束作为目标图像中的感兴趣纤维束(记为纤维束B)。同时手工的在原始目标图像空间中选取与纤维束A同一类的纤维束(记为纤维束C)。比较纤维束B和纤维束C，可以对此发明提出的配准算法做基于视觉的验证。图6显示以胼胝体为例，对该算法进行验证的结果。图6(a)显示从模板图像手动选择的胼胝体纤维束；图6(b)显示了在配准的目标图像自动标记的感兴趣纤维；图6(c)显示了从原始目标图像中手动标记的胼胝体纤维束。可以看出，手动选择和自动计算的纤维束非常相似，从而证明本发明提出的配准算法具有较好的配准精度，并且利用该算法可以实现某些神经纤维束的自动提取。权利要求1.一种基于局部快速行进模式的弥散张量图像配准方法，其特征在于步骤如下步骤1以分辨率R对DTI图像进行下采样，所述的分辨率为R＝5；步骤2当采样像素的部分各向异性值FA大于阈值时，利用快速行进模式FastMarching-FM的方法模拟水分子在其邻域内的传播模式，计算该采样像素的局部FM时间图，具体步骤为步骤a以行进速度id="icf0001"file="A2009100238490002C1.tif"wi="52"he="6"top="76"left="80"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>自点y′推进到周围邻域其他候选点y，其中vytensor(r)为主要行进速度、vyinertia(r)为平滑行进速度、η为主要行进速度和平滑行进速度之间的折中系数；所述的id="icf0002"file="A2009100238490002C2.tif"wi="28"he="4"top="107"left="55"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>w是消除在脑脊液内部快速行进的系数，ry′→y，rT为r的转置，Dy′为点y′的弥散张量；所述的w＝1/(1+exp(-α·(FA-β)))，其中α和β是调节FA对快速行进进行控制的常量；所述的id="icf0003"file="A2009100238490002C3.tif"wi="33"he="4"top="146"left="55"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>其中ry′和vy′分别是当前点y′的行进方向和速度，“.”表示两个向量点积；步骤b计算y′外部的每一个相邻点y的到达时间ty＝ty′+d/vy，其中，ty′是到达y′的时间，vy是主要行进速度；步骤c以y′外部中最小到达时间的那个相邻点y作为新的前向点，然后从步骤a开始重复进行得到从像素x向前演化的局部FM时间图；重复的范围为以像素x点为球心的，半径为16mm的一个圆球内；步骤3利用得到的局部FM时间图和部分各向异性值FA图作特征，建立DTI模板图像DT和目标图像DS之间变形场f的能量方程<mathsid="math0001"num="0001"><math><![CDATA[<mrow><mi>E</mi><mrow><mo>(</mo><mi>f</mi><mo>)</mo></mrow><mo>=</mo><mfencedopen=''close=''><mtable><mtr><mtd><munder><mi>Σ</mi><mrow><mi>x</mi><mo>∈</mo><mi>Ω</mi></mrow></munder><mo>{</mo><munder><mi>Σ</mi><mrow><mi>y</mi><mo>∈</mo><mi>N</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></munder><msup><mrow><mo>(</mo><msub><mi>t</mi><mi>S</mi></msub><mrow><mrow><mo>(</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>+</mo><mi>x</mi><mo>,</mo><mi>f</mi><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>+</mo><mi>y</mi><mo>)</mo><mo>-</mo><msub><mi>t</mi><mi>T</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mrow><mo>(</mo><mn>1</mn><mo>+</mo><mfrac><mrow><msub><mrow><mo>|</mo><mi>f</mi><mrow><mo>(</mo><mi>y</mi><mo>)</mo></mrow><mo>+</mo><mi>y</mi><mo>-</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>-</mo><mi>x</mi><mo>|</mo></mrow><mrow><mi>x</mi><mo>→</mo><mi>y</mi></mrow></msub><mo>-</mo><msub><mrow><mo>|</mo><mi>y</mi><mo>-</mo><mi>x</mi><mo>|</mo></mrow><mrow><mi>x</mi><mo>→</mo><mi>y</mi></mrow></msub></mrow><msub><mrow><mo>|</mo><mi>y</mi><mo>-</mo><mi>x</mi><mo>|</mo></mrow><mrow><mi>x</mi><mo>→</mo><mi>y</mi></mrow></msub></mfrac><mo>)</mo></mrow><mo>)</mo></mrow></mrow></mrow><mn>2</mn></msup></mtd></mtr><mtr><mtd><mo>+</mo><mi>μ</mi><msup><mrow><mo>(</mo><msub><mi>FA</mi><mi>S</mi></msub><mrow><mo>(</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>+</mo><mi>x</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>FA</mi><mi>T</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>}</mo><mo>.</mo></mtd></mtr></mtable></mfenced></mrow>]]></math></maths>其中tT(x，y)是模板图像T里面从像素x到像素y的行进时间，tS(f(x)+x，f(y)+y)是考虑变形场f后目标图s里面像素f(x)+x到像素f(y)+y的行进时间，|·|→y是沿着行进轨迹的从x到y的距离，N(x)代表像素x的邻域，FAT(x)是模板图像T里面像素x的部分各项异性值，FAS(f(x)+x)是考虑变形场f后目标图s里面像素f(x)+x的部分各项异性值，μ是折中系数；步骤4变形场f采用三次B样条配准模型来建模，对每一个三次B样条的控制点c，估计能量方程对于控制点的梯度id="icf0005"file="A2009100238490003C1.tif"wi="5"he="4"top="91"left="92"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>，用id="icf0006"file="A2009100238490003C2.tif"wi="26"he="4"top="91"left="108"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>更新控制点的变形值，并更新受这一控制点影响的变形场；对所有的控制点进行迭代，以更新变形场；步骤5重复步骤4直到收敛；完成分辨率r下配准之后，如果r＝1则转到第6步；否则进行三次B样条上采样，在一个更高分辨率上计算它的相关变形场，在获得变形场之后，采用Xu的重定向方法对张量图像进行重定向操作，并设置r＝r-1，然后转到第2步进行下一个分辨率处理；步骤6用产生的三次B样条模型计算最终的变形场；并采用Xu的重定向方法重定向目标图像到模板空间上。2.根据权利要求1所述的基于局部快速行进模式的弥散张量图像配准方法，其特征在于所述的部分各向异性值FA大于阈值为0.25。3.根据权利要求1所述的基于局部快速行进模式的弥散张量图^配准方法，其特征在于所述的7选择0.9。4.根据权利要求1所述的基于局部快速行进模式的弥散张量图像配准方法，其特征在于所述的"选择50。5.根据权利要求1所述的基于局部快速行进模式的弥散张量图像配准方法，其特征在于所述的"选择0.3。6.根据权利要求1所述的基于局部快速行进模式的弥散张量图像配准方法，其特征在于所述的/z选择1.0。全文摘要本发明涉及一种基于局部快速行进模式的弥散张量核磁共振图像配准方法，技术特征在于首先采用局部FM算法来提取弥散张量图像中某像素周边的形态特征，得到每个感兴趣像素的局部FM时间图；然后采用三次B样条配准模型建立变形场，对局部FM时间特征图进行校准；在配准时同时考虑了模板图的变形和重新定向，用修正后的局部FM算法得到的模板时间特征图和目标时间特征图的误差来修正变形场，进行目标图与模板图的配准。相对于其它方法，本发明提出的基于FM模式的配准方法可以提取更丰富、更有效的张量特征，从而更好地界定图像之间点的对应关系，可以产生更稳定、更精确的配准结果。文档编号G01R33/54GK101666865SQ200910023849公开日2010年3月10日申请日期2009年9月9日优先权日2009年9月9日发明者张德刚,海李,忠薛,雷郭申请人:西北工业大学