2019, Vol. 37
2. 中国空气动力研究与发展中心, 绵阳 621000
2. China Aerodynamics Research and Development Center, Mianyang 621000, China
模型飞行试验[1]作为空气动力学研究三大手段之一，其独特优势和重要价值历来都很受各航空航天发达国家的重视[2-3]。MF-1[4-5]是中国空气动力学研究与发展中心开展的单级固体火箭助推不分离无控高超声速空气动力学飞行试验研究项目，主要目的是通过在试验模型特定部位安装的温度、压力等传感器，对飞行试验全过程参数变化历程实现测量，为边界层转捩、激波边界层干扰机理研究提供真实飞行数据。MF-1试验飞行器全长6.23 m，最大直径0.6 m，由试验模型、级间段、助推器和尾段等组成(图 1)。
|
图 1 MF-1试验飞行器外形图 Fig.1 Sketch of MF-1 flight vehivle |
MF-1飞行试验数据表明，其真实飞行弹道与发射前设计弹道有明显差异，落点纵向射程偏近约34 km，偏右约25 km，弹道顶点高度偏低约15 km。本文首先对可能导致弹道偏差的原因进行了分析，并结合弹道仿真给出了导致该偏差的因素的量级。然后，基于飞行弹道重建技术[6-7]，利用遥测过载、角速率和外测弹道等数据，重建了包括飞行迎角、侧滑角等在内的全弹道数据，为MF-1飞行试验后的边界层转捩、激波边界层干扰机理研究提供支撑。
1 MF-1飞行弹道偏差现象及分析 1.1 弹道设计流程及实际飞行偏差现象根据飞行试验流程安排，在当日11：00放飞气球进行气象测量；13：00前提供气象测量结果给弹道设计岗位用于风修弹道设计；13：30前将弹道设计得到的发射架仰角和方位角提供给发射架岗位进行发射架角度瞄准；16：00MF-1试验飞行器准时发射。
基于11：00的气象测量数据、发动机10 ℃药温内弹道理论值，设计得到MF-1发射方位角269.82°、发射倾角62.02°。图 2给出了设计弹道与飞行弹道全程y-x曲线、x-z曲线及其起飞段放大曲线对比，x、y、z为MF-1在发射坐标系中的位置坐标，其中y-x起飞段放大曲线中，发射点附近飞行弹道的不平滑是由于发动机点火对GPS接收信号的影响造成的。图 3给出了起飞段弹道倾角θ及弹道偏角σ对比。
|
图 2 MF-1飞行弹道与设计弹道对比 Fig.2 Comparison of flight trajectory to designed trajectory |
|
图 3 飞行弹道角与设计弹道角对比 Fig.3 Path angle comparison of flight trajectory to designed trajectory |
图中，x坐标为0.5 km时对应的飞行时间约3.9 s，x坐标为5 km时对应的飞行时间约12.8 s。从图可见，在MF-1飞行的开始阶段，实际飞行弹道就明显比设计弹道倾角偏低，地面轨迹偏右，弹道倾角和弹道偏角均在刚刚起飞就迅速产生显著偏差。因此，可以断定，对于无控飞行的MF-1来说，其落点位置产生较大偏差的最主要原因就是由于弹道起飞初期就存在的偏差导致的。
1.2 飞行弹道偏差原因初步分析在MF-1发射起飞初期，弹道角在很短时间内出现偏离，一种可能是发射架的角度在瞄准调节时存在较大偏差；另一种可能是，发射离架后，其他因素使得飞行弹道迅速出现了偏离。
理论上，对于静稳定性很高、配平飞行迎角在0°附近的无控飞行器，在无风情况下，其气动力偏差和发动机推力偏差主要影响加速性能，不会迅速引起弹道角变化。在MF-1发射离架初期，由于飞行速度低、气动力小，气动不确定性偏差更是不可能迅速对飞行弹道产生显著影响。发射离架初期的固体助推火箭推力很大，加速度大于10 g，在无其它偏差条件下，即使推力存在一定的偏差，也不会在极短时间内导致弹道角出现显著偏离。但是，如果地面附近风场存在偏差，由于离架初期速度小，可能出现明显的风迎角。比如，相对于30 m/s的飞行速度，若在垂直于飞行速度方向存在风速为3 m/s的风，风迎角可能达到5.7°。由于大的气动静稳定性，飞行器的姿态会迅速改变以消除风迎角。在强大的发动机推力作用下，飞行器沿着改变姿态后的方向(不同于飞行速度方向)迅速加速，飞行速度方向即弹道角也就很快与之前发生了偏离。
因此，通过理论分析可知，导致MF-1飞行弹道在初始段就产生明显偏离的最主要潜在原因有两个，其一是发射架角度偏差，其二是风场偏差。
下面对风场偏差、发动机推力偏差、发射架角度偏差、气动力偏差等因素影响分别进行仿真，以探求多大的偏差量级，才可能导致落点位置与设计弹道之间产生那么大的偏差。基本仿真方法与弹道设计相同，参见文献[4]。
2 MF-1飞行弹道偏差影响仿真 2.1 风场偏差影响仿真MF-1飞行试验发射架的角度调整是基于发射前5 h放飞的气象气球测量结果进行弹道设计得到的，而且发射架位置距离气象气球放飞地点约30 km。如此大的时间和空间距离，气象测量必然不可能准确描述MF-1发射和飞行时所处的真实风场，特别是对弹道影响最大的刚离架时地面附近的风场。
在MF-1发射的同时，气象气球放飞点又放飞了一个气象气球进行气象测量。图 4给出了两次气象测量得到的风速Vwind和风向φwind对比图及地面附近局部放大图。
|
图 4 两次气象测结果对比 Fig.4 Comparison of the two measured atmospheric parameters |
从图 4可见，两次测量的高空风速风向具有较好的一致性，这也说明高空大气的流动由于不太受地形地貌影响平稳性较好，可提供给一定的时间空间范围内的飞行试验使用。然而在低空，特别是距离地表500 m高度以下，两次测量得到的风向差异显著，风向相差几乎180°，基本完全相反。500 m高度以下，两次测量的风速大小相差约2~6 m/s，按照4 m/s的风速计，相当于用于弹道设计的地面风速和MF-1实际飞行偏差高达8 m/s，发射离架初期的风迎角偏差可高达10°~15°。况且，地表附近的风场受地形地貌和日照辐射等影响比较大，随机性也较大，距离发射场30 km的地点其地面风场必然与发射当地风场存在差别。
发射架角度采用MF-1发射实际使用的角度，即根据11：00气象测量数据设计得到的角度，保持其它所有参数不变，仅仅将气象参数替换为16：00的测量数据，重新对MF-1发射飞行全弹道进行积分仿真，并与设计弹道和实际飞行弹道进行对比，图 5给出了对比结果。从图可见，基于11：00风场数据设计得到的发射角，采用16：00风场数据积分得到的弹道，其星下点轨迹明显比更设计弹道更靠近于飞行弹道，但仍然与实际飞行弹道有显著差异。
|
图 5 采用16:00风场数据的积分弹道与设计和飞行弹道对比 Fig.5 Comparison of the integrated trajectory with the use of 16:00 measured wind data to the designed and flight trajectory |
保持发射角不变，在采用11：00风场数据的基础上，经仿真，再叠加上速度为11.8 m/s、风向为西偏北35°的常值风，就可以使弹道积分落点与MF-1实际落点一致。图 6给出了此条件下的积分弹道与设计弹道和飞行弹道对比，从图可见，虽然积分弹道的落点和星下点轨迹与飞行弹道吻合，但积分弹道的高度轨迹明显比飞行弹道更低，这表明，虽然弹道装订采用的风场与发射时刻真实风场之间的偏差可能是导致飞行弹道偏差的重要因素，但是飞行弹道偏差并不完全是由于风场偏差导致的。在射程远大于射高情况下，同样的落点，更低的弹道必然对应更高的飞行速度和更短的飞行时间，图 7给出的这三条弹道V-x坐标对比曲线和V-T对比曲线也证明了这一点。
|
图 6 叠加常值风的积分弹道与设计和飞行弹道对比 Fig.6 Integrated trajectory by adding constant wind data compared with the designed and flight trajectory |
|
图 7 三条弹道速度历程对比 Fig.7 Velocity history comparison of the three trajectories |
如前所述，理论分析可知，导致MF-1在发射起飞初期就产生明显的弹道角偏差的原因有可能是因为发射架的角度调节不准确，初始发射角存在明显偏差导致的。当然，也因为初始发射架偏差对发射初段弹道偏差的影响与风场偏差类似，可以预计，即使通过调整初始发射架使得弹道积分落点与飞行吻合，积分弹道也会出现比实际飞行弹道更低、速度更高的现象。
经仿真，在采用11：00风场数据的基础上，将发射架的发射俯仰角设计值向小调整约5.6°，方位角向右调整约4.3°，弹道积分结果的落点就可以与实际飞行落点一致。图 8给出了此条件下的积分弹道与设计弹道和飞行弹道对比，从图可见，虽然积分弹道的落点和星下点轨迹与飞行弹道吻合，但确实与风影响类似，积分弹道的高度轨迹明显比飞行弹道更低。而且，根据发射架研制单位提供的产品技术指标，发射架仰角和方位角调节精度很高，分别为±0.05°和±0.1°。本节中4°~6°的发射架角度偏差显然偏离其指标参数太大，且在该发射架执行过的其它发射任务中，也没有证据能够证明发射架角度调节会有这么大的偏差。
|
图 8 调整发射角的积分弹道与设计和飞行弹道对比 Fig.8 Integrated trajectory by modifying initial launch angle compared with the designed and flight trajectory |
阻力偏差也可能造成航程偏差。由于MF-1静稳定度很高，飞行迎角很小，阻力近似等于轴向力，因此本节的仿真是通过改变轴向力系数的大小使得仿真的落点纵向航程与飞行一致，发射角仍然采用设计角度。对于侧向落点的调整，则是通过增加来自于正北方的侧向风实现的。经仿真，在采用11：00风场数据的基础上，将MF-1轴向力系数增大为1.125倍，再叠加5.8m/s的北风，积分仿真弹道落点将与飞行落点一致。图 9给出了此条件下的积分弹道与设计弹道和飞行弹道的飞行轨迹及速度对比。从图可见，虽然积分弹道的落点和星下点轨迹与飞行弹道吻合，正如1.1节所做的分析，由于轴向力系数变化不会改变发射起飞初期的姿态，因此积分弹道比实际飞行弹道更高，全程基本都介于设计弹道和实际飞行弹道之间。而且，由于要以较高的弹道实现与飞行弹道一致的落点，积分弹道的速度比飞行弹道更低。
|
图 9 改变轴向力系数的积分弹道与设计和飞行弹道对比 Fig.9 Integrated trajectory by modifying axial coefficient compared with the designed and flight trajectory |
MF-1飞行过程中对发动机燃烧室压力进行了测量。图 10给出了燃烧室压力飞行测量结果与0 ℃、5 ℃和10 ℃药温条件下理论压力对比。从图可见，从时间轴来说，MF-1发动机实测的燃烧室压力与5 ℃药温内弹道一致性较好，但几乎整个燃烧过程的推力都偏小，4.5 s前的大推力段实际飞行值偏低5%~7%。如果该数据可信，则必然意味着飞行过程的推力和总冲明显小于设计状态，也必将导致飞行速度和射程比设计弹道偏小。
|
图 10 燃烧室压力飞行值与不同药温理论值对比 Fig.10 Combustor pressure flight data compared with that of theoretic data under different fuel temperatures |
MF-1采用的助推器是成熟批产并进行了大量飞行的货架产品，对该助推器产品的技术参数分析发现，该系列助推器装药量的设计值最大偏差要求不大于0.6%；MF-1采用的这发助推器，实测装药量偏差小于0.2%。从该数据来看，助推器实际飞行推力与理论推力似乎不应该有显著差异。
图 11(a)给出了采用设计发射角、11：00风场和依据燃烧室压力换算得到的发动机实测推力的积分弹道与设计弹道和飞行弹道对比。可见相比于图 5，采用发动机实测推力的积分弹道射程显著向飞行弹道靠拢，但仍与实际弹道有明显差别。在此基础上，再将发动机推力进一步减小0.95%，并叠加上5.67 m/s的正北风，可以实现积分仿真落点与飞行落点一致，如图 11(b)所示，但由于发动机推力变化不会迅速改变发射起飞初期的弹道角，仿真积分弹道的高度仍然比飞行弹道更高、速度更低。
|
图 11 发动机实测推力的积分弹道与设计弹道和飞行弹道对比 Fig.11 Comparison of the integrated trajectory with the use of flight-measured thrust data to the designed and flight trajectory |
前面的仿真表明，通过改变发射架方位角或者叠加侧向风，必然可以使仿真积分弹道的星下点轨迹与飞行一致；通过改变发射架俯仰角或者叠加纵向风，都可以在MF-1发射起飞初期就可以迅速改变弹道倾角；MF-1气动力偏差和发动机推力偏差不会引起发射起飞初期的弹道角迅速改变，但可以对射程产生影响。因此，通过将能够改变初始弹道角的因素与不能改变弹道角的因素组合，理论上能够得到与飞行轨迹基本一致的仿真弹道。
组合条件1：发射架角度采用设计值；风场采用11：00测量数据再叠加速度为3.9m/s、西偏北62°常值风，轴向力系数增加8.7%，发动机推力数据采用5 ℃药温理论值。
组合条件2：发射架仰角相对设计值减小1.8°，射向朝右调整3.9°，风场采用11：00测量数据，轴向力系数增加9.3%，发动机推力数据采用5 ℃药温理论值。
组合条件3：发射架角度采用设计值；风场采用11：00测量数据再叠加速度为6.4m/s、西偏北68°常值风，轴向力系数不变，发动机推力数据采用实测推力。
组合条件4：发射架仰角相对设计值减小1.55°，射向朝右调整3.9°，风场采用11：00测量数据，轴向力系数不变，发动机推力数据采用实测推力。
图 12~图 15分别给出了这4种情况下仿真弹道与飞行弹道的纵向射面轨迹对比及飞行时间50 s前的轴向总过载对比。从图中的纵向射面轨迹曲线可见，4种组合条件下的仿真轨迹与飞行轨迹一致性都比较好。这也说明，发动机推力采用理论推力、但将地面预测轴向气动力增加9%左右对全弹道高度和射程的影响，与采用地面预测轴向气动力、但将发动机推力从理论推力减小为实测推力对全弹道高度和射程的影响相当。进一步分析图中的过载对比曲线可见，采用理论推力的仿真弹道轴向总过载与飞行测量过载的吻合度明显比实测推力仿真弹道更好。而下节的弹道重建结果将表明，基于实测过载的重建弹道与飞行弹道吻合很好，证明了飞行实测过载的偏差可忽略。因此，相对来说，MF-1实际飞行的发动机推力更接近于5 ℃药温条件下的理论设计值。至于是什么原因导致了MF-1发动机燃烧室压力飞行测量值与理论值之间的差异，后续还需要发动机、测量等专业人员进一步通过分析计算或试验开展研究。
|
图 12 组合条件1仿真弹道与飞行弹道对比 Fig.12 Integrated trajectory of condition No.1 compared with flight trajectory |
|
图 13 组合条件2仿真弹道与飞行弹道对比 Fig.13 Integrated trajectory of condition No.2 compared with flight trajectory |
|
图 14 组合条件3仿真弹道与飞行弹道对比 Fig.14 Integrated trajectory of condition No.3 compared with flight trajectory |
|
图 15 组合条件4仿真弹道与飞行弹道对比 Fig.15 Integrated trajectory of condition No.4 compared with flight trajectory |
因此，综合本节的影响仿真和下节的弹道重建，我们认为：MF-1实际飞行时的发动机推力与理论设计推力之间不存在显著差异；引起飞行弹道与设计弹道偏差的最主要因素可能包括风场偏差、气动力偏差及发射架角度偏差。
3 MF-1飞行弹道重建弹道重建，即利用所有可以利用的飞行试验测量数据计算飞行状态向量[8]，其最初的目的是重构某些在动态飞行条件下难以直接测量的变量，如速度、迎角等[9]。其后不久，弹道重建方法被扩展应用于辨识一些传感器的系统误差[10]。今天，弹道重建作为检验传感器精度和飞行试验数据相容性的手段，在飞行器飞行试验数据分析中得到广泛应用，在某种程度上已经成为气动参数辨识不可或缺的重要环节[11]。
弹道重建不需要对气动力/力矩、发动机推力/力矩以及大气特性数据进行计算，而是直接利用MF-1发射初始状态参数和飞行实测的过载及角速率数据进行积分来重建整个飞行过程的弹道参数。因为气动力、推力、风场对于飞行器弹道和姿态的影响都已体现在过载和角速率中，因此气动力、推力、风场等偏差不会对弹道重建产生影响。但是由于飞行测量数据也不可避免地存在各种误差，导致直接积分出的弹道往往与真实弹道不一致。
应用系统辨识理论[11]，以六自由度动力学模型和误差模型为系统辨识的状态方程和观测方程，对MF-1整个飞行弹道进行积分，目标函数取为位置、速度等弹道参数拟合误差最小，待辨识参数为积分弹道初值误差及角速率和过载测量数据的常值漂移误差。利用辨识结果对弹道初值及角速率和过载测量数据进行修正后，再进行弹道积分并考虑风修正，即可重建包括迎角、侧滑角等在内的飞行试验全弹道数据。具体的弹道重建方法见文献[6, 7]。
图 16给出了MF-1飞行试验过载和角速率测量数据。从图可见，过载测量精度差，分辨率约0.24 g。如此低的分辨率必然导致难以进行高精度气动参数辨识。
|
图 16 MF-1飞行试验过载和角速率测量数据 Fig.16 Flight tested overloads and angular velocities of MF-1 |
图 17给出了直接采用设计的发射初始条件和飞行实测的过载角速率数据进行六自由度积分得到的全程弹道与外测弹道对比。从图可见，在纵向射面内积分弹道总体上与外测弹道一致性较好。这也说明过载、角速率测量数据及MF-1的实际发射架仰角都没有大问题；但是积分弹道星下点轨迹与外测弹道相比，明显向左侧偏斜，而且在落点附近出现了迅速拐弯的偏差发散情况，说明用于弹道积分的发射方位角等数据可能还是存在一些偏差，需要进行弹道重建。
|
图 17 基于飞行试验测量数据的MF-1积分弹道与外测弹道对比 Fig.17 Integrated trajectory based on flight test data compared with flight trajectory |
采用弹道重建技术，对积分初值偏差及过载角速率漂移进行辨识，并对辨识结果进行显著性分析，得到滚转角速率漂移为-0.013°/s、发射架仰角相对于设计值偏小0.1°、发射方位角相对于设计值偏右1.04°，其它偏差影响不显著。在此基础上，进一步进行了仿真分析，表明积分弹道星下点在末段出现的拐弯发散主要由滚转角速率漂移引起积分误差累积所致，拐弯前的弹道偏差主要由初始发射角偏差引起。
采用修正后的发射初始状态及过载角速率数据，进行全弹道积分获得了MF-1全弹道重建结果。图 18给出了重建弹道与外测弹道对比，可见二者一致性很好。
|
图 18 MF-1重建弹道与飞行弹道对比 Fig.18 Reconstructed trajectory compared with flight trajectory |
MF-1滚转角速率测量数据的分辨率为±0.25°/s，可见弹道重建得到的滚转角速率漂移量完全在传感器分辨率误差之内。然而，与第二节MF-1飞行弹道偏差影响仿真中，将引起发射初期弹道角差异归结为发射角偏差而不是风场偏差的仿真结果相比，虽然本节弹道重建得到的发射角偏差量值明显减小，更趋合理，但仍然比发射架技术指标给出的要大，方位角偏差大了一个量级。从前面的图 17可见，在弹道初期，直接积分弹道就与实际飞行弹道的星下点在侧向出现了明显偏离。本节前面已说过，气动力、推力、风场等偏差不会对弹道重建产生影响，因此，导致在弹道初期直接积分弹道就与实际飞行弹道的星下点在侧向出现了明显偏离的最主要原因，只能是发射初始方位角出现了较大偏差，因此，本文弹道重建结果给出的发射初始方位角偏差的量级是正确的。也就是说，MF-1发射时发射架方位角偏差确实超出了技术指标给出的精度范围。
基于弹道重建结果和风场数据，即可估计得到飞行迎角和侧滑角历程。图 19给出了上升段飞行时间0~40 s、高度40 km以下未进行风修正的迎角侧滑角历程与利用16：00气象测量结果进行风修正后的迎角侧滑角历程对比，可见风的影响很大，必须予以考虑。
|
图 19 风修正前后MF-1飞行迎角及侧滑角对比 Fig.19 Flight angles of attack and slid comparison between pre- and after-wind correcting |
本文首先利用弹道仿真手段研究了导致MF-1试验飞行器飞行弹道与设计弹道之间偏差的可能原因；然后，通过弹道重建技术开展了MF-1飞行试验全弹道重建，得到了完整可靠的全弹道数据；将弹道重建结果与前面的弹道偏差影响仿真结果相结合，进一步分析了导致飞行弹道差异的原因。结果表明：
1) 引起MF-1飞行弹道与设计弹道偏差的最主要因素应该是风场测量偏差，以及一定程度的气动力偏差和发射方位角偏差；
2) 利用弹道重建技术，获得了包括迎角侧滑角在内的完整可靠的弹道数据，为MF-1飞行试验结果分析提供了依据；
3) 风修正前后的迎角侧滑角对比表明，风场对迎角侧滑角的准确估计有显著影响，因此，在无法直接实现飞行迎角测量情况下，应尽量提高风场测量的实时性。
| [1] |
张守言. 模型自由飞试验[M]. 北京: 国防工业出版社, 2002. ZHANG S Y. Modeling free flight[M]. Beijing: National Defense Industry Press, 2002. (in Chinese) |
| [2] |
战培国. 美国空气动力学飞行试验平台综述[J]. 飞行力学, 2015, 33(5): 385-389. ZHAN P G. An overview of American aerodynamic flight testbeds[J]. Flight Dynamics, 2015, 33(5): 385-389. (in Chinese) |
| [3] |
叶蕾. 美国高超声速计划发展规律探索[J]. 飞航导弹, 2008(12): 20-24. YE L. A review of American aerodynamic flight test platform[J]. Aerodynamic Missile Journal, 2008(12): 20-24. (in Chinese) |
| [4] |
和争春, 车竞, 肖涵山, 等. MF-1弹道设计与蒙特卡罗飞行仿真[J]. 空气动力学学报, 2019, 37(02): 201-206. HE Z C, CHE J, XIAO H S, et al. Trajectory design and Monte Carlo flight simulation for MF-1[J]. Acta Aerodynamica Sinica, 2019, 37(2): 201-206. DOI:10.7638/kqdlxxb-2018.0062 (in Chinese) |
| [5] |
欧朝, 吉洪亮, 肖涵山, 等. MF-1模型飞行试验结构与热防护关键问题研究[J]. 空气动力学学报, 2017, 35(5): 742-749. OU C, JI H L, XIAO H S, et al. Key problems in structure and thermal protection for MF-1 model testing flight vehicle[J]. Acta Aerodynamica Sinica, 2017, 35(5): 742-749. DOI:10.7638/kqdlxxb-2017.0170 (in Chinese) |
| [6] |
和争春.飞船返回舱气动参数辨识及其有关问题研究[D].西安: 西北工业大学, 2004. HE Z C. Research on aerodynamic parameters identification and related problems for spaceship re-entry capsule[D]. Xi'an: Northwestern Polytechnical University, 2004. (in Chinese) |
| [7] |
汪清, 和争春, 方方, 等. 返回舱弹道重建与黑障区弹道再现研究[J]. 宇航学报, 2004(6): 595-599. WANG Q, HE Z C, FANG F, et al. Trajectory reconstruction and black-out-area trajectory makeup for reentry capsules of spacecraft[J]. Journal of Astronautics, 2004(6): 595-599. DOI:10.3321/j.issn:1000-1328.2004.06.001 (in Chinese) |
| [8] |
MULDER J A, CHU Q P, SRIDHAR J K, et al. Non-linear aircraft flight path reconstruction review and new advances[J]. Progress in Aerospace Science, 1999, 35: 673-726. DOI:10.1016/S0376-0421(99)00005-6 |
| [9] |
GERLACH O H. Analysis of apossible method for the measurement of performance and stability and control characteristics in non-steady symmetrical flight[R]. Report VTH-117, Delft University of Technology, Department of Aerospace engineering, 1963.
|
| [10] |
JONKERS H L. Application of Kalman filter to flight path reconstruction from test data including estimation of instrumental bias error corrections[R]. Report VTH-162, Delft University of Technology, Department of Aerospace engineering, 1976.
|
| [11] |
蔡金狮. 飞行器系统辨识学[M]. 北京: 国防工业出版社, 2003. CAI J S. System identification of flight vehicles[M]. Beijing: National Defense Industry Press, 2003. (in Chinese) |