FAST反射面形狀調(diào)節(jié)策略研究

牟淼 孔婭妮 郭磊

摘? 要：文章綜合考慮FAST上各個(gè)主索節(jié)點(diǎn)在經(jīng)度方向、緯度方向上的位移變化和邊緣光滑度等因素，研究了反射面的變形策略，借助粒子群算法，建立了較為符合實(shí)際的理想拋物面模型;然后利用歐拉旋轉(zhuǎn)定理建立旋轉(zhuǎn)模型，得到旋轉(zhuǎn)后的理想拋物面模型，再通過(guò)旋轉(zhuǎn)基準(zhǔn)球面固定照明區(qū)域的方法，求得工作拋物面上照明區(qū)域內(nèi)692個(gè)主索節(jié)點(diǎn)編號(hào)，伸縮量以及位置坐標(biāo);最后利用降維的技巧，求出基準(zhǔn)球面和工作拋物面上的反射信號(hào)接收比。

關(guān)鍵詞：變形策略;理想拋物面;旋轉(zhuǎn)模型;接收比;粒子群算法;降維

中圖分類(lèi)號(hào)：TP181;P111.44? ? ? ? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：2096-4706（2021）16-0086-07

Research on Shape Adjustment Strategy of FAST Reflector

MOU Miao1， KONG Yani2， GUO Lei3

（1.School of Mechanical Engineering， Yangtze University， Jingzhou? 434023， China; 2.School of Information and Mathematics， Yangtze University， Jingzhou? 434023， China; 3.School of Physics and Optoelectronic Engineering， Yangtze University， Jingzhou? 434023， China）

Abstract： In this paper， considering the displacement changes and edge smoothness of each main cable node in longitude and latitude direction on FAST， the deformation strategy of reflector is studied， and an more realistic ideal paraboloid model is established with the help of particle swarm optimization algorithm; then， the rotation model is established by using Euler rotation theorem， and the ideal paraboloid model after rotation is obtained. Then， by rotating the reference sphere to fix the lighting area， the number， expansion and position coordinates of 692 main cable nodes in the lighting area on the working paraboloid are obtained; finally， using the technique of dimension reduction， the reflection signal reception ratio on the reference sphere and the working paraboloid is obtained.

Keywords： deformation strategy; ideal paraboloid; rotation model; receiving ratio; particle swarm optimization; dimension reduction

0? 引? 言

中國(guó)天眼FAST是迄今為止世界上最先進(jìn)的空間探測(cè)設(shè)備，可以接收到宇137億光年外的電磁信號(hào)，由饋源支撐塔、主動(dòng)反射面和饋源倉(cāng)組成。本文將對(duì)2021年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽題進(jìn)一步研究，建立了合理的評(píng)估指標(biāo)，利用“粒子群”搜索算法，尋找最為符合實(shí)際的理想拋物面模型，然后利用歐拉旋轉(zhuǎn)定理旋轉(zhuǎn)拋物面方程得到第一種旋轉(zhuǎn)模型，旋轉(zhuǎn)后的理想拋物面方程;再通過(guò)旋轉(zhuǎn)基準(zhǔn)球面得到第二種旋轉(zhuǎn)模型，固定照明區(qū)域，逆向旋轉(zhuǎn)基準(zhǔn)球面，得到工作拋物面上照明區(qū)域內(nèi)692個(gè)主索點(diǎn)編號(hào)，然后通過(guò)計(jì)算得到伸縮量以及位置坐標(biāo)，最終利用“降維”的技巧，求得基準(zhǔn)球面和工作拋物面上的反射信號(hào)比。

1? 天眼FAST基本介紹

1.1? FAST的組成

FAST主要由主動(dòng)反射面，信號(hào)接收系統(tǒng)以及相關(guān)的控制、測(cè)量和支承系統(tǒng)組成，如圖1所示。其中主動(dòng)反射面是由主索網(wǎng)、反射面板、下拉索、促動(dòng)器及支承結(jié)構(gòu)等主要部件構(gòu)成的一個(gè)可調(diào)節(jié)球面。主索面上每一個(gè)主索點(diǎn)處下端都通過(guò)下拉索與促動(dòng)器連接，通過(guò)調(diào)節(jié)各個(gè)主索點(diǎn)對(duì)應(yīng)的促動(dòng)器的伸縮量進(jìn)而實(shí)現(xiàn)反射面形狀的調(diào)節(jié)。

1.2? 主動(dòng)反射面的基準(zhǔn)態(tài)和工作態(tài)

主動(dòng)反射面可以分為基準(zhǔn)態(tài)和工作態(tài)兩種狀態(tài)，其中基準(zhǔn)態(tài)時(shí)反射面可以看作是半徑為300 m，口徑為500 m的球面，而當(dāng)主動(dòng)反射面處于工作態(tài)時(shí)，反射面調(diào)節(jié)成口徑為300 m的反射拋物面，圖2是FAST處于工作態(tài)時(shí)剖面示意圖，兩個(gè)同心球面的半徑差為F=0.466R，其中C點(diǎn)為基準(zhǔn)球面球心，饋源艙接受平面的中心只能位于與基準(zhǔn)球面同心的球面上，只有當(dāng)反射信號(hào)經(jīng)過(guò)以饋源艙中心為圓心，直徑為1 m的圓盤(pán)時(shí)，才算接收到信號(hào)。

2? 構(gòu)建模型思路

本文的最終目的是求解基準(zhǔn)球面和工作拋物面反射出信號(hào)的有效接收比，而為了求解工作拋物面的接收比，首先應(yīng)求解工作拋物面上照明區(qū)域各個(gè)主索點(diǎn)對(duì)應(yīng)促動(dòng)器的伸縮量，又因?yàn)樽顑?yōu)工作拋物面是在約束條件的情況下最接近理想拋物面的拋物面，所以要優(yōu)先計(jì)算理想拋物面方程，具體思維結(jié)構(gòu)如圖3所示。為了得到比較符合實(shí)際的理想拋物面方程，本文從徑向，緯向的位移變化和邊緣光滑度等因素出發(fā)，建立符合實(shí)際的評(píng)價(jià)準(zhǔn)則，利用“粒子群”優(yōu)化算法擬合最優(yōu)參數(shù)，最終求得結(jié)果。

3? 確定理想拋物面方程

3.1? 綜合評(píng)價(jià)準(zhǔn)則的確定

本文根據(jù)賽題提供的數(shù)據(jù)，進(jìn)行分析處理后，工作區(qū)域?yàn)榭趶?00米及以下的主索節(jié)點(diǎn)。為了更加形象地表示出基準(zhǔn)球面，繪出基準(zhǔn)面外形圖如圖4所示。

由于優(yōu)化模型的目標(biāo)是使主索節(jié)點(diǎn)沿徑向、經(jīng)向、緯向以及在拋物面邊緣部分最光滑，故可取權(quán)重：。

得到目標(biāo)函數(shù)：

N=aL1+bL2+cL3+dL4-ey’? ? ? ? ? ? ? ? ? ? （1）

式中：L1代表徑向位移;L2代表經(jīng)向位移;L3代表緯向位移;L4代表邊緣徑向位移;y’代表邊緣斜率。

3.2? 模型降維

主動(dòng)反射面在工作時(shí)由球面變?yōu)樗矔r(shí)拋物面，節(jié)點(diǎn)在徑向、經(jīng)向和緯度方向上產(chǎn)生位移，根據(jù)球面和拋物面具有對(duì)稱性可將此三維模型簡(jiǎn)化為二維模型，球面半徑與拋物線口徑相等R=D=300 m。以球面圓心為坐標(biāo)原點(diǎn)可得圓弧和拋物線方程，其中f為焦徑比：

（2）

（3）

3.3? 理想拋物面模型確定

3.3.1? 徑向位移

球面半徑減去節(jié)點(diǎn)到球心的距離即節(jié)點(diǎn)沿徑向方向上的位移，由此可得：

L1=R-L? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （4）

L為節(jié)點(diǎn)到球心的距離，L1（x，f，h）其中x范圍在[0，150]，h變化幅度為0.6。

3.3.2? 經(jīng)向位移

基準(zhǔn)球面變形為拋物面上節(jié)點(diǎn)到頂點(diǎn)的弧長(zhǎng)減去原基準(zhǔn)球面節(jié)點(diǎn)到頂點(diǎn)的弧長(zhǎng)即為經(jīng)向位移?；鶞?zhǔn)球面變?yōu)樗矔r(shí)拋物面時(shí)300米口徑以內(nèi)的弧長(zhǎng)會(huì)發(fā)生變化，即產(chǎn)生經(jīng)向位移，對(duì)經(jīng)向位移取微元

則有：

（5）

（6）

（7）

基準(zhǔn)球面變形為拋物面后拋物線上主索節(jié)點(diǎn)M到頂點(diǎn)弧長(zhǎng)為：

（8）

基準(zhǔn)球面的主索節(jié)點(diǎn)M’到球面頂點(diǎn)弧長(zhǎng)為：

經(jīng)向位移：

L2=sM-sM’? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（9）

3.3.3? 徑向位移

假定主索節(jié)點(diǎn)N在圓弧上的橫坐標(biāo)為x0，在拋物線上的橫坐標(biāo)x’，緯度環(huán)環(huán)長(zhǎng)變化量為l，相對(duì)基準(zhǔn)球面緯度環(huán)環(huán)長(zhǎng)變化量為：

（10）

3.3.4? 徑向位移

拋物面變形邊緣應(yīng)盡量光滑，拋物線邊緣處斜率與圓弧斜率應(yīng)盡量接近，拋物線邊緣徑向移動(dòng)位移與圓弧徑向移動(dòng)位移之差盡量小。

基準(zhǔn)球面與拋物面連接處x0=150，y’0=0.577 4，拋物面的導(dǎo)數(shù)：

（11）

邊緣徑向位移：L4=L1（149.146 2，f，h）

故該題建立的“理想曲面優(yōu)化模型”具體為：

目標(biāo)函數(shù)：MinN=aL1+bL2+cL3+dL4-ey’

3.4? 理想拋物面模型求解

3.4.1? 搜索算法的構(gòu)建

該問(wèn)題本質(zhì)上是多目標(biāo)優(yōu)化模型，通過(guò)“粒子群”智能優(yōu)化算法（圖5）的過(guò)程如下：

（1）確定各個(gè)參數(shù)的搜索區(qū)間[v1，v2];

（2）賦給每個(gè)指標(biāo)一定的權(quán)重，確定綜合評(píng)價(jià)標(biāo)準(zhǔn)為目標(biāo)函數(shù);

（3）初始化粒子位置，并計(jì)算此時(shí)綜合評(píng)價(jià)指標(biāo);

（4）重復(fù)步驟（3），直至目標(biāo)函數(shù)值逐步收斂于最優(yōu)目標(biāo)函數(shù)值;

（5）輸出最優(yōu)目標(biāo)函數(shù)對(duì)應(yīng)的各個(gè)最優(yōu)參數(shù)。

3.4.2? 理想拋物面模型的求解結(jié)果

將上述粒子群算法過(guò)程通過(guò)matlab迭代求解得出結(jié)果如圖6所示，兩個(gè)參數(shù)如下：

f=0.467 7，h=0.412 2? ? ? ? ? ? ? ? ? ? ? ? ? ? （12）

式中：f為焦徑比;h為頂點(diǎn)移動(dòng)距離。

故理想拋物面方程為：

理想拋物面與基準(zhǔn)球面的擬合曲線如圖7和圖8所示。

4? 確定工作拋物面方程和主索點(diǎn)的伸縮量

4.1? 旋轉(zhuǎn)后工作拋物面模型建立

首先選取在α=0°，β=90°時(shí)確定出的理想拋物面：

（13）

將拋物面進(jìn)行繞坐標(biāo)軸旋轉(zhuǎn)即在三維空間下曲面繞x軸旋轉(zhuǎn)一定角度，再繞y軸旋轉(zhuǎn)至對(duì)稱軸在α=36.795°，β=78.169°方向上，需要計(jì)算出繞軸旋轉(zhuǎn)的角度。以基準(zhǔn)球面球心為坐標(biāo)原點(diǎn)理想拋物面對(duì)稱軸為z軸建立直角坐標(biāo)系，將天體看作一點(diǎn)投射至坐標(biāo)面上連結(jié)各點(diǎn)，由此便得到了第一種旋轉(zhuǎn)方法，如圖9所示。

其中α為方位角，β為仰角，繞x軸旋轉(zhuǎn)的角度為∠mcp，繞y軸旋轉(zhuǎn)的角度為∠mcs。假設(shè)sc=a，則：

（14）

（15）

（16）

令，α和β已知，即確定了沿x軸，y軸的旋轉(zhuǎn)角度，則可利用旋轉(zhuǎn)矩陣對(duì)問(wèn)題一中確定的理想拋物面進(jìn)行旋轉(zhuǎn)：

假設(shè)原拋物面上一點(diǎn)（x，y，z），首先對(duì)x軸進(jìn)行旋轉(zhuǎn)得到旋轉(zhuǎn)后坐標(biāo)（x1，y1，z1）：

（17）

然后繞y軸進(jìn)行旋轉(zhuǎn)得到旋轉(zhuǎn)后坐標(biāo)（x2，y2，z2）：

（18）

（19）

記：

故可以得到旋轉(zhuǎn)后的理想拋物面，即在α=36.795°，β=78.169°時(shí)的理想拋物面方程：

（20）

4.2? 工作拋物面上主索點(diǎn)伸縮量模型建立

對(duì)基準(zhǔn)球面上所有的主索節(jié)點(diǎn)進(jìn)行旋轉(zhuǎn)，第一次旋轉(zhuǎn)繞z軸旋轉(zhuǎn)α角度，第二次旋轉(zhuǎn)繞y軸旋轉(zhuǎn)角度，即可完成旋轉(zhuǎn)，第二種旋轉(zhuǎn)模型如圖10所示。

圖10? 第二種旋轉(zhuǎn)模型

繞z軸旋轉(zhuǎn)矩陣為：

（21）

繞y軸旋轉(zhuǎn)矩陣為：

（22）

其中：α=36.795°，β=78.169°

選取其中在工作面范圍（照明區(qū)域）內(nèi)的主索節(jié)點(diǎn)，即x2+y2≤1 502。

根據(jù)這些點(diǎn)確定其位置坐標(biāo)，計(jì)算各個(gè)主索節(jié)點(diǎn)與球心連線之間的直線方程，求出直線與理想拋物線的交點(diǎn)，即理想主索的位置坐標(biāo)：

（23）

解出理想拋物面上主索節(jié)點(diǎn)的位置坐標(biāo)，且已確定基準(zhǔn)球面旋轉(zhuǎn)后的主索節(jié)點(diǎn)，即可利用兩點(diǎn)間的距離公式計(jì)算理想化的伸縮量如圖11所示。

（24）

其中Ll為主索節(jié)點(diǎn)理想的伸縮量，即不受任何約束情況下的伸縮量。

接下來(lái)考慮約束條件，主要為：

（25）

其中Ls為實(shí)際的伸縮量，xi為第i個(gè)主索節(jié)點(diǎn)的原長(zhǎng)度

根據(jù)以上約束條件，即可計(jì)算出每一個(gè)在工作區(qū)域內(nèi)的主索節(jié)點(diǎn)的伸縮量以及相應(yīng)的位置坐標(biāo)：

4.3? 工作拋物面模型的求解

利用第一種旋轉(zhuǎn)方式旋轉(zhuǎn)所得理想拋物面如圖12所示，其求解理想拋物面：

α=36.795°，β=78.169°

定義隨機(jī)變量x，y，z

N=-arctan（cos（α）cos（β））

a=zeros（3）;b=zeros（3）

帶入原理想方程可得旋轉(zhuǎn)后的理想方程：

b（3）=1/561.24×（b（1）2+b（2）2）-300.212 2

根據(jù)以上步驟，即可得出處在α=36.795°，β=78.169°時(shí)的理想拋物線方程：

4.4? 工作拋物面上主索點(diǎn)伸縮量模型求解

4.4.1? 利用第二種旋轉(zhuǎn)方式求解頂點(diǎn)坐標(biāo)

將基準(zhǔn)面的頂點(diǎn)A0（0，0，-300.212 2）進(jìn)行旋轉(zhuǎn)，令：

即得頂點(diǎn)坐標(biāo)G（-49.289，-38.866，-293.835）

4.4.2? 求解各工作區(qū)域主索節(jié)點(diǎn)的編號(hào)及對(duì)應(yīng)的位置坐標(biāo)和伸縮量

求解旋轉(zhuǎn)后主索節(jié)點(diǎn)的編號(hào)，可將整個(gè)基準(zhǔn)球面旋轉(zhuǎn)α=36.795°，β=78.169°，將工作區(qū)域移到球心的正下方，旋轉(zhuǎn)后所有主索節(jié)點(diǎn)的位置坐標(biāo)滿足：x2+y2≤1502，將滿足該條件的主索節(jié)點(diǎn)與球心相連，與第一問(wèn)中理想拋物面的交點(diǎn)即為理想位置坐標(biāo)。相鄰節(jié)點(diǎn)距離變化幅度需要不超過(guò)0.07%，節(jié)點(diǎn)徑向移動(dòng)距離要在0.06以內(nèi)，故計(jì)算基準(zhǔn)球面上相應(yīng)主索節(jié)點(diǎn)坐標(biāo)和理想位置坐標(biāo)的距離，根據(jù)約束條件進(jìn)行約束，可導(dǎo)出相應(yīng)伸縮量如表1所示。

5? 確定基準(zhǔn)球面和工作面的反射信號(hào)接收比

饋源倉(cāng)是一個(gè)面積為0.785平方米的圓盤(pán)，需要計(jì)算出調(diào)節(jié)后饋源倉(cāng)的接收比和基準(zhǔn)反射球面的接收比，并將二者進(jìn)行比較。本文先將基準(zhǔn)反射球面進(jìn)行旋轉(zhuǎn)，讓工作區(qū)域恰好與豎直照明區(qū)域一致，基準(zhǔn)球面第一次旋轉(zhuǎn)繞z軸旋轉(zhuǎn)α角度，第二次旋轉(zhuǎn)繞y軸旋轉(zhuǎn)角度。

5.1? 建立基準(zhǔn)反射球面的反射信號(hào)接收比模型

由于基準(zhǔn)反射球面具有各向同性，故可將基準(zhǔn)反射球面簡(jiǎn)化為二維曲線。又因?yàn)閳A弧關(guān)于z軸對(duì)陣，故只需考慮x∈（0，150），如圖13所示。

在照明區(qū)域圓弧上，從0°至150°，取x步長(zhǎng)為0.1的1 500個(gè)點(diǎn)，然后根據(jù)圓弧的方程曲線：

x2+y2=3002? ? ? ? ? ? ? ? ? ? ? ? ? ? （26）

算出對(duì)應(yīng)點(diǎn)的y值，即可得到每個(gè)點(diǎn)的坐標(biāo)。圖13中

饋源倉(cāng)距離球心（1-0.466）R，故：

L=[y-（1-0.466）R]·tan（2θ）? ? ? ? ? ? ? ? ? （27）

若|L|≥0.5則不能有效接收，反之則能有效接收。

計(jì)算1 500個(gè)點(diǎn)中，能被接收到的點(diǎn)有多少個(gè)，即可計(jì)算出接收率。

5.2? 建立工作拋物面上反射信號(hào)接收比模型

根據(jù)第二問(wèn)求出的頂點(diǎn)坐標(biāo)，計(jì)算出對(duì)稱軸的直線方程。

天體發(fā)出的電磁波視為直線傳播，平行對(duì)稱軸射入，設(shè)入射電磁波方程為：

ax+by+cz+di=0? ? ? ? ? ? ? ? ? ? ? ? ?（28）

式中：a，b，c均已知，與對(duì)稱軸方程的系數(shù)相同。

根據(jù)每個(gè)主索節(jié)點(diǎn)的坐標(biāo)，即可計(jì)算出每個(gè)主索節(jié)點(diǎn)對(duì)應(yīng)的di。

主索節(jié)點(diǎn)與球心的連線方程為：

（29）

式中：xi，yi，zi為各個(gè)主索節(jié)點(diǎn)的坐標(biāo)。

將（28）繞（29）旋轉(zhuǎn)180度即可得到新的直線方程：

Ax+By+Cz+D=0? ? ? ? ? ? ? ? ? ? ? ? ?（30）

又因?yàn)轲佋磁摰挠行Ы邮諈^(qū)域容易取得，即為滿足如下條件：

（31）

本文即可認(rèn)為此主索節(jié)點(diǎn)上的信號(hào)可以有效接收到。

記共有n個(gè)主索節(jié)點(diǎn)在工作區(qū)域，有m個(gè)主索節(jié)點(diǎn)可以被有效接收到，所以可得接受率為：

（32）

式中：ξ為有效接收率。

5.3? 求解基準(zhǔn)球面和工作拋物面上反射信號(hào)接收比模型

根據(jù)物理模型球面反射，由于天體與地球的距離足夠可以將天體發(fā)出的電磁波看作平行波直線射入，需要得到最優(yōu)的信號(hào)之比。饋源倉(cāng)可視為一個(gè)小圓盤(pán)，接收反射的電磁波，由于實(shí)際情況下會(huì)有部分射入部分落在了范圍之外，射在曲面上的坐標(biāo)與反射到達(dá)焦面的坐標(biāo)是非線性關(guān)系，故本文可對(duì)其進(jìn)行機(jī)器學(xué)習(xí)的數(shù)據(jù)處理，使其在更高維度呈線性關(guān)系。

本文采用多項(xiàng)式核函數(shù)：

K（x，xi）=[p（xt·x）+s]q? ? ? ? ? ? ? ? ? （33）

將輸入集射入電磁波與曲面的交點(diǎn)坐標(biāo)F（f1，f2，f3…fm）通過(guò)選取的多項(xiàng)式核函數(shù)映射到高維空間得到X（x1，x2，x3…xn），本文將輸出的表達(dá)式表示為：

（34）

其中θ為本文要求的參數(shù)權(quán)值，根據(jù)損失函數(shù)計(jì)算公式：

（35）

結(jié)合選取的多項(xiàng)式核函數(shù)可以得到兩種損失函數(shù)：cost0，cost1，如圖14和圖15所示。

對(duì)于基準(zhǔn)球面接收比的求解為：8.51%;對(duì)于調(diào)整后拋物面接收比的求解為：83.20%。

6? 結(jié)? 論

本文針對(duì)FAST反射面形狀調(diào)節(jié)進(jìn)行了研究，首先通過(guò)機(jī)理建模的方式，從徑向，緯向，經(jīng)向方向上位移的變化以及邊緣光滑度等因素出發(fā)，建立了綜合評(píng)價(jià)準(zhǔn)則，利用“粒子群”智能算法，搜索出符合各個(gè)因素的最優(yōu)參數(shù)，從而得到符合實(shí)際的理想拋物面方程。在理想拋物面確定的情況下，利用兩種旋轉(zhuǎn)方式分別求出旋轉(zhuǎn)后理想拋物面方程以及旋轉(zhuǎn)后照明區(qū)域內(nèi)主索點(diǎn)的伸縮變化量和理想拋物面頂點(diǎn)坐標(biāo)，從而得到在符合約束條件的工作拋物面模型，最后再分別計(jì)算工作拋物面模型和基準(zhǔn)球面模型的反射信號(hào)接收比。

在實(shí)際中，本文建出拋物面最佳變形策略模型，對(duì)于家用電視天線的形狀設(shè)計(jì)，聲、波接收面的優(yōu)化等具有較強(qiáng)的指導(dǎo)意義。同時(shí)給FAST天眼提供了較好的變形策略，有利于更好的接收信號(hào)。

參考文獻(xiàn)：

[1] 李明輝，朱麗春.FAST瞬時(shí)拋物面變形策略優(yōu)化分析 [J].貴州大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，29（6）：24-28+43.

[2]朱麗春.500米口徑球面射電望遠(yuǎn)鏡（FAST）主動(dòng)反射面整網(wǎng)變形控制 [J].科研信息化技術(shù)與應(yīng)用，2012，3（4）：67-75.

[3] 王璐.大型陣面位姿調(diào)整系統(tǒng)的布局優(yōu)化和協(xié)調(diào)控制研究 [D].哈爾濱：哈爾濱工程大學(xué)，2019.

[4] 卓金武，李必文，魏永生，等.MATLAB在數(shù)學(xué)建模中的應(yīng)用：第2版 [M].北京：北京航空航天大學(xué)出版社，2013.

[5] 天工在線.MATLAB從入門(mén)到精通 [M].北京：中國(guó)水利水電出版社，2018.

作者簡(jiǎn)介：牟淼（2000.09—），男，漢族，山東濟(jì)南人，本科在讀，研究方向：機(jī)械工程;孔婭妮（2001.11—），女，漢族，湖北咸寧人，本科在讀，研究方向：數(shù)學(xué)與應(yīng)用數(shù)學(xué);郭磊（2000.08—），男，漢族，山東濟(jì)南人，本科在讀，研究方向：應(yīng)用物理。